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 / (R∞TZ) | |
= | p / (R∞T) |
Z contains the local real gas effects, R/R∞.
rho0 - Stagnation density
For frozen and finite-rate chemistry,
ρ0 | = | p0 / (R∞T0Z) |
For a perfect gas,
ρ0 | = | ρ [ 1 + (γ∞−1) M 2/2 ] 1 / (&gamma∞−1) |
p - Pressure
p | = | ρR∞TZ | |
= | ρR∞T | ||
= | (β−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,T0 − Sfuni,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 | = | (p−p∞) / q∞ |
deltap - Delta pressure
Δp | = | p − p∞ |
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 | = | (p0 − p∞) / q∞ |
q - Dynamic pressure
q | = | ρV 2/2 |
T - Temperature
T | = | p / (ρR∞Z) | |
= | 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 − ρ(g ⋅ r)/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 | = | (γ∞R∞TZ)1/2 | |
= | (γ∞R∞T)1/2 |
s - Entropy
For frozen and finite-rate chemistry,
s − s∞ | = | ∑nsi=1 Ci (Ru / MWi) [Sfuni,T − Sfuni,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,
s − s∞ | = | (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 |
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 = κ ∇T ⋅ n = κ (∂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 − (V ⋅ n) 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 | = | τ / (ρ∞V∞2/2) |
Similarly, the Cartesian components of the skin friction coefficient are defined as
(Cf )x | = | τx / (ρ∞V∞2/2) | |
(Cf )y | = | τy / (ρ∞V∞2/2) | |
(Cf )z | = | τz / (ρ∞V∞2/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( Sij − Skkδ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 |
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 | = | cp − R |
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 |
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/T − Ru 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.
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.
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 (Vi ⋅ Ai) ] / [ ∑ ρi (Vi ⋅ Ai) ] |
Arithmetic Mean
f | = | (1/Nc) ∑ fi |
σ | = | { [ 1/(Nc−1) ] ∑ (fi − f )2 }1/2 | |
= | { [ 1/(Nc−1) ] [ (∑ fi2) − Ncf 2 ] }1/2 |
Mass Flux
m | = | ρ (V⋅n) A |
Momentum Flux
Fmomentum | = | ρV (V⋅n) A |
Pressure Flux
Fpressure | = | pA n |
Total Thrust
Ftotal | = | Fmomentum + Fpressure |
Gross Thrust
FG | = | (p − p∞) A n |
F | = | Fpressure + Fviscous |
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 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+1 −
fi−1) / 2
Central difference
∂f/∂ξ ≈ Δf =
fi+1 − fi
Forward difference
∂f/∂ξ ≈ ∇f =
fi − fi−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+1 − fi
First order
∂f/∂ξ ≈
(−fi+2 +
4fi+1 −
3fi) / 4
Second order
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 i − PAVLOWRing 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 = (PFAV − PAVRing i) / PFAV
where
PFAV = Area-weighted face-average total pressure
Circumferential Intensity
Intensity = (PAVRing i − PMINRing 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