Background
This page summarizes the formulation and the numerical methods of the default, steady Euler model. The RANS model, time integration and the characteristic far fields build on it and are described in RANS, time integration and far fields. For how to use them, see the user guide.
Governing equations
The flow is governed by the steady compressible Euler equations for the conservative state \(U = (\rho, \rho\mathbf{u}, \rho E)\):
The freestream flows in the \(+x\) direction at angle of attack \(\alpha\): its direction is \((\cos\alpha, \sin\alpha)\) in 2D (x chordwise, y vertical) and \((\cos\alpha, 0, \sin\alpha)\) in 3D (x chordwise, y spanwise, z vertical). The freestream state is set by the user’s \(\rho_\infty\), \(p_\infty\) and Mach number \(M_\infty\).
Discretization
The state is approximated by piecewise-constant (p=0) DG functions. For every test function \(v\) the residual is
At p=0 the volume term vanishes and the scheme is a cell-centered finite-volume method.
The interior numerical flux \(\hat{H}\) is HLL [HLvL83] by default; HLLC
[TSS94] and local Lax-Friedrichs are available (flux_function of
CompressibleEulerModel). The integrals are assembled with FEniCSx, and the Jacobian
\(\partial R/\partial U\) is obtained by automatic differentiation of the weak form.
Boundary conditions
Boundary conditions are imposed weakly. The boundary flux is the normal flux of the interior state plus a local Lax-Friedrichs penalty toward a boundary state \(U_b\): \(\hat{H}_b = \mathbf{F}(U)\cdot\mathbf{n} + \lambda\,(U - U_b)\).
Boundary |
\(U_b\) |
|---|---|
Far field |
characteristic state from the Riemann invariants of the interior and the freestream (see far fields) |
Slip wall |
interior state with the normal velocity reflected, \(\mathbf{u}_b = \mathbf{u} - 2(\mathbf{u}\cdot\mathbf{n})\mathbf{n}\) |
Symmetry plane (3D, \(y=0\)) |
as the slip wall |
The whole outer boundary is one far field; no inflow/outflow division is made, so the angle of attack can change freely.
Viscous (RANS) models, pseudo-transient continuation and unsteady time stepping, and the characteristic far field are described in RANS, time integration and far fields.
Nonlinear and linear solvers
\(R(U) = 0\) is solved with Newton’s method (PETSc SNES, full steps). Each update \(\delta U\) is scaled by a
step length \(\theta \le 1\) that keeps density and pressure above a floor (\(10^{-4}\)) in every cell: \(\theta\)
starts at \(\min\left(1, (\|\delta U\| / \max(\|U\|, 1))^{-1/2}\right)\) and is halved until the updated state is
admissible on all MPI ranks. The solve stops when the residual norm drops below euler_residual_conv_limit.
Each Newton system is solved inexactly with GMRES (relative tolerance \(10^{-3}\)), preconditioned with additive
Schwarz (overlap asm_overlap) and ILU(ilu_levels) on the subdomains. The adjoint and tangent systems use the
same preconditioner with tight tolerances.
Optionally, a reduced-order model replaces the full-order solve when it is accurate enough: a least-squares Petrov-Galerkin solve [CBMF11] in a POD basis of earlier solutions, see Reduced-order modeling.
Shape parameterization
The wall is embedded in a B-spline free-form deformation (FFD) block [SP86], built with lsdo_geo. Wall node positions are a linear B-spline map of the FFD control points, so design variables that move control points, either individually or through sectional operations, deform the wall smoothly. See Shape parameterization.
Mesh deformation
The volume mesh follows the wall with IDWarp [LCB12]: each volume node moves by an
inverse-distance-weighted average of the surface node displacements and rotations. The far-field
boundary is included as a fixed surface, so the far field stays in place. The weighted sums are evaluated with
a kd-tree (the idwarp_jax package) and pruned with a relative tolerance err_tol (default \(10^{-6}\)).
Outputs
The pressure force on the wall is \(\mathbf{F}_p = \int_{wall} p\,\mathbf{n}\,ds\). Drag \(D\) and lift \(L\) are its
components along and normal to the freestream, and each moment \(M = \mathbf{d}\cdot(\mathbf{r}\times\mathbf{F}_p)\) is taken about a user-given point along a
user-given direction moments=[(name, point, direction), ...] (the pitching moment about a point, positive
nose-up, is direction "pitch"; the root bending moment of a half-wing is direction \((1, 0, 0)\) about a point
on the root). The coefficients divide these by \(\tfrac{1}{2}\rho_\infty|\mathbf{u}_\infty|^2 S_{ref}\) (and a
reference length for the moments, the reference chord by default), with \(S_{ref} = 1\) in 2D and the planform area of the half-wing in 3D, unless set by
the user.
Derivatives
CSDL composes the model as
and computes total derivatives in reverse mode. The flow solve is an implicit operation: for an output \(f\) it solves the adjoint system \((\partial R/\partial U)^T \psi = (\partial f/\partial U)^T\) and propagates \(-\psi^T \partial R/\partial \mathbf{x}_{mesh}\) and \(-\psi^T \partial R/\partial \alpha\) upstream through the mesh deformation and the FFD. Derivatives are provided for \(L\), \(D\) and \(M\); the coefficients are monitoring outputs only.
Bibliography
Kevin Carlberg, Charbel Bou-Mosleh, and Charbel Farhat. Efficient non-linear model reduction via a least-squares Petrov–Galerkin projection and compressive tensor approximations. International Journal for Numerical Methods in Engineering, 86(2):155–181, 2011. doi:10.1002/nme.3050.
Amiram Harten, Peter D. Lax, and Bram van Leer. On upstream differencing and Godunov-type schemes for hyperbolic conservation laws. SIAM Review, 25(1):35–61, 1983.
Edward Luke, Eric Collins, and Eric Blades. A fast mesh deformation method using explicit interpolation. Journal of Computational Physics, 231(2):586–601, 2012.
Thomas W. Sederberg and Scott R. Parry. Free-form deformation of solid geometric models. ACM SIGGRAPH Computer Graphics, 20(4):151–160, 1986.
Eleuterio F. Toro, Michael Spruce, and William Speares. Restoration of the contact surface in the HLL-Riemann solver. Shock Waves, 4(1):25–34, 1994.