Skip navigation links
(Wind-US Documentation Home Page) (Wind-US User's Guide) (GMAN User's Guide) (MADCAP User's Guide) (CFPOST User's Guide) (Wind-US Utilities) (Common File User's Guide) (Wind-US Installation Guide) (Wind-US Developer's Reference) (Guidelines Documents)

(General Overview) (General Purpose Commands) (File Specification Commands) (Domain Selection Commands) (Variable Selection and Unit Control Commands) (File and Report Generation Commands) (Plotting Commands) (GENPLOT File Format) (rake file Format) (Equations Used by CFPOST)

Equations Used by CFPOST

Variable Derivations

The following equations are used by CFPOST for deriving a requested variable that is not present in the solution file. If a variable cannot be derived by one of the following equations or by the simple multiplication or division by ρ, then the requested variable must exist within the solution file. Note that multiple derivations may be used to derive a variable. For example, when deriving the total temperature T0, the temperature T and Mach number M may themselves be derived.

The derivation of variables that include the gas constant R or the ratio of specific heats γ will take real gas effects into account if sufficient information is available. For frozen and finite-rate chemistry flows, species properties are read from a chemistry data (.chm) file specified with the chemistry command. Otherwise, real gas effects are modeled by multiplying R by the compressibility factor Z, and using the effective specific heat ratio β in place of γ.

Note that some variables require the specification of a surface, or normal direction into the flow field, and others require knowledge of "up" and "down" axes orientation. Variables denoted with (s) indicate that the normal direction must be specified properly using the subset or surface command. Variables denoted with (o) are affected by the axes orientation set by the orientation command.

Several variables may be computed in multiple ways, depending on the information available in the solution file. The equations are shown below in the order that CFPOST attempts to use them. If any of the needed quantities are missing from the solution file and cannot be derived, CFPOST will try using the subsequent equation, etc. Note though, that for frozen and finite-rate chemistry flows, if the necessary information isn't available CFPOST will abort with an error message, rather than use the perfect gas equations.

rho - Density

    ρ  =  p / (RTZ)
 =  p / (RT)

Z contains the local real gas effects, R/R.

rho0 - Stagnation density

For frozen and finite-rate chemistry,

    ρ0  =  p0 / (RT0Z)

For a perfect gas,

    ρ0  =  ρ [ 1 + (γ−1) M 2/2 ] 1 / (&gamma−1)

p - Pressure

    p  =  ρRTZ
 =  ρRT
 =  (β−1) (ρe0ρV 2/2− ρhf)
 =  (γ−1) (ρe0ρV 2/2− ρhf)

Z contains the local real gas effects, R/R.

p0 - Stagnation pressure

For frozen and finite-rate chemistry,

    p0  =  p exp{(1/Rmix) ∑nsi=1 Ci (Ru / MWi) [Sfuni,T0Sfuni,T]}

where Sfun is a function used to compute the entropy for a species as a function of temperature, from the curve fit data in the .chm file.

For a perfect gas,

    p0  =  p [ 1 + (γ−1) M 2/2 ] &gamma / (&gamma−1)

For a value in a rotating reference frame (i.e., p0r), the total temperature T0 and Mach number M in the above equations are replaced by the corresponding values in the rotating reference frame.

Cp - Pressure coefficient

    Cp  =  (pp) / q

where q = ρ(Ma)2/2

deltap - Delta pressure

    Δp  =  pp

pp - Pitot pressure

For M ≤ 1,

    pp  =  p [ 1 + (γ−1) M 2/2 ] &gamma / (&gamma−1)

For M > 1,

    pp  =  p [ (γ+1) M 2/2 ] &gamma / (&gamma−1) / [ 2γ M 2/(γ+1) − (γ−1)/(γ+1) ] 1 / (γ−1)

Cpt - Total pressure coefficient

    Cpt  =  (p0p) / q

where q = ρ(Ma)2/2

q - Dynamic pressure

    q  =  ρV 2/2

T - Temperature

    T  =  p / (ρRZ)
 =  p / (ρR)

Z contains the local real gas effects, R/R.

T0 - Stagnation temperature

For frozen and finite-rate chemistry, iterate on T0 until

    hmix,T0  =  h0,mix

    nsi=1 Ci (Ru / MWi) Hfuni,T0  =  nsi=1 Ci (Ru / MWi) Hfuni,T + V 2/2

where Hfun is a function used to compute the enthalpy for a species as a function of temperature, from the curve fit data in the .chm file.

For a perfect gas,

    T0  =  [ 1 + (γ−1) M 2/2 ]

For a value in a rotating reference frame (i.e., T0r), the values of T0, h0,mix, V 2, and M in the above equations are replaced by the corresponding values in the rotating reference frame.

h - Sensible enthalpy

    h  =  βei
 =  γei

h0 - Sensible stagnation enthalpy

    h0  =  βei + ek
 =  γei + ek

ei - Sensible internal energy per unit volume

    ei  =  p / [ρ (β−1)]
 =  p / [ρ (γ−1)]

ek - Kinetic energy per unit volume

    ek  =  V 2/2

e0 - Absolute stagnation energy per unit volume

    e0  =  (ρe0) / ρ
 =  [ ρ(e0)r + ρ(ω × r)2/2 − ρ(gr)/2 ] / ρ
 =  ei + ek + hf

Ma - Equivalent isentropic Mach number

    Ma  =  { [2/(β−1)] [ [1 + (β−1)M 2 / 2] [1 + βCpM 2 / 2](1−β)/β − 1 ] }1/2
 =  { [2/(γ−1)] [ [1 + (γ−1)M 2 / 2] [1 + γCpM 2 / 2](1−γ)/γ − 1 ] }1/2

where Cp is the pressure coefficient.

u - x velocity

In a non-rotating reference frame,

    u  =  (ρu)/ρ
 =  ur + (ω × r) ⋅ ex
 =  0

In a rotating reference frame,

    ur  =  (ρur)/ρ
 =  0

v - y velocity

In a non-rotating reference frame,

    v  =  (ρv)/ρ
 =  vr + (ω × r) ⋅ ey
 =  0

In a rotating reference frame,

    vr  =  (ρvr)/ρ
 =  0

w - z velocity

In a non-rotating reference frame,

    w  =  (ρw)/ρ
 =  wr + (ω × r) ⋅ ez
 =  0

In a rotating reference frame,

    wr  =  (ρwr)/ρ
 =  0

V - Velocity magnitude

    V  =  (V 2)1/2

For a value in a rotating reference frame (i.e., Vr), V 2 in the above equation is replaced by the corresponding value in the rotating reference frame.

The value of V 2 in a non-rotating reference frame is computed using one of the following equations. (As noted earlier, the equation used depends on the information available in the solution file. If any of the needed quantities in the first equation listed below are missing from the solution file and cannot be derived, CFPOST will try using the subsequent equation, etc.)

    V 2  =  |V2|
 =  u2 + v2 + w2
 =  (ρ |V|)2 / ρ2 = [(ρu)2 + (ρv)2 + (ρw)2] / ρ2
 =  (|Vr| + |ω × r|)2
 =  (ρ |Vr| + ρ |ω × r|)2 / ρ2

In a rotating reference frame, Vr2 is computed using one of the following equations.

    Vr 2  =  |Vr|2
 =  ur2 + vr2 + wr2
 =  (ρ |Vr|)2 / ρ2 = [(ρur)2 + (ρvr)2 + (ρwr)2] / ρ2

Vxy, Yyz, Vxz - Crossflow velocity magnitudes

    Vxy  =  (u2 + v2)1/2
Vyz  =  (v2 + w2)1/2
Vxz  =  (u2 + w2)1/2

M - Mach number

    M  =  V / a

For a value in a rotating reference frame (i.e., Mr), V in the above equation is replaced by the corresponding value in the rotating reference frame.

a - Speed of sound

    a  =  (γRTZ)1/2
 =  (γRT)1/2

s - Entropy

For frozen and finite-rate chemistry,

    ss  =  nsi=1 Ci (Ru / MWi) [Sfuni,TSfuni,T] − Rmix ln (p/p)

where Sfun is a function used to compute the entropy for a species as a function of temperature, from the curve fit data in the .chm file.

For a perfect gas,

    ss  =  (cv) ln [(p/p) / (ρ/ρ)γ]

where (cv) = R/(γ−1). Here s may be assumed to be zero.

S - Strain rate magnitude

    S  =  (2 Sij Sij)1/2

Sxx, Sxy, Sxz, Syx, Syy, Syz, Syz, Szx, Szy, Szz - Strain rate components

    Sij  =  (∂ui/∂xj + ∂uj/∂xi )/2

W - Rotation rate magnitude or vorticity magnitude

    Ω  =  (2 Wij Wij)1/2  =  (ωx2 + ωy2 + ωz2)1/2

The magnitude of the rotation rate tensor is equivalent to the vorticity magnitude.

Wxx, Wxy, Wxz, Wyx, Wyy, Wyz, Wyz, Wzx, Wzy, Wzz - Rotation rate components

    Wij  =  (∂ui/∂xj∂uj/∂xi )/2

wx, wy, wz - Vorticity components

    ωx  =  ∂w/∂y∂v/∂z
ωy  =  ∂u/∂z∂w/∂x
ωz  =  ∂v/∂x∂u/∂y

helicity - Helicity

    helicity  =  ΩV = ωxu + ωyv + ωzw

swirl - Swirl

    swirl  =  ΩV / (ρV 2) = (ωxu + ωyv + ωzw) / (ρV 2)

dila - Dilatation

    dila  =  ∇ ⋅ V = ∂u/∂x + ∂v/∂y + ∂w/∂z

localpha - Local α (o)

Let vds be the velocity in the "downstream" direction and vup be the velocity in the "upward" direction as defined by the aerodynamic axes set in the grid file or specified with the orientation command. Local α is then defined as:

    α  =  arctan (vup / vds)

locbeta - Local β (o)

Let vss be the velocity in the "sidestream" direction as defined by the aerodynamic axes set in the grid file or specified with the orientation command and let V be the magnitude of the velocity vector. Local β is then defined as:

    β  =  arcsin (vss / V)

mu - Total viscosity coefficient

    μ  =  μl + μt

mul - Laminar viscosity coefficient

For frozen and finite-rate chemistry,

    μl  =  nsi=1 [(Xi μi) / (∑nsj=1 Xj φij)]

where

    μi  =  μ0 (T/Tμ)3/2 [(Tμ + Sμ) / (T + Sμ)]
    φij  =  [ 1 + (μi / μj)1/2 (MWj / MWi)1/4 ]2 / ( 8 + 8MWi / MWj )1/2

The molecular weight of each species is obtained from the thermodynamic data in the chemistry (.chm) file, and the values of μ0, Tμ, and Sμ for each species are obtained from the transport data in the .chm file.

For a perfect gas,

    μi  =  μ0 (T/T0)3/2 [(T0 + S1) / (T + S1)]

where [Pratt and Whitney Aircraft Company "Aeronautical Vestpocket Handbook," 19th. ed., Part Number 79500]:

    μ0  =  0.37445×10−6 slug/ft-sec = 1.79287×10−5 kg/m-sec
    T0  =  518.67 °R = 288.15 K
    S0  =  216.0 °R = 120 K

kappa - Thermal conductivity coefficient

For frozen and finite-rate chemistry,

    κ  =  nsi=1 [(Xi κi) / (∑nsj=1 Xj φij)]

where

    κ  =  κ0 (T/Tκ)3/2 [(Tκ + Sκ) / (T + Sκ)]
    μi  =  μ0 (T/Tμ)3/2 [(Tμ + Sμ) / (T + Sμ)]
    φij  =  [ 1 + (μi / μj)1/2 (MWj / MWi)1/4 ]2 / ( 8 + 8MWi / MWj )1/2

The molecular weight of each species is obtained from the thermodynamic data in the chemistry (.chm) file, and the values of κ0, Tκ, Sκ, μ0, Tμ, and Sμ for each species are obtained from the transport data in the .chm file.

For a perfect gas,

    κ  =  μl cp / Pr

where cp = /(γ−1).

delta - Boundary layer thickness (s)

The edge of the boundary layer is determined according to the criteria established by the boundary layer edge command.

delta1 or delta* - Boundary layer displacement thickness (s)

    δ1  =  0δ [ 1 − (ρV) / (ρδVδ) ] ds = 0δ ds − 1 / (ρδVδ) 0δ ρV ds

where V is the total velocity magnitude.

Starting at the wall and advancing one point at a time along the grid line intersecting the wall, the integrals in the second equation are accumulated. At each grid point, the values of the subscripted variables are assumed to be the value at the current grid point. The value of the equation is then computed and compared to the value at the previous iteration. If the change meets the criteria specified by the boundary layer edge command, the iteration stops and the thickness will be set to the accumulated arc length along the grid line (i.e., the value of the first integral). (Algorithm compliments of R. H. Bush.)

delta2 or THETA - Boundary layer momentum thickness (s)

    δ2  =  0δ (ρV) / (ρδVδ) ( 1 − V / Vδ ) ds = 1 / (ρδVδ) 0δ ρV ds − 1 / (ρδVδ2) 0δ ρV 2 ds

The method described for the boundary layer displacement thickness delta1 is used to evaluate the integral.

delta3 - Boundary layer energy thickness (s)

    δ3  =  0δ (ρV) / (ρδVδ) ( 1 − V 2 / Vδ 2 ) ds = 1 / (ρδVδ) 0δ ρV ds − 1 / (ρδVδ3) 0δ ρV 3 ds

The method described for the boundary layer displacement thickness delta1 is used to evaluate the integral.

deltau - Boundary layer velocity thickness (s)

    δu  =  0δ [ 1 − V / Vδ ] ds = 0δ ds − 1 / Vδ 0δ V ds

The method described for the boundary layer displacement thickness delta1 is used to evaluate the integral. Note that the velocity thickness is the same as the incompressible displacement thickness.

THETAinc - Boundary layer incompressible momentum thickness (s)

    θinc  =  0δ V / Vδ [ 1 − V / Vδ ] ds = 0δ V / Vδ ds − 1 / Vδ2 0δ V2 ds

The method described for the boundary layer displacement thickness delta1 is used to evaluate the integral. The incompressible momentum thickness is often used to compute an incompressible shape factor for evaluating flow control devices.

Q - Heat transfer rate (s)

    Q  =  κ ∂T/∂n = κTn = κ (∂T/∂x ∂x/∂n + ∂T/∂y ∂y/∂n + ∂T/∂z ∂z/∂n)

The equation is evaluated using the same process described in the derivation of tau.

tau, taux, tauy, tauz - Total (laminar+turbulent) shear stress (s)

    τ  =  μ Vtan / ∂n

where μ is the total (laminar+turbulent) viscosity

    μ  =  μl + μt

and Vtan is the tangential velocity vector determined as follows:

    Vtan  =  V − (Vn) n

The magnitude of the shear stress is computed as

    τ  =  μ [(∂utan / ∂n)2 + (∂vtan / ∂n)2 + (∂wtan / ∂n)2]1/2

where the partial derivatives are evaluated using the chain rule:

    ∂f / ∂n  =  (∂f / ∂x) (∂x / ∂n) + (∂f / ∂y) (∂y / ∂n) + (∂f / ∂z) (∂z / ∂n)

The partial derivatives are computed using the equations presented in the section on "Derivatives of Variables With Respect to Cartesian Coordinates".

The shear stress acts in the direction of the velocity, so the components are:

    τx  =  τ utan / Vt
    τy  =  τ vtan / Vt
    τz  =  τ wtan / Vt

If wall function boundary conditions were used in the flow solver, CFPOST will also use wall functions when computing the shear stress.

Cf, Cfx, Cfy, Cfz - Skin friction coefficient (s)

The total skin friction coefficient is the shear stress magnitude normalized by the freestream dynamic pressure.

    Cf  =  τ / (ρV2/2)

Similarly, the Cartesian components of the skin friction coefficient are defined as

    (Cf )x  =  τx / (ρV2/2)
    (Cf )y  =  τy / (ρV2/2)
    (Cf )z  =  τz / (ρV2/2)

Note that the freestream values used are the reference values specified with the Wind-US FREESTREAM keyword and may differ substantially from those at the local boundary layer edge.

uu, uv, uw, vu, vv, vw, wu, wv, ww - Turbulent stress components (s)

    ui"uj"  =  {ρui"uj"} / ρ  =  −[2μt( SijSkkδij/3 ) − 2/3*ρkδij]

Turbulent stresses can be computed for any of the structured grid turbulence models. The formula above represents the Boussinesq approximation used by most models. For the Spalart-Allmaras turbulence model, the turbulent kinetic energy is estimated by means of a generalized form of Bradshaw's relation:

    ρk  =  μt (2SijSij)1/2 / 0.31
Algebraic Reynolds stress models use a more complex form for the stresses that involves non-linear expressions of the strain and rotation rate tensors.

y+ - Non-dimensional boundary layer wall coordinate (s)

    y+  =  y uτ / νwall

where

    uτ  =  (τwall / ρwall)1/2
    νwall  =  μwall / ρwall

and here y is the distance from the wall, measured along the grid line intersecting the wall.

u+ - Non-dimensional boundary layer velocity (s)

    u+  =  V / uτ

where uτ = (τwall / ρwall)1/2.

k+ - Non-dimensional turbulent kinetic energy (s)

    k+  =  k / uτ2

where uτ = (τwall / ρwall)1/2.

epsilon+ - Non-dimensional turbulent dissipation rate (s)

    ε+  =  ε μl / (ρ uτ4)

where uτ = (τwall / ρwall)1/2.

uu+, uv+, uw+, vu+, vv+, vw+, wu+, wv+, ww+ - Non-dimensional turbulent stress components (s)

    (ui"uj")+  =  {ρui"uj"} / (ρuτ2)

where uτ = (τwall / ρwall)1/2.

Mt - Turbulent Mach number

    Mt  =  (2 k / a2)1/2

Rt - Turbulent Reynolds number

[Ladd, J. A., and Kral, L. D. (1992) "Development and Application of a Zonal k-ε Turbulence Model for Complex 3-D Flowfields," AIAA/SAE/ASME/ASEE 28th. Joint Propulsion Conference and Exhibit, AIAA Paper 92-3176, p. 3]

    Rt  =  ρk2 / (μl ε)

Ry - Turbulent Reynolds number based on y (s)

[Ladd, J. A., and Kral, L. D. (1992) "Development and Application of a Zonal k-ε Turbulence Model for Complex 3-D Flowfields," AIAA/SAE/ASME/ASEE 28th. Joint Propulsion Conference and Exhibit, AIAA Paper 92-3176, p. 3]

    Ry  =  yρk1/2/μl

where here y is the distance from the wall, measured along the grid line intersecting the wall.

fmuJL, fmuCH, fmuSP, fmuLB - Turbulence model damping functions (s)

[Ladd, J. A., and Kral, L. D. (1992) "Development and Application of a Zonal k-ε Turbulence Model for Complex 3-D Flowfields," AIAA/SAE/ASME/ASEE 28th. Joint Propulsion Conference and Exhibit, AIAA Paper 92-3176, p. 3]

    (fμ)JL  =  e−2.5/(1 + Rt / 50)
    (fμ)CH  =  1 − e−0.0115y+
    (fμ)SP  =  (1 + 3.45 / Rt1/2) tanh (y+/70)
    (fμ)LB  =  (1 + 20.5 / Rt) (1 − e−0.0165Ry)2

C-species - Mass fraction of species

Using O2 as an example for species,

    C-O2  =  ρO2 / ρ

X-species - Mole fraction of species

Using O2 as an example for species,

    X-O2  =  (C-O2/MWO2) ∑nsi=1 MWi

Veff - Effective collision frequency

    Veff  =  C (p/T) Te1/2

where C is a function of the gas, p is the pressure in torr (i.e., millimeters of mercury), T is the temperature in K, and Te is the electron temperature. In CFPOST, C is assumed to be 7.3×109, the value for air, and Te is assumed equal to T.

MW - Molecular weight

For frozen and finite-rate chemistry,

    MW  =  1 / (∑nsi=1 Ci / MWi)

where Ci is the mass fraction of species i. The molecular weights for all the species may be listed using the show chemistry command.

For a perfect gas,

    MW  =  Ru / R

R - Gas constant

For frozen and finite-rate chemistry,

    R  =  Ru / MW

For a perfect gas,

    R  =  R

cp - Specific heat at constant pressure

For frozen and finite-rate chemistry,

    cp  =  nsi=1 Ci (Ru / MWi) Cpfuni,T

where Cpfun is a function used to compute the specific heat at constant pressure for a species as a function of temperature, from the curve fit data in the .chm file.

For a perfect gas,

    cp  =  γR / (γ − 1)

cv - Specific heat at constant volume

    cv  =  cpR

gamma - Ratio of specific heats

For frozen and finite-rate chemistry,

    γ  =  cp / cv

For a perfect gas,

    γ  =  γ

hf - Heat of formation

For frozen and finite-rate chemistry,

    hf  =  nsi=1 Ci [hf(298)]i

For a perfect gas,

    hf  =  0

Thermodynamic Functions

The functions Cpfun, Hfun, and Sfun appearing in some of the above equations for frozen and finite-rate chemistry flows are used to compute thermodynamic properties for a species as a function of temperature. Properties are computed using curve fit equations, with the coefficients supplied in a chemistry data (.chm) file. The .chm file to be used must be specified using the CFPOST chemistry command.

Three different (but related) formats are available for the type of curve fits being specified in the .chm file - SPARKCRV, WINDNASA, and NASA3287.

Cpfun - Molar specific heat at constant pressure

Function Cpfun computes the specific heat per mole at constant pressure for a species. For the SPARKCRV and WINDNASA formats, the curve fit formula is

    Cp / Ru  =  a1 + a2T + a3T2 + a4T3 + a5T4

For the NASA3287 format, the formula is

    Cp / Ru  =  a1T−2 + a2T−1 + a3 + a4T + a5T2 + a6T3 + a7T4 + a8T5

Note that the value returned by Cpfun is dimensionless.

Hfun - Molar absolute enthalpy

Function Hfun computes the enthalpy per mole for a species. For the SPARKCRV and WINDNASA formats, the curve fit formula is

    H / Ru  =  a1T + a2T2/2 + a3T3/3 + a4T4/4 + a5T5/5 + b1

For the NASA3287 format, the formula is

    H / Ru  =  a1T−1 + a2 ln T + a3T + a4T2/2 + a5T3/3 + a6T4/4 + a7T5/5 + a8T6/6 + b1

For the SPARKCRV format, the enthalpy reference state is at T = 0 K. I.e., the curve fit for H/Ru computes

    H / Ru  =  (Hf(0) + ΔH(0→T)) / Ru

where Hf(0) is the heat of formation at 0 K, and ΔH(0→T) is the change in enthalpy between 0 K and T.

However, for the WINDNASA and NASA3287 formats the enthalpy reference state is at T = 298.15 K. I.e., for these formats the curve fit computes

    H' / Ru  =  (Hf(298) + ΔH(298→T)) / Ru

where Hf(298) is the heat of formation at 298.15 K, and ΔH(298→T) is the change in enthalpy between 298.15 K and T.

Thus, for the WINDNASA and NASA3287 formats, function Hfun adjusts the value to shift the reference state to 0 K. I.e., the actual returned value is

    H / Ru  =  H' / Ru + H(shift) / Ru
     =  (Hf(298) + ΔH(298→T)) / Ru + (Hf(0)Hf(298) + ΔH(0→298)) / Ru
     =  (Hf(0) + ΔH(0→T)) / Ru

Note that the value returned by Hfun has units of K.

Sfun - Molar entropy

Function Sfun computes the entropy per mole for a species. For the SPARKCRV and WINDNASA formats, the curve fit formula is

    S / Ru  =  a1 ln T + a2T + a3T2/2 + a4T3/3 + a5T4/4 + b2

For the NASA3287 format, the formula is

    S / Ru  =  a1T−2/2 − a2T−1 + a3 ln T + a4T + a5T2/2 + a6T3/3 + a7T4/4 + a8T5/5 + b2

It should be noted that the entropy is defined by the equation

    dS  =  Cp dT/TRu dp/p

Thus

    S / Ru  =  (1/Ru) dS
     =  (Cp / Ru) dT/T dp/p

and Sfun computes the temperature portion of the above equation, ∫ (Cp / Ru) dT/T.

Note also that the value returned by Sfun is dimensionless.

Integration Formulas

The area to be integrated is considered to be a collection of individual polygons. A surface may be defined directly by a subset or surface command. If the dimensions of the surface is i points by j points, there will be (i−1) × (j−1) polygons that define the surface. The perimeter of the surface is described by the line segments that join adjacent grid points at the extreme bounds of the subset. A surface may also be defined as the intersection of a cutting plane defined by a cut command with a computational volume defined by a subset command. The polygons created by passing the cutting plane through the cells that make up the volume will be those processed by the integration.

For surfaces defined directly, each polygon will have exactly four verticies defined by the grid points that make up the polygon. The area will be determined by computing the magnitude of the cross products of the diagonals. For surfaces defined by a cutting plane, each polygon will have from three to six verticies defined by the intersection of the cutting plane with the cell edges. The area will be determined by decomposing the polygon into triangles and adding their individual areas. The unit normal vector for each polygon will be defined to be the unit vector in the direction of the vector sum of the normal vectors at the vertices.

For any cell, the value of a particular property for that cell (pressure, temperature, density, etc.) will be defined to be the average of the values at the verticies for that cell; i.e., for a cell with Nv vertices,

    fcell  =  (1/Nv) ∑ fi

where the summation is from i = 1 to Nv.

Area and Mass Average Integration

Averages over a group of Nc cells are computed as shown below. In these formulas the summations are from i = 1 to Nc, and the variables being summed are cell averages as defined above.

Area Average

    farea  =  (∑ fi Ai) / (∑ Ai)

Mass Average

    fmass  =  (∑ ρi fi Ai) / (∑ ρi Ai)

Mass Flux Average

    fmassflux  =  [ ∑ ρi fi (ViAi) ] / [ ∑ ρi (ViAi) ]

Arithmetic Mean

    f  =  (1/Nc) ∑ fi

    σ  =  { [ 1/(Nc−1) ] ∑ (fif )2 }1/2
     =  { [ 1/(Nc−1) ] [ (∑ fi2) − Ncf 2 ] }1/2

Flux Integration

Mass Flux

    m  =  ρ (Vn) A

Momentum Flux

    Fmomentum  =  ρV (Vn) A

Pressure Flux

    Fpressure  =  pA n

Total Thrust

    Ftotal  =  Fmomentum + Fpressure

Gross Thrust

    FG  =  (pp) A n

Force and Moment Integration (not yet done!)

    F  =  Fpressure + Fviscous

Derivatives of Variables With Respect to Cartesian Coordinates

Computing partial derivatives of variables with respect to Cartesian coordinates requires the calculation of the metric coefficients. Using subscripts to denote partial differentiation, the Jacobian is defined as:

    J = xξ (yηzζyζzη) + xη (yζzξyξzζ) + xζ (yξzηyηzξ)

The metric coefficients are then computed as follows:

    ξx = J −1 (yηzζyζzη)
    ξy = J −1 (xζzηxηzζ)
    ξz = J −1 (xηyζxζyη)
    ηx = J −1 (yζzξyξzζ)
    ηy = J −1 (xξzζxζzξ)
    ηz = J −1 (xζyξxξyζ)
    ζx = J −1 (yξzηyηzξ)
    ζy = J −1 (xηzξxξzη)
    ζz = J −1 (xξyηxηyξ)

Given the above, it is then possible to compute the partial derivative of a variable with respect to a Cartesian coordinate. Using the chain rule:

    fx = fξ ξx + fη ηx + fζ ζx

It is often necessary to know the vector normal to computational surface. First we define vectors tangent to a line in which only one of the computational coordinates vary:

    Vξ = xξ i + yξ j + zξ k
    Vη = xη i + yη j + zη k
    Vζ = xζ i + yζ j + zζ k

The vector normal to a constant computational surface is then the cross product of the two vectors whose computational coordinates vary on the surface:

    Nξ = Vη × Vζ
    Nη = Vζ × Vξ
    Nζ = Vξ × Vη

Given a normal vector, one can then simply form a unit normal vector and directional derivatives.

    n = N / |N| = (∂x/∂n) i + (∂y/∂n) j + (∂z/∂n) k

Derivatives of Variables With Respect to the Computational Coordinates

Derivatives of values with respect to computational coordinates may be either first or second order accurate depending on the need. Where normal vectors are not of a concern, second order central differences are used except where a boundary is involved and then first order forward or backward differences are used:

    ∂f/∂ξ ≈ δf = (fi+1fi−1) / 2     Central difference
    ∂f/∂ξ ≈ Δf = fi+1fi     Forward difference
    ∂f/∂ξ ≈ ∇f = fifi−1     Backward difference

For normal derivatives, first or second order may be selected by the first order or second order command or parameter on the integrate force command. Because of the internal structure of CFPOST, these derivatives are always forward derivatives:

    ∂f/∂ξfi+1fi     First order
    ∂f/∂ξ ≈ (−fi+2 + 4fi+1 − 3fi) / 4     Second order

Engine Analysis Calculations

Distortion - SAE Definition [Society of Automotive Engineers (1978) "Gas Turbine Engine Inlet Flow Distortion Guidelines," SAE-ARP-1420]

A low-pressure region is defined as a range of rake probes on a particular ring whose total pressures fall below the average total pressure for that ring.

Circumferential Intensity

    Intensity = (ΔPC / P)Ring i = (PAVRing iPAVLOWRing i) / PAVRing i

where

    PAVRing i = Ring-averaged total pressure
    PAVLOWRing i = Ring-averaged total pressure over all low-pressure regions

Circumferential Extent

    Extent = Cumulative angular extent of all low-pressure regions, in degrees

Multiple-per-revolution (MPR)

    MPRRing i = ∑Qk=1 [ (ΔPC / P)i k θik ] / max [ (ΔPC / P)i k θik ]

where

    Q = Number of low-pressure regions
    &thetaik = Extent of low-pressure region k on the i'th ring

Radial Intensity

    Intensity = (ΔPR / P)Ring i = (PFAVPAVRing i) / PFAV

where

    PFAV = Area-weighted face-average total pressure

Distortion - GE Definition

Circumferential Intensity

    Intensity = (PAVRing iPMINRing i) / PFAV

where

    PAVRing i = Area-averaged total pressure of ring i
    PMINRing i = Minimum total pressure of ring i
    PAV = Area-averaged total pressure over the complete face

Circumferential Extent

The circumferential extent is the same as the SAE definition.

Radial Intensity

The radial intensity is the same as the SAE definition.

Shape Factor

    Shape Factor = (GE ring circumferential intensity) / (SAE ring circumferential intensity)


Last updated 2 June 2016