The basic field equation to be solved for the grid dynamics is as follows:
\begin{equation} \ddot{s} + 2\left(\zeta\omega - \zeta^\prime\omega^\prime \nabla^2_\xi \right)\dot{s} + \left(\omega^2 - \omega^{\prime 2} \nabla^2_\xi \right)s = 2\zeta\omega\dot{s}_s + \omega^2 s_s \end{equation}Our discretization of this equation will be cell-centered to match the discretization scheme of DPLR, which is where we ultimately expect to implement this algorithm. Furthermore, since we do not require high accuracy numerics or or strong conservation properties, we will use second order finite differences to approximate the spatial derivatives:
\begin{equation} \nabla^2_\xi s = s_{i-1,j} + s_{i+1,j} + s_{i,j-1} + s_{i,j+1} - 4s_{i,j} = -\left(4s_i - \sum_{j \in \mathcal{J}_i} s_j \right) \end{equation}Where $\set{J}_i$ is the set cells neighboring the $i$-th cell in the mesh. Applying the discretization, we arrive at the following set of coupled, ordinary differential equations:
\begin{equation} \ddot{s}_i + 2\left(\zeta\omega + 4\zeta^\prime\omega^\prime\right) \dot{s}_i + \left(\omega^2 + 4\omega^{\prime 2} \right) s_i - 2\zeta^\prime\omega^\prime \sum \dot{s}_j - \omega^{\prime 2}\sum s_j = 2\zeta\omega\dot{s}_{s,i} + \omega^2 s_{s,i} \end{equation}Let $\vec{S}_i = [ s_i, \dot{s}_i ]^T $. With this notation, we can re-write our second order differential equation as system of first order equations.
\begin{equation} \dot{\vec{S}}_i + \left(\mat{U} + \mat{K} + 4\mat{K}^\prime \right)\vec{S}_i - \sum \mat{K}^\prime \vec{S}_j = \mat{K}\vec{S}_s \end{equation}Where
\begin{align} \mat{U} & = \left[\begin{array}{cc} 0 & 1 \\ 0 & 0 \\ \end{array}\right] \\ \mat{K} & = \left[\begin{array}{cc} 0 & 0 \\ \omega^2 & 2\zeta\omega \\ \end{array}\right] \\ \mat{K}^\prime & = \left[\begin{array}{cc} 0 & 0 \\ \omega^{\prime 2} & 2\zeta^\prime\omega^\prime \\ \end{array}\right] \\ \end{align}DLPR uses implicit time integration to circumvent the CFL restrictions associated with grids that have very fine near-wall spacing. The integration engine uses the 1st-order accurate Backward Euler scheme with the option of modifying the residual to achieve 2nd order accuracy. Since our grid dynamics equation is linear with constant coefficients, deriving the implicit time advancement equation is trivial.
\begin{equation} \frac{1}{\Delta t} \Delta \vec{S}_i + \left( \mat{U} + \mat{K} + 4\mat{K}^\prime \right) \left( \vec{S}_i + \Delta\vec{S}_i \right) - \sum \mat{K}^\prime \left( \vec{S}_j + \Delta \vec{S}_j \right) = \mat{K} \left( \vec{S}_s + \Delta \vec{S}_s \right) \end{equation}\begin{equation} \left(\frac{\mat{I}}{\Delta t} + \mat{U} + \mat{K} + 4\mat{K}^\prime \right) \Delta \vec{S}_i - \sum \mat{K}^\prime \Delta \vec{S}_j - \mat{K} \Delta \vec{S}_s = \vec{F}_i\left(\vec{S},\vec{S}_s\right) \end{equation}Where
\begin{equation} \vec{F}_i\left(\vec{S},\vec{S}_s\right) = \mat{K}\vec{S}_s + \sum \mat{K}^\prime \vec{S}_j - \left(\mat{U} + \mat{K} - 4\mat{K}^\prime \right) \vec{S}_i \end{equation}In the case of a stationary shock where $\Delta\vec{S}_s = 0$, the equation above degenerates into a block-pentadiagonal linear system. However, while iterating to steady state, the shock will in fact move and change velocity depending on the local flow solution near the shock front. Therefore this is (one of) the terms whereby the linear solve for the grid dynamics couples into the linear solve for the flow dynamics.
The fluid dynamics of for the problem are governed by the Navier Stokes equations, which can be written in integral form for an arbitrarily deforming control volume as follows:
\begin{equation} \frac{d}{dt} \int_{V(t)} \vec{U} dV + \int_{\partial V(t)} \left( \mat{F} - \vec{U}\dot{\vec{x}}^T \right) d\vec{A} = 0 \end{equation}Where, the flux matrix $\mat{F}$ is
\begin{equation} \mat{F} = \left[ \vec{F}_z, \vec{F}_y, \vec{F}_z \right] \end{equation}A typical finite volume code will solve this conservation equation for each cell in the computational grid. Consider the $i$-th cell in a regular hexahedral mesh with the volume-average state $\vec{U}_i$. In this situation, the integral of the flux across the boundary of the control volume will be approximated using a midpoint quadrature on each of the six faces as follows:
\begin{equation} \frac{d}{dt}\vec{U}_iV_i + \sum_{j\in\set{J}_i} \left( \mat{F}_j - \tilde{\vec{U}}_j\dot{\vec{x}}_j^T \right) \vec{A}_j = 0 \end{equation}As with the grid dynamics, we will apply Backward Euler time advancement to the discretized fluid equations.
\begin{equation} \frac{1}{\Delta t} \left( \vec{U}_i^{n+1}V_i^{n+1} - \vec{U}_i^nV_i^n \right) + \sum_{j \in \set{J}_i} \left( \mat{F}^{n+1}_j - \tilde{\vec{U}}_j^{n+1} {\vec{x}}_j^{n+1,T} \right) \vec{A}_j^{n+1} = 0 \end{equation}Introducing a linear approximatation for all $n+1$ terms, we arrive at the following coupled time advancement equation for the flow state in each computational cells:
\begin{equation} \frac{1}{\Delta t} \left( \Delta\vec{U}_i V_i + \vec{U}_i\Delta V_i \right) + \sum_{j \in \set{J}_i} \left( \Delta\mat{F}_j - \Delta \left( \tilde{\vec{U}}_j {\dot{\vec{x}}_j^T} \right) \right) \vec{A}_j + \left( \mat{F}_j - \tilde{\vec{U}}_j {\dot{\vec{x}}_j^T} \right) \Delta\vec{A}_j = \vec{R}_i(\vec{U},\vec{x}) \end{equation}Where
\begin{equation} \vec{R}_i(\vec{U},\vec{x}) = -\sum_{j \in \set{J}_i} \left( \mat{F}_j - \tilde{\vec{U}}_j \dot{\vec{x}}_j^{T} \right) \vec{A}_j \end{equation}In the above equation, three terms appear in the summation over the faces. The first term represents the standard linearization of the flux matrix with respect to the state variables. This term is already implemented in DPLR and would be present in most CFD codes that use implicit time advancement. The second term arises due to changes in the grid velocity, while the third arises from changes in the shape of each grid cell.
In the general case, the flux matrix is a non-linear function of the reconstructed fluid state on either side of the cell face, $\vec{U}_j^L,\vec{U}_j^R$ and the reconstructed gradients, $\nabla\vec{U}_j^L,\nabla\vec{U}_j^R$. These reconstructions are in turn computed from a non-linear combination of the cell-averaged data in neighboring cells and may take into account the geometry of the grid. To begin, consider the variation in one column of the flux matrix. Via the chain rule, we have:
\begin{align} \Delta \vec{F} &= \frac{\partial \vec{F}}{\partial \vec{U}^L} \left[ \sum_{s \in \set{S}_L} \frac{\partial\vec{U}^L}{\partial \vec{U}_s} \Delta \vec{U}_s + \sum_{v \in \set{V}_L} \frac{\partial\vec{U}^L}{\partial \vec{x}_v} \Delta \vec{x}_v \right] \\ &+ \frac{\partial \vec{F}}{\partial \grad\vec{U}^L} \left[ \sum_{t \in \set{T}_L} \frac{\partial\grad\vec{U}^L}{\partial \vec{U}_t} \Delta \vec{U}_t + \sum_{w \in \set{W}_L} \frac{\partial\grad\vec{U}^L}{\partial \vec{x}_w} \Delta \vec{x}_w \right] \\ &+ \frac{\partial \vec{F}}{\partial \vec{U}^R} \left[ \sum_{s \in \set{S}_R} \frac{\partial\vec{U}^R}{\partial \vec{U}_s} \Delta \vec{U}_s + \sum_{v \in \set{V}_R} \frac{\partial\vec{U}^R}{\partial \vec{x}_v} \Delta \vec{x}_v \right] \\ &+ \frac{\partial \vec{F}}{\partial \grad\vec{U}^R} \left[ \sum_{t \in \set{T}_R} \frac{\partial\grad\vec{U}^R}{\partial \vec{U}_t} \Delta \vec{U}_t + \sum_{v \in \set{W}_R} \frac{\partial\grad\vec{U}^R}{\partial \vec{x}_w} \Delta \vec{x}_w \right] \\ \end{align}| $\;$Variable$\;$ | Description |
|---|---|
| $\set{S}_L, \set{S}_R$ | Set of state vectors used to reconstruct $\vec{U}^L,\vec{U}^R$ |
| $\set{V}_L, \set{V}_R$ | Set of grid vertices used to reconstruct $\vec{U}^L,\vec{U}^R$ |
| $\set{T}_L, \set{T}_R$ | Set of state vectors used to reconstruct $\grad\vec{U}^L,\grad\vec{U}^R$ |
| $\set{W}_L, \set{W}_R$ | Set of grid vertices used to reconstruct $\grad\vec{U}^L,\grad\vec{U}^R$ |
Typically, a single reconstruction scheme will be used for all states in $\vec{U}$; similarly, the same scheme will typically be used for all gradients in $\grad\vec{U}$. In this scenario, the matrices $\partial\vec{U}^L/\partial\vec{U}_s$, $\partial\vec{U}^L/\partial\vec{x}_v$, etc. will collapse to scalar values. However, that's still a lot of terms that must be evaluated to exactly linearize the flux matrix.
Having noted that a precise linearization is possible for an arbirary flux matrix and reconstruction scheme, I'm going to focus on linearizing the Euler (inviscid) flux matrix assuming 1st order accurate state reconstruction. That's all I will need for my initial proof-of-concept implementation. With this approach $\vec{U}^L=\vec{U}_i$ and $\vec{U}^R=\vec{U}_{j}$ at the $j$-th face; there is no functional dependence on $\grad\vec{U}^L,\grad\vec{U}^R$. Therefore,
\begin{equation} \Delta\vec{F}_j = \frac{\partial \vec{F}}{\partial \vec{U}^L} \Delta\vec{U}_i + \frac{\partial \vec{F}}{\partial \vec{U}^R} \Delta\vec{U}_j \end{equation}Next let us consider the linearization of the flux due to grid motion, $\vec{U}_j\dot{\vec{x}}_j$. To evalute this term we take the interface state to be an average of the left and right states:
\begin{equation} \tilde{\vec{U}}_j = \frac{1}{2}\left(\vec{U}^L + \vec{U}^R\right) = \frac{1}{2}\left(\vec{U}_i + \vec{U}_j\right) \end{equation}Next, to evaluate the grid velocity, we will assume that each quadrilateral face is defined using bi-linear interpolation of it's vertices. Assuming $\vec{t} \in [-1,1] \times [-1,1]$:
\begin{align} 4\vec{x}_j(\vec{t}) &= (1 - t_1)(1 - t_2)\;\vec{x}_{j1} \\ &+ (1 + t_1)(1 - t_2)\;\vec{x}_{j2} \\ &+ (1 + t_1)(1 + t_2)\;\vec{x}_{j3} \\ &+ (1 - t_1)(1 + t_2)\;\vec{x}_{j4} \end{align}Taking the derviative of this expression and evaluating at the midpoint of the face yields the following expression for $\dot{\vec{x}}_j$:
\begin{equation} \dot{\vec{x}}_j = \frac{1}{4} \sum_{k=1}^4 \dot{\vec{x}}_{jk} \end{equation}The linearization of the product $\vec{U}_j\dot{\vec{x}}_j$ may then be written as:
\begin{equation} \Delta\left(\tilde{\vec{U}}_j\dot{\vec{x}}_j\right) = \frac{1}{2}\left( \Delta\vec{U}_i + \Delta\vec{U}_j \right)\dot{\vec{x}}_j + \frac{1}{4}\tilde{\vec{U}_j} \sum_{k=1}^4 \Delta\dot{\vec{x}}_{jk} \end{equation}Formally, the directed face area, $\vec{A}_j$, is is defined as
\begin{equation} \vec{A}_j = \int_{\set{F}_j} \vec{n} dA = \int_{-1}^1 \int_{-1}^1 \vec{n} |J| dt_1 dt_2 = \int_{-1}^1 \int_{-1}^1 \frac{\partial \vec{x}_j}{\partial t_1} \times \frac{\partial \vec{x}_j}{\partial t_2} dt_1 dt_2 \end{equation}Differentiating the coordinate transformation in the previous section
\begin{equation} 4\frac{\partial \vec{x}_j}{\partial t_1} = (1-t_2)(\vec{x}_{j2} - \vec{x}_{j1}) + (1+t_2)(\vec{x}_{j3} - \vec{x}_{j4}) \end{equation}\begin{equation} 4\frac{\partial \vec{x}_j}{\partial t_2} = (1-t_1)(\vec{x}_{j4} - \vec{x}_{j1}) + (1+t_1)(\vec{x}_{j3} - \vec{x}_{j2}) \end{equation}With some algebra, the cross product of these two vectors can be written as:
\begin{align} 16 \frac{\partial \vec{x}_j}{\partial t_1} \times \frac{\partial \vec{x}_j}{\partial t_2} &= (1 - t_1)(1 - t_2)\;\vec{n}_1 + (1 + t_1)(1 - t_2)\;\vec{n}_2 \\ &+ (1 + t_1)(1 + t_2)\;\vec{n}_3 + (1 - t_1)(1 + t_2)\;\vec{n}_4 \end{align}Which is to say that the normal at any point on the face is obtained by interpolating the normal vector at each vertex. Integrating this function over the face yields:
\begin{equation} 4\vec{A}_j = \sum_{k=1}^4 \vec{n}_k = \sum_{k=1}^4 (\vec{x}_{k+1} - \vec{x}_k) \times (\vec{x}_{k-1} - \vec{x}_k) \end{equation}Consider the variation of the sum with respect to vertex $\vec{x}_v$:
\begin{equation} 4 \Delta\vec{A}_j = \Delta\vec{x}_{ v } \times (\vec{x}_{v-2} - \vec{x}_{v-1}) - \Delta\vec{x}_{ v } \times (\vec{x}_{v-1} - \vec{x}_{ v }) - (\vec{x}_{v+1} - \vec{x}_{ v }) \times \Delta\vec{x}_{ v } + (\vec{x}_{v+2} - \vec{x}_{v+1}) \times \Delta\vec{x}_{ v } \end{equation}This expression can be simplified by exchanging the order of the first two cross products and collecting terms. If one then sums the linear sensitivity for all four vertices, the following linearization is obtained:
\begin{equation} 2 \Delta\vec{A}_j = \sum_{k=1}^4 (\vec{x}_{k-1} - \vec{x}_{k+1}) \times \Delta\vec{x}_k = \sum_{k=1}^4 \set{C}(\vec{x}_{k-1} - \vec{x}_{k+1}) \Delta\vec{x}_k \end{equation}In the above equation, the $\set{C}(\vec{x})$ operator returns the cross product equivilent matrix of $\vec{x}$, i.e. a matrix with the elements of $\vec{x}$ distributed such that $\set{C}(\vec{x})\,\vec{y} = \vec{x}\times\vec{y}$.
The final term that we need to linearize is the change in volume. To do this, condsider the application of our discrete flow equation to a static, spatially uniform flow field. In this cases, the flux matrix is identical at all quadrature points and $\vec{U}_i = \tilde{\vec{U}}_j$. Since the integral of a constant flux over a closed surface is zero, we obtain the Geometric Conservation Law, or GCL:
\begin{equation} \frac{d}{dt}V_i = \sum_{j\in\set{J}_i} \dot{\vec{x}}_j^T \vec{A}_j \end{equation}There is a considerable amount of published literature on the GCL and how best to deal with it. This literature is generally focused on time-accurate calculations with the goal of updating the mesh in such a way as to retain second-order time accuracy. However, second order accuracy is not required in our application, so we will approximate the change in volume using the standard Backward Euler update:
\begin{align} \frac{1}{\Delta t}\Delta V_i &= \sum_{j\in\set{J}_i} \dot{\vec{x}}_j^T \vec{A}_j + \Delta\dot{\vec{x}}_j^T \vec{A}_j + \dot{\vec{x}}_j^T \Delta\vec{A}_j \\ &= \sum_{j\in\set{J}_i} \left[ \dot{\vec{x}}_j^T \vec{A}_j + \frac{1}{4} \sum_{k=1}^4 \Delta\dot{\vec{x}}_{jk}^T \vec{A}_j + 2\dot{\vec{x}}_j^T \set{C}(\vec{x}_{j,k-1} - \vec{x}_{j,k+1}) \Delta\vec{x}_{jk} \right] \end{align}In the preceeding section, all linearizations were written in terms of the change in position and velocity of the vertices defining a computational cell. However, the motion of the cell vertices are derived from the motion of the $\xi_3=N_s$ isosurface and the clustering function used to distribute grid nodes along each $\zeta_1,\zeta_3$-gridline.
Let $s_j,\;j \in 1...4$ be the curvilinear position of the $\xi_3=N_s$ isosurface at the midpoint of the four cell centers surrounding the $i$-th gridline of the mesh (these are the unknowns of the grid motion equation above). Furthermore, define
Note that even though arc length is measured inward from the outer boundary, grid points are ordered such that grid points radiate outward from the surface of the vehicle.
The distribution function, $\set{D}$, takes as input the grid point location in computational space, $\xi_{3,k}$, and return the position of the point along the line. This position is normalized by the distance between the vehicle surface and the $\xi_3$-isosurface such that $\set{D} = 0$ yields a point on the vehicle surface:
\begin{equation} s_k = s_{MAX} - \left(s_{MAX} - s_{MIN} \right) \set{D}(\xi_{3,k}) \end{equation}We have a few choices for the grid distribution function: uniform, geometric, hyperbolic tangent, etc. Since our proof of concept will assume an inviscid flow field, we can utilize a uniform grid distribution as shown below.
\begin{equation} \set{D}(\xi_3) = \frac{\xi_3-1}{N_s} \end{equation}The nodal coodinates for each point in the computational grid are then computed via arc-length interpolation of the background grid:
\begin{equation} \vec{x}_{ik} = \set{X}_i(s_k) = \set{X}_i\left( s_{MAX} - \left(s_{MAX} - s_{MIN} \right) \set{D}(\xi_{3,k}) \right) \end{equation}This relationship provides the basis for coupling the linear update equations for the grid and flowfield. Linearizing our relationship, we find that:
\begin{equation} \Delta \vec{x}_{ik} = \frac{\partial \set{X}_i}{\partial s} \frac{\partial s}{\partial s_{MIN}} \Delta s_{MIN} = \vec{t}_k \set{D}_k \Delta s_{MIN} \end{equation}Where
| Symbol | Definition |
|---|---|
| $\vec{t}_{ik}$ | Tangent vector for gridline $i$ at $\vec{x}_k$ |
| $\set{D}_k$ | Value of the grid distribution function at $\vec{x}_k$ |
Since the linear update equation for the $\xi_3=N_s$ isosurface is written in terms of cell-centered values, we must interpolate these updates from the cell centers to the grid lines to compute $\Delta s_{MIN}$. A simple linear interpolation in the computational domain is sufficient:
\begin{equation} \Delta s_{MIN} = \frac{1}{4}\sum_j \Delta s_j \end{equation}Combining the above equations:
\begin{equation} \Delta \vec{x}_{ik} = \vec{t}_k \set{D}_k \sum_j \Delta s_j \end{equation}Consider again the mapping from computational space to physical space via interpolation of the background grid.
\begin{equation} \vec{x}_{ik} = \set{X}_i\left( s_{MAX} - \left(s_{MAX} - s_{MIN} \right) \set{D}(\xi_{3,k}) \right) \end{equation}This can be differntiated to find the velocity of each grid point in terms of the rate of change of the $s_{MIN}$:
\begin{align} \dot{\vec{x}}_{ik} &= \frac{\partial \set{X}_i}{\partial s} \frac{\partial s}{\partial s_{MIN}} \dot{s}_{MIN} \\ &= \vec{t}_k \set{D}_k \dot{s}_{MIN} \\ &= \vec{t}_k \set{D}_k \sum_j \dot{s}_j \end{align}When linearizing this relationship, we must take into consideration the fact that our gridlines may be curved, and therefore $\vec{t}_k$ will change as $s_{MIN}$ changes. Therefore, we have:
\begin{equation} \Delta \dot{\vec{x}}_{ik} = \vec{\kappa}_k\set{D}_k \dot{s}_{MIN} \sum_j \Delta{s}_j + \vec{t}_k \set{D}_k \sum_j \Delta \dot{s}_j \end{equation}With the paramter $\vec{\kappa}_k$ representing the curvature of the gridlines at the location of the $k$-th gridpoint.
The final details required to implement this method relate to the manner in which we interpolate the background grid. We have already established that we will be performing arc-length interpolation and will require tangent and curvature vectors for each gridline.
First, let us establish how to compute arc length along a given gridline. If we assume that the background grid is defined using a mesh of discrete points (as opposed to some analytical mapping), that means that a spline of some kind is required to create a continuous, differentiable mapping from $\zeta_3$ to physical space. There are several strategies that can be employed here:
Catmul-Rom Cubic Spline: This $C^1$ continuous spline is easy to compute on the fly using finite differences to estimate the nodal derivatives. However, since we are not $C^2$ continuous, jumps in the curvature parameter will occur (which may or may not be a problem). Furthermore, there if there are any kinks in the background grid, undesireable local extrema may occur in the interpolant.
Akima Cubic Spline: A more elaborate $C^1$ spline that is robust to kinks and discontinuities in the grid. Like Catmul-Rom, can still be computed on-the-fly and will exhibit jumps in the curvature.
Classic Cubic Spline: $C^2$ continuous, so guarantees smooth curvature metrics. However, computing the derivatives requires a tridiagonal solve so this method is not amenable to on-the-fly computation. Rather, we would precompute the slopes and store them, doubling the required storage. However, the background mesh can be fairly coarse in the wall normal direction, so this may not be a deal breaker.
Catmul-Rom Quintic Spline: $C^2$ continuous by construction and easy to compute on the fly. If memory is a concern, this is a nice option, but would be even more prone to local extrema on non-smooth grids. An interesting option, but would complicate the lookup process.
Generally speaking, I think any of these method will work. The Akima