Gyredocs
Depth Practice
API

Gyre/Mathematics & Numerical Methods

Governing equations

In one minuteGyre solves the compressible Euler and Navier–Stokes equations in conservation form, non-dimensionalised by the free stream, and uses their entropy structure as the organising principle of the discretization.

Written at L2 Practice goes to L4 Theory 716 words Markdown

Gyre models a calorically perfect gas in d=2d = 2 or 33 space dimensions. The unknowns are the conserved variables

U=(ρρuρE),p=(γ−1)(ρE−12ρ ∣u∣2),T=pρR,\mathbf U = \begin{pmatrix} \rho \\ \rho \mathbf u \\ \rho E \end{pmatrix}, \qquad p = (\gamma - 1)\left(\rho E - \tfrac12 \rho\,|\mathbf u|^2\right), \qquad T = \frac{p}{\rho R},

and the governing equations are a system of conservation laws with a source s\mathbf s that is zero except in manufactured-solution studies:

∂tU+∇⋅F(U)−∇⋅Fv(U,∇W)=s.\partial_t \mathbf U + \nabla\cdot\mathbf F(\mathbf U) - \nabla\cdot\mathbf F_v(\mathbf U, \nabla \mathbf W) = \mathbf s .

For turbulent flow the same system is read as the Reynolds-averaged equations: the state gains one transported variable, ρν~\rho\tilde\nu, and the viscosity and conductivity gain an eddy contribution. That closure has its own page, Turbulence closure; everything below describes the laminar system it extends.

The inviscid flux in Cartesian direction bb and the viscous flux are

Fb=(ρubρu ub+p eb(ρE+p) ub),Fv,b=(0τ⋅bu⋅τ⋅b+κ ∂bT),\mathbf F_b = \begin{pmatrix} \rho u_b \\ \rho \mathbf u\, u_b + p\,\mathbf e_b \\ (\rho E + p)\, u_b \end{pmatrix}, \qquad \mathbf F_{v,b} = \begin{pmatrix} 0 \\ \boldsymbol\tau_{\cdot b} \\ \mathbf u\cdot\boldsymbol\tau_{\cdot b} + \kappa\, \partial_b T \end{pmatrix},

with the Newtonian stress under Stokes’ hypothesis and Fourier conduction,

τ=μ(∇u+∇uT−23(∇⋅u) I),κ=μ cpPr,cp=γRγ−1.\boldsymbol\tau = \mu\left(\nabla\mathbf u + \nabla\mathbf u^{\mathsf T} - \tfrac23 (\nabla\cdot\mathbf u)\,\mathbf I\right), \qquad \kappa = \frac{\mu\, c_p}{Pr}, \quad c_p = \frac{\gamma R}{\gamma - 1}.

The viscous flux depends on the gradient of only d+1d+1 quantities, the gradient variables W=(u,T)\mathbf W = (\mathbf u, T). The discretization lifts and differentiates exactly these.

Reference density, velocity and temperature are the free-stream values, and the reference length is one mesh unit (the chord of the built-in aerofoils, the diameter-independent unit of the cylinder mesh). Consequently

ρ∞=1,∣u∞∣=1,T∞=1,p∞=1γM∞2,R=1γM∞2,μ=1Re.\rho_\infty = 1,\quad |\mathbf u_\infty| = 1,\quad T_\infty = 1,\quad p_\infty = \frac{1}{\gamma M_\infty^2},\quad R = \frac{1}{\gamma M_\infty^2},\quad \mu = \frac{1}{Re}.

Force coefficients are therefore C=2F/SrefC = 2F / S_\text{ref}, and the pressure coefficient is cp=2 (p−p∞)c_p = 2\,(p - p_\infty). The Mach number, the angle of attack α\alpha (which sets the direction of u∞\mathbf u_\infty and of the lift/drag axes), the Reynolds number, γ\gamma and PrPr are entries of one flat parameter dictionary th; because it is a JAX pytree, each entry is a differentiable input without any further code.

The discretization is organised around the mathematical entropy

η(U)=−ρsγ−1,s=ln⁡p−γln⁡ρ,\eta(\mathbf U) = -\frac{\rho s}{\gamma - 1}, \qquad s = \ln p - \gamma \ln \rho ,

which is strictly convex for ρ,p>0\rho, p > 0. Its gradient defines the entropy variables

w=∂η∂U=(γ−sγ−1−ρ∣u∣22p,  ρup,  −ρp)T.\mathbf w = \frac{\partial \eta}{\partial \mathbf U} = \left( \frac{\gamma - s}{\gamma - 1} - \frac{\rho |\mathbf u|^2}{2p},\; \frac{\rho \mathbf u}{p},\; -\frac{\rho}{p} \right)^{\mathsf T}.

Contracting the Euler equations with w\mathbf w gives a scalar conservation law for smooth solutions, ∂tη+∇⋅(η u)=0\partial_t \eta + \nabla\cdot(\eta\,\mathbf u) = 0, because wT∂UFb=∂U(ηub)\mathbf w^{\mathsf T} \partial_{\mathbf U}\mathbf F_b = \partial_{\mathbf U}(\eta u_b). The entropy potential ψb=wTFb−ηub=ρub\psi_b = \mathbf w^{\mathsf T}\mathbf F_b - \eta u_b = \rho u_b is what the entropy-conservation condition for two-point fluxes is written in. With viscosity and heat conduction the entropy can only decay: ddt∫η≤0\tfrac{\mathrm d}{\mathrm dt}\int \eta \le 0 up to boundary terms. A scheme that reproduces this inequality discretely cannot blow up through the mechanism that destabilises under-resolved high-order methods (aliasing-driven entropy production of the wrong sign). Which parts of this Gyre reproduces discretely, and under which conditions, is stated with the interface fluxes: the inviscid volume terms always, the interface terms provably only with es_roe or llf, the viscous terms not at all. The unit test test_entropy_variables_are_gradient_of_entropy checks the formula above against jax.grad of η\eta.

key meaning default
gamma ratio of specific heats 1.4
mach free-stream Mach number 0.5
alpha angle of attack in radians (API accepts alpha_deg) 0
re Reynolds number per unit length; 103010^{30} when inviscid —
pr Prandtl number 0.72
p_back, pt_in, tt_in, t_wall boundary data relative to free-stream values case dependent
suth Sutherland constant S/T∞S/T_\infty; absent means constant viscosity absent

    ↑ ↓ moveEnter open[ ] change depthResults show section and depth level