2D BOUNDARY LAYER MODELLING

The analysis of an aerofoil's boundary layer can be done as an extension of the results for a simple flat plate. Effects due to surface curvature and angular acceleration around leading edge can be lumped into a single parameter, the input definition of pressure gradient distribution along the surface. The geometric model is simplified to a flat surface with the pressure gradient and velocity distribution that has been obtained for section inviscid flow analysis.

subsonic5_html_2f7b2b2e

The height of the boundary layer δ can be estimated as the point where u = 0.99 x Uδ. Here Uδ is the velocity of the inviscid flow field outside the viscous boundary layer. In this analysis Uδ is obtained from the previous potential flow panel method solution as the tangential surface velocities on the surface of each of the panels.

partial P / partial x is the pressure gradient along the aerofoil surface (assuming static pressure P is constant through the depth of the boundary layer).

tau_w is the shear stress on the surface due to friction with the air flow.

P = P_o - {1 over 2} rho U_delta^2 for the inviscid external flow where Po is the stagnation pressure of the flow.

For incompressible flows the stagnation pressure is constant so that

{partial P over partial x} = - rho U_delta {partial U_delta over partial x}

The governing equations for this simple one-dimensional flow model can be found by applying the laws of conservation of mass and momentum to a control volume covering the boundary layer.

subsonic5_html_bl1
Conservation of mass : Sigma dot m = 0

Conservation of momentum : Sigma dot m u = Sigma F_x

Mass and Momentum flow will occur at the three free surfaces (1) , (2) and (3) of the control volume.

dot m_1 = int_0^delta rho u dydot m_2 = - int_0^delta ( rho u + { partial rho u over partial x} .dx ) dydot m_3 = - rho v_delta .dx Thus

int_0^delta {partial rho u over partial x}}.dx.dy = - rho v_delta dx

The momentum balance will be,

Sigma dot m u = -( dot m u)_1 +( dot m u)_2 + ( dot m u)_3

Sigma dot m u = int_0^delta ( rho u^2 + {partial rho u^2 over partial x}.dx ).dy - int_0^delta rho u^2.dy + rho v_delta U_delta dx
The horizontal forces on the control volume will come from pressure on side 1 and side 2 and shear stress from the surface.

Sigma F_x = P dot delta - (P + {partial P over partial x}.dx)cdot delta - tau_w .dx

Equating momentum change to applied forces give

rho v_delta U_delta dx + int_0^delta {partial rho u^2 over partial x}.dx.dy = - {partial P over partial x}.delta.dx - tau_w .dx

= rho U_delta {partial U_delta over partial x}.delta.dx - tau_w .dx

Using the conservation of mass result to substitute for v_delta and cancelling element length dx gives,

- u_delta int_0^delta {partial rho u over partial x}.dy + int_0^delta {partial rho u^2 over partial x }.dy = rho U_delta {partial U_delta over partial x} delta - tau_w = rho U_delta {partial U_delta over partial x} int_0^delta 1.dy - tau_w
The velocity profile of the layer can be further processed by defining terms specific to its size, displacement effect and momentum loss.subsonic5_html_22de9660Boundary layer thickness, delta = int_0^delta dy

Displacement effect or lost mass flow, delta^* = int_0^delta (1 - { rho u over rho_delta U_delta}) dy

Momentum loss in profile, theta = int_0^delta { rho u over rho_delta U_delta}(1 - { rho u over rho_delta U_delta}) dy

and Wall shear stress, tau_w = ( mu {partial u over partial x} )_{y to 0}

These definitions are applied to low speed incompressible flow where ρ is constant through the boundary layer. If the above momentum equation is divided by external momentum, rho u_delta^2 and the above profile definitions are substituted in, then the following form is obtained,

{partial theta over partial x} + {(2 theta + delta^*) over U_delta}.{partial U_delta over partial x}- {tau_w over rho U_delta} = 0
In practice a boundary layer shape factor, H = delta^* / theta , can be defined. Its values do not vary greatly and this can be used to simplify the solution methods.

{partial theta over partial x} + (2 + H) {theta over U_delta}.{partial U_delta over partial x}- {tau_w over rho U_delta} = 0
The solution of this approximately linear differential equation can be found by starting with initial conditions at the stagnation point and then integrating along the surface in small steps, dx, to find the momentum loss at any point. To simplify calculations, curvature of the aerofoil surface is ignored and assumed polynomial shapes for the boundary layer profile are used to predict local expressions for δ, δ*, θ and τw.

Laminar boundary layer.

Initially the flow will be laminar. An exact solution with zero pressure gradient has been predicted by Blassius, so the expected shape for the aerofoil flow will be based on a polynomial fit to this solution, plus an addition term to handle the pressure gradient effect.

{u over U_delta} = 2 eta- 2 eta^2+ eta^4+{eta over 6}(1-eta)^3 Lambda
where eta = y /delta and Lambda = delta^2 frac{rho}{mu} frac{ partial U_delta}{partial x}
with Λ = 7.052 at the stagnation point, Λ > 0 for accelerating flow, Λ = 0 when there is no pressure gradient (Blassius solution), Λ < 0 for decelerating flow and Λ = -12 where the boundary layer flow starts to separate.

Boundary layer shapes will then appear as shown in the following figure and will depend on their location along the section and the local pressure gradient (velocity gradient).

subsonic5_html_fd5ac16

By integration and differentiation of the polynomials for this family of profiles, the boundary layer parameters can be found,

{delta^* over delta}=0.3-{Lambda over 120}
{theta over delta}={37 over 315}-{Lambda over 945}}- { Lambda^2 over 9072}
{tau_w over rho U_delta^2}={mu over rho U_delta delta} (2 + {Lambda over 6} )
Substituting these expressions into the boundary layer momentum equation and assembling terms with θ as the primary variable leads to,

{rho U_delta over mu} .{partial theta^2 over partial x}+ 6 {rho theta^2 over mu}.{partial U_delta over partial x} - 0.47 =0
multiplying the equation by U_delta^5 gives,

{rho over mu} .U_delta^6.{partial theta^2 over partial x}+ 6.U_delta^5.{rho theta^2 over mu}.{partial U_delta over partial x} = 0.47 U_delta^5
or

{rho over mu} .{partial theta^2 U_delta^6 over partial x} = 0.47 U_delta^5
which can be integrated with respect to x to obtain a solution along the boundary layer.

{rho over mu} int d(U_delta^6 theta^2) = 0.47 int U_delta^5 .dx
This integration can be performed from the stagnation point to any position on the aerofoil surface. The conditions at the stagnation point are Uδ=0 and θ=0, so the result of the integration becomes,

{rho over mu} U_delta^6 theta^2 = 0.47 int_0^x U_delta^5 .dx
or

theta^2 = 0.47 {mu over rho} {1 over U_delta^6 } int_0^x U_delta^5 .dx
and in non-dimensional form,

R_{e theta}^2 = R_e {0.47 over (U_delta / V_infty)^4} int_0^{x/c} ({U_delta over V_infty})^5 .d{x over c}
where Reynolds number based on momentum loss is

R_{e theta} = {rho theta U_delta over mu}
Applying the previous assumptions of profile shape, including some curve fitting of the variables gives the other boundary layer parameters at this point (x),

H = {delta^* over theta} = 2.61 - 3.75K + 5.24K^2 for 0<K<0.1
H = 2.088+ {0.731 over {K+0.14}} for -0.1<K<0
where

K = {rho over mu} theta^2 {partial U_delta over partial x} = ( {theta over delta} )^2.Lambda
In most cases the boundary layer does not remain laminar all the way to the trailing edge, so this solution only gives momentum loss up to the point where instability sets in and the boundary layer becomes turbulent.

Transition

At low angles of attack, transition is usually caused by the growth of natural instabilities at high Reynolds number and by the onset of an adverse pressure gradient.

An empirical prediction of this effect has been obtained by Michel, such that transition occurs when,

R_{e theta} > 1.174 ( 1+ {22400 over R_{ex}}).R_{ex}^{0.46}
At higher angles, with sudden adverse pressure gradients, the laminar boundary layer will start to separate causing a small bubble. The shear flow over this bubble becomes unstable leading to turbulence and a reattached boundary layer. This condition can be predicted by observing when the boundary layer reaches a value of Λ=-12.

However it is difficult to assess whether a laminar bubble transition occurs or whether the laminar flow completely separates and stall occurs. This prediction is getting very close to the limit of application of this simple theory and care should be taken in its use.

Turbulent Boundary Layers

Once the surface layer becomes unstable the flow breaks up into random turbulence. The profile at a particular location is no longer fixed but varies with time. In this case profiles can only be analysed as time averages with an associated level of turbulence at each height.

subsonic5_html_202a1bf4
Instantaneous velocities in turbulent layer at three separate times.
subsonic5_html_m6ec26d5d
Time averaged velocity profile also showing variation of turbulence level through the layer.

The time average profile is itself rather complex with three distinct regions. At the bottom very close to the surface, is a laminar sublayer which is stable due to the low velocity but has a high shear. At the edge near the invisicid external flow is turbulent decay layer or wake like region and joining these two in the middle is a highly turbulent mixing layer.

For the purposes of our momentum solution we use an approximate polynomial fit to this complex shape,

{}u over U_delta} = ( {y over delta} )^{ {H-1} over 2}
where the standard shape factor H= delta^* / theta is used. H is preset based on experimental data obtained from boundary layer tests with various pressure gradients.

subsonic5_html_m764e7d9d
For zero pressure gradient H = 9/7 so that a standard 1/7th power law profile is obtained. Using only small H variations gives reasonable estimate of profiles for accelerating, decelerating or near separation (H ≈ 2.6) turbulent boundary layers.

Treating H as a constant value and using the following experimental prediction for shear stress in a turbulent layer,

{tau_w over rho U_delta^2} =0.246 times 10^{-0.678H} R_{e theta}^{-0.268}
applying the integration of the polynomial approximations for a turbulent layer, leads to a governing momentum equation that is then in a similar format to the one found for laminar flow,

{partial theta ({rho over mu} theta U_delta)^{0.2} over partial x} = 0.0106 - {4 over U_delta}{partial U_delta over partial x} ({rho over mu} theta U_delta)^{0.2}
As in the case of the laminar layer, this turbulent layer equation can be solved by combining the partial derivatives and then integrating.

R_{e theta}^{6/5} = R_e . {0.0106 over (U_delta / V_infty)^3} int_{transition}^{x/c} ({U_delta over V_infty})^4 d {x over c}}
By adding the momentum loss of the turbulent layer to that already found for the initial laminar layer, a final value of momentum loss can be estimated for the flow over the surface of the aerofoil. This will include the sum of both upper and lower surface momentum losses. The total momentum loss over the surface of the section to the trailing edge will be called , θte at x/c = 1.

Turbulent boundary layer separation is predicted when H reaches or exceeds approximately 2.6. If separation hasn't occurred then nominally the H value at the trailing edge is 2.0, so that δ*=2θ at the trailing edge.

Drag Coefficient

If momentum loss over the whole section is known then drag of the section can be predicted.subsonic5_html_2dbd8e9cChange of flow momentum in-to/out-of the control volume is equal to the sum of the applied forces. The control volume is extended downstream where the pressure influence of the aerofoil is negligible so that outlet will have same stream static pressure as the inlet. Thus the only applied force component will be the support force to hold the section in place and thus will be equivalent to the drag (D) on the section.

The momentum balance will be,

- int_{inlet} rho u^2 dA + int_{outlet} rho u^2 dA = - D
also conservation of mass must be obeyed so,

int_{inlet} rho u. dA = int_{outlet} rho u. dA
Combining these two results gives and setting the inlet velocity to a constant stream velocity,

D = V_infty int_{inlet} rho V_infty .dA - int_{outlet} rho u^2 dA
D = V_infty int_{outlet} rho u .dA - int_{outlet} rho u^2 dA
assuming incompressible 2D flow so that dA = dy.1 (/per unit span),


D = rho V_infty^2 int_{outlet} ({u over V_infty} - ( {u over V_infty })^2 ) .dy
D = rho V_infty^2 int_{outlet} {u over V_infty} (1 - ( {u over V_infty }) ) .dy
non-dimensionalising to obtain drag coefficient gives,

C_D = { D over {1 over 2} rho V_infty^2 S} = 2 int_{outlet} {u over V_infty} (1 - ( {u over V_infty }) ) .d ({y over c} )
C_D =2 {{theta }_{wake} over c}
The section boundary layer prediction only gave a value of momentum loss up to the trailing edge. There is a small deficit between trailing edge and downstream wake due to the pressure imbalance at the trailing edge. The difference between trailing edge momentum loss and downstream wake momentum loss can be found by continuing the previous boundary momentum integration on downstream into the wake. The integral momentum equation is,

{partial theta over partial x}+(2+H){theta over U_delta}  {partial U_delta over partial x}- {tau_w over rho U_delta} = 0
in the wake τw=0 and H can be taken as approximately an average of the trailing edge value (H=2), and the far downstream value (H=1), the momentum equation then becomes,

{partial theta over partial x}+{ 7 over 2}{theta over U_delta}  {partial U_delta over partial x}} = 0
rearranging and integrating gives,

{ 1 over theta }  {partial theta over partial x} = -{7 over 2} { 1 over U_delta}  {partial U_delta over partial x }
int_{theta_te}^{theta_wake} {1 over theta } d theta  = -{7 over 2} int_{u_{te}}^{V_infty}{ 1 over U_delta}. d U_delta
ln (theta_{wake} )- ln (theta_{te})  = -{7 over 2} ( ln (V_infty)- ln (u_{te}))
theta_{wake} = theta_{te} ( { u_{te} over V_infty} )^{7/2}
The 7/2 power is dependent on our assumptions of the final shape of the turbulent boundary layer at the trailing edge and the reader may wish to investigate more detailed options to improve the accuracy of the final drag coefficient prediction. The final drag coefficient prediction for the aerofoil section will be

C_D  = 2 ( {theta_{te} over c} ) ( { u_{te} over V_infty} )^{7/2}