21. Geometric Conservation Law for Moving Meshes#

21.1. Introduction#

For cases where a moving mesh is employed, it is essential to ensure that the numerical scheme adheres to the Geometric Conservation Law (GCL). The GCL guarantees that uniform flow conditions are preserved when the mesh moves, preventing the introduction of spurious numerical artifacts.

21.2. Overview#

This appendix documents the GCL enforcement used by Stream’s moving-mesh implementation.

At a high level, the approach is:

  • Use an ALE (Arbitrary Lagrangian–Eulerian) finite-volume form, where the convective fluxes use the relative (fluid minus grid) velocity.

  • Compute a GCL-consistent grid volumetric flux (\(\phi_g\)) on every face from the old and new mesh geometry (i.e., from the swept volume).

  • Use the same time-averaged face area vectors in the pressure–velocity coupling (momentum interpolation) so that the moving-mesh correction is consistent with the rest of the discretization.

21.3. ALE Form and the Geometric Conservation Law#

Consider a conservation law for state (\(\vec{U}\)) on a moving control volume V(t) with boundary \(\partial V(t)\):

(21.3.1)#\[\frac{d}{dt}\int_{V(t)} \vec{U}\,dV + \oint_{\partial V(t)} \left(\vec{u} \vec{U} - \vec{w} \vec{U} + \vec{G}\right)\cdot d\vec{S} = \int_{V(t)} \vec{Q}\,dV\]

where (\(\vec{w}\)) is the grid velocity, \(\vec{u}\vec{U} + \vec{G}\), is the physical flux, and \(\vec{S}\) is the outward area vector. This can be written more compactly as

(21.3.2)#\[\vec{F}(\vec{U}) = \vec{u}\vec{U} + \vec{G}\]

And the area vector differential as \(d\vec{S} = \vec{n} dA\), where \(\vec{n}\) is the outward unit normal.

Then Eq. (21.3.1) becomes

(21.3.3)#\[\frac{d}{dt}\int_{V(t)} \vec{U}\,dV + \oint_{\partial V(t)} (\vec{F}(\vec{U}) - \vec{w}\vec{U}) \cdot \vec{n}\,dA = \int_{V(t)} \vec{Q}\,dV .\]

For a uniform (constant) solution (\(\vec{U}=\vec{U}_0\)) with no sources, the physical flux is constant and the scheme should not generate spurious changes purely due to mesh motion. Enforcing this leads to the geometric conservation law:

(21.3.4)#\[\frac{dV}{dt} = \oint_{\partial V(t)} \vec{w}\cdot \vec{n} dA\]

In a face-based finite-volume discretization, Eq. (21.3.4) becomes a requirement that the discrete volume change is exactly matched by the sum of discrete grid volumetric fluxes on the cell faces.

21.4. Discrete GCL Statement#

21.4.1. Continuous-to-discrete form#

If we consider a polyhedral control volume with volume \(V(t)\) and discrete boundary faces, the geometric conservation law in Eq. (21.3.4) can be written as:

(21.4.1.1)#\[\frac{dV}{dt} = \sum_{faces} \int_{face} \vec{w} \cdot \, d\vec{S} = \sum_{faces} G_{f}\]

Where \(V\) is the volume of the control volume, \( \vec{w} \) is the velocity of the control volume surface, and \(d\vec{S}\) is the outward normal area vector on the surface. The area vector is defined such that \( d\vec{S} = \hat{n} \, dA \) where \( \hat{n} \) is the outward normal unit vector and \( dA \) is the differential area. This is a vector that points normal to the surface with magnitude equal to the area of the surface element (\(\vec{S} = \vec{n} \cdot A\)). The \(G_{f}\) is the volumetric flow rate (sometimes called “flux” even though it isn’t a per-unit-area quantity) through face \(f\) due to the motion of the control volume surface (units are volume per time i.e. \(m^3/s\)). For clarity, the definition of \(G_{f}\) is:

(21.4.1.2)#\[G_{f} = \int_{face} \vec{w} \cdot \, d\vec{S}\]

If we discretize the time derivative of the control volume using a BDF (backward Differentiation Formula) scheme, we have:

(21.4.1.3)#\[\frac{1}{\Delta t} \sum_{j=0}^{N} \alpha_j V^{n+1-j} = \sum_{faces} \Delta V_{face}\]

where \( \Delta t \) is the timestep size, \( N \) is the order of the BDF scheme, and \( \alpha_j \) are the BDF coefficients. This is the generic BDF discretization formula.

To connect the RHS to a “midpoint” flux, integrate the continuous GCL over the step:

(21.4.1.4)#\[V^{n+1}-V^{n} \;=\; \int_{t^n}^{t^{n+1}} \sum_{f} G_f(t)\,dt \;=\; \sum_{f} \Delta V_{f},\]

where \(\Delta V_{f} \equiv \int_{t^n}^{t^{n+1}} G_f(t)\,dt\) is the swept volume through face \(f\) in one step. Dividing by \(\Delta t\) defines the time-averaged face grid flux

(21.4.1.5)#\[G_f^{\,n+1/2} \;\equiv\; \frac{1}{\Delta t}\int_{t^n}^{t^{n+1}} G_f(t)\,dt .\]

If one applies the midpoint rule (or if \(G_f(t)\) is approximately linear in time), then \(G_f^{\,n+1/2} \approx G_f(t^{n+1/2})\). In Stream, \(G_f^{\,n+1/2}\) is computed from the swept geometry, so the notation denotes the time average rather than a literal point evaluation.

The discrete GCL condition for a first-order (BDF1) time step (\(N=1\) and \(\alpha_0 = 1\), \(\alpha_1 = -1\)) is:

(21.4.1.6)#\[\frac{V^{n+1}-V^{n}}{\Delta t} = \sum_{f} G_{f}^{\,n+1/2} ,\]

where (\(G_{f}^{\,n+1/2}\)) is a time-step averaged grid flux over \([t^n,t^{n+1}]\).

For BDF2 with constant \(\Delta t\), the time derivative of volume is

(21.4.1.7)#\[\left.\frac{dV}{dt}\right|^{n+1} \approx \frac{\tfrac{3}{2}V^{n+1}-2V^n+\tfrac{1}{2}V^{n-1}}{\Delta t} \;=\; \frac{3}{2}\frac{V^{n+1}-V^{n}}{\Delta t} - \frac{1}{2}\frac{V^{n}-V^{n-1}}{\Delta t}.\]

In Stream, the per-step grid flux computed from geometry is stored at each timestep. Therefore, for BDF2 the face grid flux used in the solver is formed as:

(21.4.1.8)#\[G_{f}^{\,\text{BDF2}} =\ \frac{3}{2}\,G_{f}^{\,n+1/2} - \frac{1}{2}\,G_{f}^{\,n-1/2} \qquad(\text{constant }\Delta t),\]

21.5. Face-Based GCL: Facet Decomposition and Time-Averaged Area Vectors#

21.5.1. Why the face is decomposed into facets#

On an unstructured mesh, faces are generally polygons. A common pitfall is to compute (\(\phi_{g,f}\)) using a single face-normal/area and a single representative face velocity, which can violate the discrete GCL when the face deforms.

Each polygonal face is decomposed into triangular facets formed by the face center and adjacent node pairs on a polyhedral face. For a face with nodes \(\{\vec{P}_i\}_{i=0}^{N-1}\) and center \(\vec{C}\), facet \(i\) uses the triangle \((\vec{C},\vec{P}_i,\vec{P}_{i+1})\) (cyclic indexing).

21.5.2. Time-averaged facet area vector (matches facetNormal)#

Assume that, over a time step, each vertex position varies linearly in time (equivalently: vertex velocity is constant over the step).

For one triangular facet, define the edge vectors (following the ordering used in src/movingMesh.loci):

(21.5.2.1)#\[\vec{u}(t) = \vec{P}_1(t)-\vec{C}(t),\qquad \vec{v}(t) = \vec{P}_0(t)-\vec{C}(t).\]

The instantaneous facet area vector is

(21.5.2.2)#\[\vec{S}(t) = \frac{1}{2}\,\vec{u}(t)\times\vec{v}(t).\]

With linear-in-time motion, \(\vec{u}(t)\) and \(\vec{v}(t)\) are linear in time, so \(\vec{S}(t)\) is quadratic. The time-averaged facet area vector over the step is

(21.5.2.3)#\[\bar{\vec{S}} \;=\; \frac{1}{\Delta t}\int_{t^n}^{t^{n+1}} \vec{S}(t)\,dt.\]

Evaluating Eq. (21.5.2.3) in closed form (using only the old and new vertex coordinates) yields the expression implemented as facetNormal:

(21.5.2.4)#\[\bar{\vec{S}} \;=\; \frac{1}{6}\Big( \vec{u}^n\times\vec{v}^n \;+\; \vec{u}^{n+1}\times\vec{v}^{n+1} \;+\; \frac{1}{2}\left( \vec{u}^n\times\vec{v}^{n+1} \;+\; \vec{u}^{n+1}\times\vec{v}^{n} \right) \Big).\]

When the mesh is stationary \((\vec{u}^{n+1}=\vec{u}^n, \vec{v}^{n+1}=\vec{v}^n)\), Eq. (21.5.2.4) reduces to the standard triangle area vector \(\bar{\vec{S}}=\vec{S}\).

Summing \(\bar{\vec{S}}\) over all facets on a face gives the GCL-consistent time-averaged face area vector, stored in facetNormalSum{...}.

21.5.3. Grid volumetric flux on a face (matches facetVDotNormalSum)#

We compute the vertex velocities for the step as finite differences:

(21.5.3.1)#\[\vec{w}_{\vec{P}} = \frac{\vec{P}^{n+1}-\vec{P}^{n}}{\Delta t},\qquad \vec{w}_{\vec{C}} = \frac{\vec{C}^{n+1}-\vec{C}^{n}}{\Delta t}.\]

In the code, \(\vec{w}_{\vec{C}}\) is computed from the actual stored face-center coordinates (facecenter{n} and facecenter{n,it}), rather than by averaging nodal velocities, to remain consistent with how \(\vec{C}\) is defined.

For each facet, a centroid velocity is formed as

(21.5.3.2)#\[\vec{w}_\text{facet} = \frac{\vec{w}_{\vec{P}_0}+\vec{w}_{\vec{P}_1}+\vec{w}_{\vec{C}}}{3}.\]

The face grid volumetric flux is then computed as the facet sum

(21.5.3.3)#\[\phi_{g,f}^{\,n+1/2} \;=\; \sum_{\text{facets }k\in f} \vec{w}_{\text{facet},k}\cdot \bar{\vec{S}}_k .\]

This value is stored as facetVDotNormalSum{...} and exported as movingMesh::gridVolFlux{...}.

21.6. How This Appears in src/movingMesh.loci#

Key quantities (naming as in the code):

  • Mesh positions: pos{n} (start of time step) and pos{n,it} (current mesh at nonlinear iteration).

  • Nodal mesh velocity: node_v{n,it} = (pos{n,it}-pos{n})/dt{n}.

  • Time-averaged face area vector: facetNormalSum{n,it} (sum of facet (\bar{\vec{S}})).

  • Grid volumetric flux: facetVDotNormalSum{n,it} (sum of (\vec{w}_\text{facet}\cdot\bar{\vec{S}})).

  • ALE usage: movingMesh::gridVolFlux is subtracted from the fluid volumetric flux to form the relative/ALE mass flux throughout the solver (e.g., netMassFlux = rho_f*(volFluxCorrected_c - gridVolFlux) in the transport equations).

For BDF2 (constant (\Delta t)), the solver forms second-order grid terms using stored history:

  • movingMesh::gridVolFlux{n,it} = 1.5*facetVDotNormalSum{n,it} - 0.5*facetVDotNormalSum{n} (after the first step).

  • The same idea is applied to area vectors inside the momentum-interpolation flux corrections by replacing the instantaneous face area vector with (1.5,\texttt{facetNormalSum}^{n+1/2}-0.5,\texttt{facetNormalSum}^{n-1/2}) in the relevant volFluxStar{n,it} rules.

21.6.1. Restart requirements#

For BDF2 restarts, src/movingMesh.loci reads:

  • restart/volPrev_hdf5.* to seed vol{n=-1}.

  • restart/facetNormalSum_hdf5.* and restart/facetVDotNormalSum_hdf5.* to seed the BDF2 history of the GCL facet sums.

If these facet-history files are missing on a BDF2 restart, the code aborts (because a second-order history is required).

21.7. Assumptions and Limitations#

  • Linear-in-time mesh motion per step: the derivation of Eq. (21.5.2.4) assumes vertex trajectories are linear over each (\Delta t) (constant vertex velocities over the step).

  • Facet-based representation: each polygon face is assumed to have a consistent node ordering and a well-defined face center so that the fan triangulation ((\vec{C},\vec{P}i,\vec{P}{i+1})) is valid.

  • BDF2 constant time step (current implementation): the BDF2 grid-flux and area-vector combinations in src/movingMesh.loci are noted in the source as “currently only good for constant timestep”; variable-(\Delta t) coefficients would need to be generalized.

21.8. References#