How can the zero-lift drag coefficient (parasitic drag) be calculated? Consider a 3-D wing made from an arbitrary airfoil, say a NACA0012 airfoil. The wing has a trapezoidal shape, with a fixed span, root chord, and tip chord. Also, assume that the wing loading is known as well. I am trying to calculate the speed at minimum drag of this wing (assume that there are no other parts of the aircraft, just the wing!) My thought process is as follows:

We know that, to a reasonable degree of accuracy, that there are two types of drag on the wing in steady, level, flight: parasitic drag and lift-induced drag. This can be shown mathematically as: 

$$ C_D = C_{D_0} + C_{D_i} = C_{D_0} + \frac{C_L^2}{\pi e AR}$$ 

Also, assume that the AR and efficiency factor are known. Now, for minimum drag to occur, there has to me a maximum lift to drag ratio. The formula for drag is 

$$ D = \frac{1}{2} \rho V^2 S C_D = \frac{1}{2} \rho V^2 S \Big(C_{D_0} + \frac{C_L^2}{\pi e AR}\Big) = \frac{1}{2} \rho V^2 S C_{D_0} + \frac{\rho V^2 S}{2\pi e AR} C_L^2$$

Lift has a similar formula to drag, and in steady, level flight, is equal to the weight of the aircraft. Lift is related to the lift coefficient as $L = \frac{1}{2} \rho V^2 S C_L$. So we solve for the lift coefficient as follows.

$$ C_L = \frac{2L}{\rho V^2 S} = \frac{2W}{\rho V^2 S}$$. 

Plugging into our original formula, we obtain

$$ D = \frac{1}{2} \rho V^2 S C_{D_0} + \frac{\rho V^2 S}{2\pi e AR} \cdot \frac{4W^2}{\rho^2 V^4 S^2} = \frac{1}{2} \rho S C_{D_0} V^2  + \frac{2W^2}{\pi e AR \rho S}\frac{1}{V^2}$$ 

This is great for us, because now we have a relationship between drag and lift, and to find the speed at minimum drag, all we have to do is take the derivative and set it equal to 0. I've done this and the resulting answer comes out to be 

$$V_{md} = \Bigg( \frac{4W^2}{\rho^2 S^2 \pi e AR C_{D_0}} \Bigg)^{1/4},$$
where the 'md' stands for minimum drag. My problem arises because I cannot for the life of me figure out how to analytically calculate $C_{D_0}$. It can also be shown that at minimum drag, $C_{D_0} = C_{D_i}$ so that the total drag coefficient becomes $C_D \equiv C_{D_0} + C_{D_i} = 2C_{D_i} = \frac{2C_L^2}{\pi e AR}$, but then we are back to our starting point, which confuses me again. 

My last resort was to read some papers that said that there is a method to finding $C_{D_0}$ using the skin friction coefficient, because at subsonic speeds, a large part of the parasitic drag is due to skin friction (and a little due to pressure drag). Anyways, this led me to the formula $C_{D_0} = C_{fe}\frac{S_{wetted}}{S_{ref}},$ where you use an equivalent skin friction and wetted area. Now I don't understand what a wetted surface area is, since in this example we are only dealing with one wing (would it just be twice the regular area??) As you can see, I am greatly confused. How do you find this zero-lift drag, and subsequently the minimum flight speed.
Yes, the wetted area is roughly twice the reference area. Now details depend on how well the reference area captures the exposed area of the wing - dihedral will already increase the wetted area by a factor proportional to the inverse of the cosine of the dihedral angle.

But there is more. Airfoil thickness means that the air has to flow around the airfoil. This displacement effect causes the flow around a thick airfoil to speed up more than around an equivalent but thinner airfoil. The thicker airfoil pushes the air aside and around itself more, causing the flow to accelerate and create more friction than the slower flow around a thinner airfoil. This effect is normally approximated with an additional term in the friction drag formula which is proportional to relative thickness.

Next, the type of boundary layer flow needs to be known. Rough surfaces or high sweep angles will provoke an early transition from laminar to turbulent flow. Read this answer for a more detailed discussion.

Another correction is needed for the Mach number, even in subsonic flow. Of course, once flow becomes transsonic or supersonic, wave drag needs to be added, too.

First you need to calculate the friction coefficient which depends on the Reynolds and Mach numbers of your airfoil flow and the relative mean roughness R:
$$c_f = \frac{\frac{0.43}{log(100/R)^{2.56}}-\frac{1700}{100/R}}{\sqrt{1+0.14\cdot Ma^2}}$$

Next, you approximate airfoil drag as explained above:
$$c_{d0} = c_f\cdot \left(2 + 4\cdot\delta + 120\cdot\left(\frac{1}{\sqrt{1-Ma^2}}\right)^3\cdot\delta^4 - 0.09\cdot Ma^2\right)$$
where $\delta$ is the relative thickness of your airfoil.

The term $\frac{1700}{100/R}$ in the friction drag equation allows for the initially laminar boundary layer. Change the factor 1700 depending on how much laminarity your airfoil offers. This answer shows a graph with the possible range. In the airfoil zero lift drag formula you see first the factor 2 which accounts for the fact that the wing has two sides. To that you add the thickness summand to allow for the displacement effect. The third term with the Prandtl-Glauert factor shows that the formula only works well for Mach &lt; 1, and both the third and fourth term are empirical factors to improve accuracy over Mach.