SUBSONIC COMPRESSIBILITY CORRECTIONS

Potential flow solutions in two and three dimensions give accurate results for aerofoil and wing analysis provided the flow Mach number is less that 0.4. In this low subsonic region the flow is incompressible so that no density variations need to be considered in the governing equations. For higher Mach numbers the density does change, so an extra variable is introduced into the governing equations. The usual conservation of mass equation for two dimensional flow then becomes,

{partial rho u over partial x} + {partial rho v over partial y} = 0

With an extra unknown, an extra governing equation must be invoked to obtain a solution, so the conservation of momentum equation (Euler equ.) is added,

u {partial u over partial x} + v {partial u over partial y} =  { -1 over rho }   {partial P over partial x}
u {partial v over partial x} + v {partial v over partial y} =  { -1 over rho }   {partial P over partial y}
As momentum is a directional quantity, this produces two gradient equations. The flow is still assumed to be inviscid and these equations represent the change of flow momentum in a given direction due to the pressure forces in that direction.

This system of equations can be solved using complex numerical schemes (CFD) but as high speed aerofoil and wing flow is dominated by the stream direction with only small perturbations from this direction, simplifying assumptions can be made to convert the equations to a simple incompressible flow correction system.

The properties of air can be assumed to be a perfect, isentropic gas, P=rho R T (perfect gas law) and with a disturbance speed of a= sqrt{ gamma R T} (speed of sound in air). From the application of these isentropic flow conditions , P / rho^gamma = constant, the following gas dynamic relation (see section on Gasdynamics - Speed of Sound) can be applied,

{partial P over partial rho} = a^2
The terms of the continuity equation and the pressure gradient components of the Euler equation can be expanded as,

{partial P over partial x} = {partial Pover partial rho}.{partial rho over partial x}=a^2 {partial rho over partial x}
{partial P over partial y} = {partial Pover partial rho}.{partial rho over partial y}=a^2 {partial rho over partial y}
{partial rho u over partial x} = u {partial  rho over partial x} + rho {partial u over partial x}
{partial rho v over partial y} = v {partial  rho over partial y} + rho {partial v over partial y}
Substituting these terms into the continuity equation,and dividing by density, gives,

{partial u over partial x} - {u^2 over a^2}.{partial u over partial x}-{ u v over a^2}.{partial u over partial y}+{partial v over partial y}-{v^2 over a^2}{partial v over partial y}- {uv over a^2}{partial v over partial x} =0
as the flow is inviscid and hence irrotational,

{partial v over partial x}-{partial u over partial y} =0
hence continuity simplifies to,

{partial u over partial x}(1 - {u^2 over a^2} )-2{uv over a^2}{partial u over partial y}+{partial v over partial y} (1 -{v^2 over a^2} )  =0
Assuming that the stream velocity,V_infty, dominates the flow, then this allows the definition of perturbation velocities, u' and v' such that

u = V_infty + u^' v= v^'
where the prime values represent only small directional disturbances to the stream velocity. Substituting this definition into the continuity equation gives,

{partial u^' over partial x}(1 - {(V_infty+u^')^2 over a^2} )-2{(V_infty+u^')v^' over a^2}{partial u^' over partial y}+{partial v^' over partial y} (1 -{v^{'2} over a^2} )  =0.
If only terms of primary magnitude are kept,

{partial u^' over partial x}(1 - {V_infty^2 over a^2} ) + { partial v^' over partial y} =0
and if the analysis is extended to three dimensions,

{partial u^' over partial x} (1 - {M_infty^2  ) + { partial v^' over partial y} + {partial w^' over partial z}=0.
This is similar in form to the original governing continuity equation for incompressible flow but with an additional multiplier on the x-direction component. Its solution can be obtained in a similar fashion to the incompressible flow solution procedure.

A perturbation velocity potential phi^'can be defined such that,

{ partial phi^' over partial x} = u^' { partial phi^' over partial y} = v^' { partial phi^' over partial z} = w^'
and on substitution the governing equation becomes,

(1-M_infty^2) {partial^2 phi^' over partial x^2} + {partial^2 phi^' overpartial y^2}+{ partial^2 phi^' over partial z^2} =0.
Without the constant multiplier for the x-direction gradient, this equation would be exactly the same as the potential flow governing equation. A simple solution technique is to apply a geometric scaling transform to the x-direction so that the multiplier is canceled in the second flow field and the solution for the new geometry will be carried out using just potential flow.

The transform scaling multiplier will be beta = (1-M_infty^2) . This will be applied to x-dirn only. The transformed geometry will be x^' = x/beta, y^' = y and z^' =z. This will lead to a new geometry with the same thickness and span but a longer chord as the scaling factor is always less than 1.

Substituting into the governing compressible flow equation gives,

beta^2 {partial^2 phi^' over partial (beta x^')^2} + {partial^2 phi^' over partial y^{'2}} + {partial^2 phi^' over partial z^{'2}} = 0
 {partial^2 phi^' over partial x^{'2}} + {partial^2 phi^' over partial y^{'2}} + {partial^2 phi^' over partial z^{'2}} = 0.
This standard Laplace equation in perturbation potential for the transformed flow field can be solved using the usual potential flow techniques for aerofoil sections and wings as shown previously in this chapter.

Note that for two dimensional, section solutions, the transformed geometry is just a smaller thickness to chord ratio version of the original geometry. In most cases this can be assumed to have the same incompressible “thin aerofoil” solutions as the original geometry. For three dimensional wings, the new geometry has a reduced aspect ratio so these solutions will need to be predicted separately for this geometry.

In all cases, inviscid, incompressible flow solutions for velocities u^' , v^' , w^' and then pressures and forces C_p^' , C_L^'...etc, can be easily found. These then need to be transformed back (inverse transform) to find their equivalent values for the original compressible flow field. For velocity,

u^'(x,y,z) = {partial phi^' over partial x}
u^'(x^',y^',z^') = {partial phi^' over partial x^'}={partial phi^' over partial beta x}
thus the horizontal perturbation velocity in the original flow field is scaled by the factor beta compared to the component in the incompressible flow field. The other component velocities are unchanged.

u^'(x,y,z) = { 1 over beta} u^'(x^',y^',z^')
subsonic8_html_4436d701
To determine the pressure scaling, a compressible version of the Bernoulli equation (based on conservation of momentum) can be applied along a stream tube flowing over the aerofoil.

{ V_infty^2 over 2 } - { V^2 over 2} = int_{P infty}^P   { 1 over rho} . dP
Assuming small perturbation properties,

V^2 = (V_infty + u^')^2 + v^{'2} + w^{'2}   approx  V_infty^2 +2 V_infty u^'
and rho = rho_infty + rho^'   approx   rho_infty

Substituting and expanding the integration gives,

{V_infty^2 over 2}- ( { V_infty^2 over 2} - {2 V_infty u^' over 2} ) = {P - Pinfty over  rho_infty}
Hence the pressure coefficient for the compressible flow will be,

C_P = {P - P_infty over {1 over 2} rho_infty V_infty^2 }  =  - {2 u^' over V_infty }
Comparing this between the two parts of the transformed flow gives the mapping of pressure coefficient between the two.

C_P =  - {2 u^'(x,y,z) over V_infty } = -{1 over beta} {2 u^'(x^',y^',z^') over V_infty}  = { C_P^' over beta}
The compressible flow field coefficients can be found by scaling the incompressible flow field coefficients by the factor, { 1 over beta} = {1 over sqrt{ 1- M_infty^2}}.

The force coefficients will scale in the same manner as the pressure coefficients. Although the lengths are different between the two flow fields, the non-dimensionalising of the length scales will correct for this. In fact, this potential flow result could be generalised to all force coefficients for a subsonic compressible flow, if the incompressible results are known.

C_L = {C_L(incompressible) over sqrt{1-M_infty^2 }}
C_D = {C_D(incompressible) over sqrt{1-M_infty^2 }}
C_M = {C_M(incompressible) over sqrt{1-M_infty^2 }}

Limiting Case, Critical Mach Number.

This simple compressible flow correction theory works reasonably well up to the point where supersonic flow starts to appear near the surface of the aerofoil section. This will happen before the free-stream becomes supersonic due to the acceleration of the air in the vicinity of the aerofoil. Supersonic flow obeys different physical rules as the flow is moving faster than disturbances can propagate within the gas, (speed of sound) and the assumptions of a simple perturbation theory are no longer are valid. The free-stream Mach number for which supersonic flow first occurs on the wing or section is called the critical Mach number. Not only is this an important limit for theory but also marks the start of transonic flow and the likely-hood of a significant drag rise for the section.

In order to predict the critical Mach number, the above compressibility correction factor can be combined with standard gas dynamic equations. This is not a linear first order problem so analysis may be a bit difficult.

Sample pressure fields for NACA 0012 aerofoil section at 0 deg angle of incidence and a range of subsonic and transonic Mach numbers are shown in the following figure. High pressure is shown as red, low pressure as blue.

Critical Mach number is approx 0.7 for this case. A significant difference in flow structure can be observed between high subsonic flow (0-->0.65) and transonic flow (0.75-->1.5).

subsonic8_transonic
Images produced by DSMC 2-D numerical experiment © 2008 Auld.
Critical Mach number will vary for aerofoils and wings depending on their angle and camber as this determines the minimum pressure coefficient point, the point at which the local flow is fastest. This will be the point at which the flow first becomes supersonic.

Minimum pressure coefficient for incompressible flow, C_P^' can be found using potential flow techniques. By applying the subsonic correction factor, minimum pressure coefficient for compressible flow can be found,

C_P(min) = { C_P^'(min) over sqrt{1 -M_infty^2}}
However as speed where M_infty= M_{crit} is still unknown, this leads to a complex second order problem.

The variation of pressure along the stream tube can be used to provide another estimate of CP. The following is the full gas dynamics equation (see Compressible Flow Section – Isentropic Relations.)

subsonic8_html_10e0ef7d
{P over P_infty} =( {1+ {gamma-1 over 2}M_infty^2 over 1+{gamma-1 over 2}M^2} )^{gamma over gamma-1}

.For the point where P=P(min) at M=1, then this corresponds to the condition when free-stream Mach number is Mcrit ,

{P_{min} over P_infty} =( {1+ {gamma-1 over 2}M_infty^2 over 1+{gamma-1 over 2}} )^{gamma over gamma-1}Thus a second expression for CP(min) is,

C_P (min) = { P_{min} - P_infty over {1 over 2} rho_infty V_infty^2} = { P_{min} - P_infty over {1 over 2} gamma P_infty M_infty^2}
C_P (min) = {2 over gamma M_{crit}^2} ({ P_{min} over P_infty} -1 )
Substituting the compressibility correction factor and the above expression for pressure ratio, leads to the following equation.

{C_P^' (min) over sqrt{1 - M_{crit}^2}}= {2 over gamma M_{crit}^2} ( ({2 + (gamma-1) M_{crit}^2 over gamma+1} )^{gamma over gamma-1} -1 ).
Given a value of the incompressible flow minimum pressure coefficient for a wing or section at set free-stream flow angle, it should theoretically be possible to solve for the critical Mach number of the section at that condition. However as the equation is highly non-linear and difficult to invert, the solution is obtained either by an iterative process or by plotting minimum pressure coefficient versus critical Mach number and using empirical data lookup or curve fitting.

The following is a plot of minimum incompressible pressure coefficient versus critical Mach number.

subsonic8_html_32894a25