2.1 Existing unsteady CFD meshing methods

Historically, the simulation of unsteady aerodynamic problems has been restricted to existing techniques such as mesh motion, Chimera grids or immersed boundary methods amongst others, that allow volume meshes to accommodate surface mesh motions. These methodologies must often be used along with an arbitrary Lagrangian-Eulerian (ALE) formulation of the fluid equations (3) and come with limitations regarding the type of motions they can cope with. In general, mesh deformation techniques can only deal with small movements if a good quality mesh is to be retained after the deformation process. The fact that no cells can appear or disappear due to a fixed connectivity between cells at two consecutive time levels yields distorted cells with high aspect ratios when large body motions are involved. In such cases re-meshing, i.e. generating a completely new mesh from the geometry at the current time level, would be a suitable solution to the problem of low-quality meshes. However, this implies using an interpolation method to relate the fluid variables in the new mesh with those in the previous one, introducing a computational burden. Interpolation in this manner is not a trivial task and can be non-conservative.

When complex rotational parts or relative motions [Reference Dougherty, Benek, Steger, Dougherty, Benek, Steger and Research Center4] are involved, Chimera or overset grids become a more reasonable alternative to mesh motion techniques thanks to the use of a separate body-fitted mesh for each of the moving parts. In addition, simple yet powerful high-quality structured grids can be used around each of the moving parts, which translates into more efficient and faster fluid solvers and mesh generators [Reference Meakin and Meakin5, Reference Renzoni, D’Alascio, Kroll, Peshkin, Hounjet, Boniface, Vigevano, Allen, Badcock, Mottura, Schöll and Kokkalis6]. Boundary motions are very much simplified and only a rotation and/or translation of the existing grids is required before the intersection process happens again, hence saving computational effort. Despite these benefits, interpolation algorithms needed at the boundaries of two overlapping grids are usually costly and complex [Reference Renzoni, D’Alascio, Kroll, Peshkin, Hounjet, Boniface, Vigevano, Allen, Badcock, Mottura, Schöll and Kokkalis6], and can introduce numerical errors unless special care is taken to minimise them. They still cannot deal with arbitrary motions such as aeroelastic problems [Reference Pomin and Wagner7], where a mesh deformation technique is still required in addition, or situations involving topological changes with appearing/disappearing cells, which, once again, rely on interpolations of the solution.

As opposed to Chimera, sliding grid planes are based on grids whose boundaries fit together without any overlapping at all, and slide past each other when there exists a relative motion. An interpolation method must still be put into place in order to communicate flow variables at both sides of the interface [Reference Rumsey8, Reference Steijl and Barakos9]. A method for the study of helicopter rotor-fuselage interaction is proposed by Steijl et al. [Reference Steijl and Barakos9] proving its accuracy and efficiency, provided the mesh size is not too big, but performing poorly under parallel computations. Moreover, limitations regarding the allowable timestep are also important [Reference Rumsey8, Reference Fenwick and Allen10].

Immersed boundary methods [Reference Lai and Peskin11, Reference Peskin12, Reference Kidron, Mor-Yossef and Levy13] or Cartesian cut-cell grids [Reference Kidron, Mor-Yossef and Levy13, Reference Meinke, Schneiders, GÜnther and SchrÖder14] can also be a feasible alternative to deal with mesh deformation in unsteady aerodynamics. Immersed boundary methods often fail to ensure conservation of mass, momentum or energy in cells cut by a solid boundary and it is therefore not a popular technique across the aerospace industry where compressible aerodynamics demand a good quality representation of boundaries. Moreover, very thin boundary layers cannot be captured unless the cell count is high, meaning large and expensive simulations, or anisotropic refinement is used [Reference de Tullio, De Palma, Iaccarino, Pascazio and Napolitano15].

Meshless methods [Reference Katz16, Reference Oñate, Idelsohn, Zienkiewicz and Taylor17, Reference Monaghan and Gingold18, Reference Katz16, Reference Netuzhylov and Zilian19], in an attempt to circumvent the bottleneck of automated mesh generation for complex geometries with sharp edges [Reference Katz16], replace the traditionally used grid by a dense cloud of points based on which conservation laws can be discretised [Reference Ortega, Onate, Idelsohn and Flores20]. While connectivity information is inherently lost, more programming effort is needed compared to traditional mesh-based methods since there is still the need for finding neighbours residing within the domain of influence of each node [Reference Monaghan and Gingold18], which can be a time consuming task.

As outlined before, re-meshing can deal with arbitrary motion, preserving good-quality meshes at all times [Reference Formaggia, Peraire and Morgan21, Reference Murman, Aftosmis and Berger22, Reference Probert, Hassan, Morgan and Peraire23]. Compared to the above methods it is likely to be the most demanding approach in terms of computational effort. Being capable of dealing with structured grids, it is with unstructured meshes where the greater gains in efficiency are achieved given the ability to modify only certain regions of the domain. Nevertheless, it is always necessary to work out the connectivity relationship between cells at consecutive time levels and an appropriate interpolation method has to be derived in the event of topological changes such as appearing/disappearing cells, which may introduce numerical errors across the solution.

The spacetime framework explored here offers an alternative conservative simulation approach even with topological changes and variable real time-steps, and if appropriately implemented preserves time accuracy.

2.2 Development of spacetime method

The concept of considering meshes in time as well as space is an old one; it might easily be argued the approach is not a departure from existing techniques, but rather a return to the natural dimensional setting of the differential system, which is then immediately discretised. As such, there are numerous published examples of its application, although its application to the most relevant industrial cases is yet to be realised.

Early spacetime results came from Giles [Reference Giles24] in 1988. The imposition of periodic boundary conditions in turbomachinery flows is a complex task, especially when the rotor and stator have different pitch values (i.e. distance between blades). By inclining the computational time plane Giles circumvents this issue and transforms the Euler equations so that any stator-rotor pair can be treated with a pitch ratio of 1. At the same time, Hughes et al. [Reference Hughes and Hulbert25] apply a spacetime technique to classical elastodynamics problems via the use of a finite element approach with a discontinuous Galerkin (DG) formulation in the time direction. Lowrie et al. [Reference Lowrie, Roe and van Leer26] build on the previous work for a much more general problem involving hyperbolic conservation laws. Again, they use a discontinuous Galerkin formulation to create a higher order scheme within the spacetime framework. However, this results in a computationally expensive method compared to other conventional approaches. Moreover, their work implicitly assumes some regularity on the structure of the grids used, hence invalidating a more general boundary motion approach. Thompson et al. [Reference Thompson and He27] and Ray [Reference Ray28] use a DG formulation to give their own interpretation of the spacetime method. In particular, Thompson et al. [Reference Thompson and He27] present an adaptive spacetime technique that allows refinement and coarsening of the grid and which they define as robust. This robustness comes at a sacrifice of the general applicability of the method since they retain orthogonal planes in the time direction leaving the time integration fully decoupled from the space integration.

Tsuei et al. [Reference Tsuei, Liou and Yu29] successfully apply the spacetime method developed earlier by Chang [Reference Chang30, Reference Chang, Wang and Chow31, Reference Chang32] at NASA to blade row interaction problems in turbomachinery flows. This space-time conservation element and solution element method (CE/SE) is able to predict unsteady flows without any previous assumptions imposed on the solver, by simply considering fluxes both in space and time. They argue that the scope of this new method is large and that a wide range of applications can benefit from it, however their work is focused on turbomachinery applications and no attempts are made towards a more general and arbitrary motion. Perhaps the most general implementation of the spacetime method are the works by Hixon [Reference Hixon33, Reference Hixon34] and Golubev et al. [Reference Golubev, Mankbadi and Hixon35, Reference Golubev and Rohde36]. Their method allows for a variable timestep size across the fluid domain and no decoupling is made between temporal and spatial integrations, allowing for higher-order schemes to be used in the time integration. Zwart et al. [Reference Zwart, Raithby and Raw37] apply the spacetime formulation to the solution of a breaking dam, although their implementation lacks a general spacetime mesh with varying timestep sizes across the spatial domain. Similarly, Van der Ven [Reference van der Ven38] applies a conservative adaptive multigrid algorithm under the spacetime framework to investigate an oscillating two-dimensional aerofoil, demonstrating the potential of the method. Rendall et al. [Reference Rendall, Allen and Power2, Reference Rendall and Allen39, Reference Rendall and Allen40, Reference Power, Rendall and Allen41] use a general formulation of the spacetime method and show its ability to simulate complex moving geometries in one and two-dimensional unsteady cases. They observe a slight difference in the pressure distribution with respect to a conventional solver and explain it with information propagated backwards in time as a consequence of the central-difference stencil used, as noted previously.

Parallel to the work presented in this paper, Wang et al. [Reference Wang and Persson42] develop a high-order discontinuous Galerkin spacetime formulation for fully unstructured meshes. In fact, like Thompson et al. [Reference Thompson and He27], they still retain orthogonal planes in the time direction but they manage to generate a fully unstructured spacetime mesh between two consecutive time slabs, being able to effectively simulate complex motions and topology changes. They successfully solve the compressible Navier-Stokes equations for a spinning cross, a pair of NACA 0012 aerofoils pitching in tandem and a spoiler case.

Recent work by Nishikawa and Padway has developed a Jacobian free Newton-Krylov implicit approach for solving the full, discrete unsteady system [Reference Nishikawa and Padway43], using a generalised conjugate gradient method for the linear solve, and extended this to viscous cases [Reference Nishikawa and Padway44] including a shedding cylinder in cross-flow and a boundary layer. In their work this was successfully coupled to mesh adaptation on a tetrahedral grid for both a conventional second order scheme, and a low-cost third order method. Results were also presented for a vortex, moving cylinder and shock problems, illustrating excellent results. It is important to appreciate that use of mesh refinement implies an adaptive, locally varying and solution dependent real timestep, alongside the usual spatial adaptivity. This is a particularly strong point of using meshes in space and time; they are not only able to handle any motion type, but also permit large changes in space and time refinement within a single framework. Indeed, mesh adaptation in four dimensions has been been explored by Caplan et al. [Reference Caplan, Haimes, Darmofal and Galbraith45] and indicates the feasibility of such an approach for three-dimensional unsteady problems.

4.1 Two-dimensional upwind formulation of inviscid fluxes

At any inclined face in a spacetime grid the total inviscid flux may be split into space and time contributions. In a two-dimensional case there will be three terms: two contained within the spatial plane $\left( x,y \right)$ and one along the time direction t. Equation (10) can be cast to

(14) \begin{equation}\begin{split}\widetilde{\mathbf{R}} = \sum_\textrm{faces} \widetilde{\mathbf{F}} \widetilde{\mathbf{n}} S & = \sum_\textrm{faces} \left( \mathbf{f}_t n_t + \mathbf{f}_x n_x + \mathbf{f}_y n_y\right) S = \\& = \sum_\textrm{faces} \mathbf{f}_t S_t + \mathbf{f}_x S_x + \mathbf{f}_y S_y,\end{split}\end{equation}

where $\mathbf{f}_t$ is the first column vector in the two-dimensional version of Equation (11), which represents the time fluxes

(15) \begin{equation}\mathbf{f}_t=\begin{Bmatrix}\rho \\\rho u \\\rho v \\\rho E\end{Bmatrix},\end{equation}

and $\mathbf{f}_x$ and $\mathbf{f}_y$ are the second and third column vectors in the two-dimensional version of Equation (11), which correspond to the space fluxes

(16) \begin{equation}\mathbf{f}_x=\begin{Bmatrix}\rho u \\\rho u^2 + p \\\rho u v \\\left(\rho E + p \right) u\end{Bmatrix} \qquad\qquad \mathbf{f}_y=\begin{Bmatrix}\rho v \\\rho u v \\\rho v^2 + p \\\left(\rho E + p\right) v\end{Bmatrix}.\end{equation}

As implied by a characteristic analysis, a pseudo-velocity that is constant and equal to one can be defined in the time direction, i.e. $u^*=\dfrac{dt}{dt^*}=1$ . This means that, for a physically meaningful solution, information in time is always convected from the upstream (or previous in time) cell which corresponds to that with the smallest t coordinate. Therefore time fluxes must be calculated using the primitive variables evaluated only at the upstream side of the face in Equation (15), i.e.

(17) \begin{equation}\mathbf{f}_t \left( \mathbf{W}^+ \right) = \begin{Bmatrix}\rho^+ \\\rho^+ u^+ \\\rho^+ v^+ \\\rho^+ E^+\end{Bmatrix},\end{equation}

where superscript $+$ denotes upstream and the column vector $\mathbf{W}^+$ of primitive variables is evaluated at the upstream side of the face

(18) \begin{equation}\mathbf{W}^+=\begin{Bmatrix}\rho^+\\u^+\\v^+\\p^+\end{Bmatrix}.\end{equation}

The sign of $n_t = \widetilde{\mathbf{n}}\cdot\mathbf{e}_t$ at each face may be used to discern between upstream and downstream cells in time (or past and future cells).

For fluxes in space information may be convected from both sides, upstream and downstream, if the flow is subsonic. This is consistent with the flux-vector splitting method developed by Van Leer [Reference Van Leer53] which works out the contribution of each side of the face to the total flux.

Since only the component of the velocity normal to the face will give non-zero fluxes in the momentum equation, the local normal Mach number $M_n=u_n/\sqrt{\gamma p/\rho}$ is used. There are two possible cases. If the normal flow is supersonic $\left( \left| M_n \right| \geq 1 \right)$ fluxes are given by the properties at the upstream side only

(19) \begin{equation}\mathbf{f}_n\left(\mathbf{W}^+_n\right) = \begin{Bmatrix}\rho^+ u^+_n \\[2mm]\rho^+ \left( u^+_n \right)^2 + p^+ \\[2mm]\rho^+ u^+_n u^+_{t_1} \\[2mm]\left(\rho^+ E^+ + p^+\right)u^+_n\end{Bmatrix}.\end{equation}

In the case of subsonic flow $\left( \left| M_n \right| < 1 \right)$ fluxes are formed by contributions from both sides, upstream ( $+$ ) and downstream (−)

(20) \begin{equation}\mathbf{f}_n\left(\mathbf{W}_n\right) = \mathbf{f}^+_n\left(\mathbf{W}^+_n\right)+\mathbf{f}^-_n\left(\mathbf{W}^-_n\right),\end{equation}

where the normal flux functions $\mathbf{f}^+_n$ and $\mathbf{f}^-_n$ are given [Reference Van Leer53] in Equation (21), as follows

(21) \begin{equation}\mathbf{f}^{\pm}_n\left(\mathbf{W}^{\pm}_n\right) =\pm \dfrac{\rho^{\pm} c^{\pm}}{4} \left(M_n^{\pm} \pm 1\right)^2 \begin{Bmatrix}1 \\[3mm]\dfrac{c^{\pm}}{\gamma}\left[\left(\gamma-1\right)M_n^{\pm} \pm 2\right] \\[3mm]u_{t_1}^{\pm} \\[3mm]u_{t_2}^{\pm} \\[3mm]\dfrac{\left(c^{\pm}\right)^2\left[ \left(\gamma-1\right)M_n^{\pm} \pm 2 \right]^2}{2\left(\gamma^2-1\right)} + \dfrac{\left(u_{t_1}^{\pm}\right)^2+\left(u_{t_2}^{\pm}\right)^2}{2}\end{Bmatrix},\end{equation}

and the column vectors of normal primitive variables, $\mathbf{W}^+_n$ and $\mathbf{W}^-_n$ , at the upstream and downstream sides of the face, respectively, are

(22) \begin{equation}\mathbf{W}^+_n=\begin{Bmatrix}\rho^+\\u_n^+\\u_{t_1}^+\\u_{t_2}^+\\p^+\end{Bmatrix}\qquad \qquad \mathbf{W}^-_n=\begin{Bmatrix}\rho^-\\u_n^-\\u_{t_1}^-\\u_{t_2}^-\\p^-\end{Bmatrix}.\end{equation}

Notice here the difference between functions $\mathbf{f}_n\left(\mathbf{W}_n^\pm\right)$ and $\mathbf{f}_n^\pm\left(\mathbf{W}_n^\pm\right)$ . Also, $u_n$ , $u_{t_1}$ and $u_{t_2}$ , given by Equation (23), are the components of the real velocity $\mathbf{v}=\left\lbrace 0,u,v \right\rbrace^T$ projected onto a spacetime coordinate system defined locally at each face such that axis n is normal to the face and the other two, $t_1$ and $t_2$ , are tangential. Depending on the face orientation the new coordinates may not be purely spatial or purely temporal but a combination of spatial and temporal coordinates. In other words, the projection of the flow velocity $\mathbf{v}$ , strictly defined in space $\left( x,y \right)$ , onto the local coordinate system yields components in spacetime, as follows,

(23) \begin{equation}\begin{Bmatrix}u_n \\u_{t_1} \\u_{t_2} \end{Bmatrix} = \begin{bmatrix} \\ & & & P^T & & & \\ \\ \end{bmatrix}\begin{Bmatrix}0 \\u \\v \end{Bmatrix} = \begin{Bmatrix}un_x+vn_y \\\dfrac{vn_x-un_y}{\sqrt{n_x^2+n_y^2}} \\[7mm]\dfrac{n_t\left(un_x+vn_y\right)}{\sqrt{n_x^2+n_y^2}}\end{Bmatrix},\end{equation}

where transformation matrix $P\;:\;\left(t,x,y\right) \mapsto \left(n,t_1,t_2\right)$ , which maps the global space+time coordinate system to the spacetime coordinate system locally defined at each face, is given later by Equation (27). Note here the intended distinction between a space+time ( $\mathbb{R}^2 \cup \mathbb{R}$ ) and a spacetime ( $\mathbb{R}^3$ ) coordinate system. The former is a concatenation of space and time into one single frame while still keeping space and time coordinates separate. The latter, however, constitutes a coupling between space and time coordinates such that an increment in one of the coordinates yields increments both in space and time. At this point, many different local spacetime ( $\mathbb{R}^3$ ) coordinate systems may be defined at each face in the mesh. However, since the purpose of this projection is the application of Van Leer’s flux-splitting method to the spatial fluxes, the normal vector $\widetilde{\mathbf{n}}$ to the face is taken as the first local direction

(24) \begin{equation}\mathbf{e}_1 = \widetilde{\mathbf{n}} = \begin{Bmatrix}n_t\\[4mm]n_x\\[4mm]n_y\end{Bmatrix}.\end{equation}

For the second and third directions there exists an infinite number of possible vectors perpendicular to $\mathbf{e}_1$ and between each other, i.e. such that $\mathbf{e}_i \cdot \mathbf{e}_j = 0$ for $i \neq j$ , as required for the coordinate system to be orthogonal. Nevertheless, for the sake of simplicity, $\mathbf{e}_2$ is chosen such that it has a null component in time t. After normalisation ( $\vert\vert \mathbf{e}_2 \vert\vert = 1$ ) it can be written as follows

(25) \begin{equation}\mathbf{e}_2 = \mathbf{t}_1 = \begin{Bmatrix}0 \\\dfrac{-n_y}{\sqrt{n_x^2+n_y^2}}\\[7mm]\dfrac{n_x}{\sqrt{n_x^2+n_y^2}}\end{Bmatrix}.\end{equation}

The third direction may be defined such that the orientation of the new coordinate system is right-handed or positive, i.e.

(26) \begin{equation}\mathbf{e}_3 = \mathbf{e}_1 \times \mathbf{e}_2 \quad \Longrightarrow \quad \mathbf{e}_3 = \mathbf{t}_2 = \begin{Bmatrix}-\sqrt{n_x^2+n_y^2} \\\dfrac{n_xn_t}{\sqrt{n_x^2+n_y^2}}\\[7mm]\dfrac{n_yn_t}{\sqrt{n_x^2+n_y^2}}\end{Bmatrix},\end{equation}

where, again, a normalisation has been applied so that $\vert\vert \mathbf{e}_3 \vert\vert = 1$ . As in any coordinate system transformation column vectors in matrix P are each of the unit vectors of the new base $\left( \mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3 \right)$ written in terms of the old base’s $\left( \mathbf{e}_t,\mathbf{e}_x,\mathbf{e}_y \right)$ , i.e.

(27) \begin{equation}P = \begin{bmatrix} \begin{Bmatrix} \\ \\ \mathbf{e}_1 \\ \\ \\ \end{Bmatrix} \begin{Bmatrix}\\ \\ \mathbf{e}_2 \\ \\ \\ \end{Bmatrix} \begin{Bmatrix}\\ \\ \mathbf{e}_3 \\ \\ \\ \end{Bmatrix}\\ \end{bmatrix} = \begin{bmatrix}n_t\;\; & 0\;\; & -\sqrt{n_x^2+n_y^2} \\n_x\;\; & \dfrac{-n_y}{\sqrt{n_x^2+n_y^2}}\;\; & \dfrac{n_xn_t}{\sqrt{n_x^2+n_y^2}} \\[7mm]n_y\;\; & \dfrac{n_x}{\sqrt{n_x^2+n_y^2}}\;\; & \dfrac{n_yn_t}{\sqrt{n_x^2+n_y^2}}\end{bmatrix}.\end{equation}

The inverse transformation matrix yields

(28) \begin{equation}P^{-1}=P^T=\begin{bmatrix}n_t\;\; & n_x\;\; & n_y \\[3mm]0\;\; & \dfrac{-n_y}{\sqrt{n_x^2+n_y^2}}\;\; & \dfrac{n_x}{\sqrt{n_x^2+n_y^2}} \\[7mm]-\sqrt{n_x^2+n_y^2}\;\; & \dfrac{n_xn_t}{\sqrt{n_x^2+n_y^2}}\;\; & \dfrac{n_yn_t}{\sqrt{n_x^2+n_y^2}}\end{bmatrix}.\end{equation}

Space fluxes calculated through Equations (19)–(23) need to be projected back onto the global space coordinates, x and y. Second and third rows in matrix P, Equation (27), are used for this purpose in Equation (29)

(29) \begin{equation}\mathbf{f}_x n_x + \mathbf{f}_y n_y = \begin{bmatrix}1\;\; & 0\;\; & 0\;\; & 0\;\; & 0\\[2mm]0\;\; & n_x\;\; & \dfrac{-n_y}{\sqrt{n_x^2+n_y^2}}\;\; & \dfrac{n_xn_t}{\sqrt{n_x^2+n_y^2}}\;\; & 0\\0\;\; & n_y\;\; & \dfrac{n_x}{\sqrt{n_x^2+n_y^2}}\;\; & \dfrac{n_yn_t}{\sqrt{n_x^2+n_y^2}}\;\; & 0\\[2mm]0\;\; & 0\;\; & 0\;\; & 0\;\; & 1\end{bmatrix} \begin{Bmatrix}\rho u_n \\[2mm]\rho u_n^2 + p \\[2mm]\rho u_n u_{t_1} \\[2mm]\rho u_n u_{t_2} \\[2mm]\left(\rho E + p\right)u_n\end{Bmatrix},\end{equation}

which can be shown to yield

(30) \begin{equation}\begin{Bmatrix}\rho u \\\rho u^2 + p \\\rho u v \\\left(\rho E + p\right)u\end{Bmatrix} n_x + \begin{Bmatrix}\rho v \\\rho v u \\\rho v^2 + p \\\left(\rho E + p\right)v\end{Bmatrix} n_y.\end{equation}

4.2 Extrapolation to face values

In a central-difference scheme the value of the primitive variables at each face is worked out as an average of the values at the neighbouring cells. In an upwind scheme the direction of propagation of some quantities determine whether the values used at the cell interface are taken from the upstream or downstream neighbour cell. First-order methods use the cell-centred value of the neighbour cell as the value at the face. However, a second-order correction term may be added to the cell-centred values when extra accuracy is required. In order to avoid spurious oscillations resulting from these second-order correction terms, the so-called slope limiters $\varphi_t$ , $\varphi_x$ and $\varphi_y$ are introduced in the current two-dimensional spacetime formulation (the superscripts f and c denote the face and cell centres, respectively)

(31) \begin{equation}\mathbf{W}^\textrm{f} = \mathbf{W}^\textrm{c} + \underbrace{\nabla \mathbf{W}^\textrm{c} \mathbf{\Phi} \left( \textbf{q}^\textrm{f} - \textbf{q}^\textrm{c}\right).}_\textrm{second-order correction term}\end{equation}

Where gradients are found using a Green-Gauss integration. $\mathbf{q}$ is the column vector of spacetime coordinates

(32) \begin{equation}\mathbf{q} = \begin{Bmatrix}q_1\\q_2\\q_3\end{Bmatrix} = \begin{Bmatrix}t\\x\\y\end{Bmatrix},\end{equation}

the matrix of slope limiters $\mathbf{\Phi}$ is defined as

(33) \begin{equation}\mathbf{\Phi} = \begin{bmatrix}\varphi_t\;\;\;\;\; & 0\;\;\;\;\; & 0\\0\;\;\;\;\; & \varphi_x\;\;\;\;\; & 0\\0\;\;\;\;\; & 0\;\;\;\;\; & \varphi_y\end{bmatrix},\end{equation}

and the matrix of gradients of the primitive variables $\nabla \mathbf{W}^\textrm{c} = \dfrac{\partial \mathbf{W}^\textrm{c}}{\partial \mathbf{q}}$ is given by

(34) \begin{align}\nabla \mathbf{W}^\textrm{c} = \begin{bmatrix} \\[10pt] & & & & \dfrac{\partial W_i^\textrm{c}}{\partial q_j} & & & & \\ \\[10pt] \end{bmatrix} = \begin{bmatrix} \left. \begin{matrix} \\ \\ \dfrac{\partial \mathbf{W}^\textrm{c}}{\partial t} \\ \\ \\ \end{matrix} \right\vert & \left. \begin{matrix} \\ \\ \dfrac{\partial \mathbf{W}^\textrm{c}}{\partial x}\\ \\ \\ \end{matrix} \right\vert & \dfrac{\partial \mathbf{W}^\textrm{c}}{\partial y} \end{bmatrix} = \begin{bmatrix} \dfrac{\partial \rho}{\partial t} & \dfrac{\partial \rho}{\partial x} & \dfrac{\partial \rho}{\partial y} \\[4mm]\dfrac{\partial u}{\partial t} & \dfrac{\partial u}{\partial x} & \dfrac{\partial u}{\partial y}\\[4mm]\dfrac{\partial v}{\partial t} & \dfrac{\partial v}{\partial x} & \dfrac{\partial v}{\partial y}\\[4mm]\dfrac{\partial p}{\partial t} & \dfrac{\partial p}{\partial x} & \dfrac{\partial p}{\partial y} \end{bmatrix}.\end{align}

The limiters have been chosen to comply with TVD (Total Variation Diminishing) conditions [Reference Sohn54] and depend upon the changes of the fluid variables in the nearby of the cell. To account for these changes one can define the monitor at cell j as the ratio

(35) \begin{equation}r^j = \min\left\lbrace\dfrac{\nabla W^1}{\nabla W^j},\dots,\dfrac{\nabla W^{n_\textrm{ne}}}{\nabla W^j}\right\rbrace,\end{equation}

where $n_\textrm{ne}$ is the number of neighbouring cells. In the specific literature there are many available methods for the computation of the limiters. One of the most well-known slope limiters is due to Van Leer [Reference Van Leer53], Equation (36), and is the one used in this work.

(36) \begin{equation}\varphi^\textrm{VL} \left( r \right) = \dfrac{r+\left|r\right|}{1+\left|r\right|}.\end{equation}

It was further hypothesised that a useful step in the development of a spacetime framework would be taking advantage of a central-difference approach in space whilst still upwinding in time. The idea underpinning this hybrid (CSUT, namely Central-difference in Space, Upwind in Time) formulation would allow the strength and robustness of the JST dissipation scheme to be retained and, at the same time, achieve more time accurate solutions, comparable to those obtained through the upwind formulation, as a consequence of the time stencil. Moreover it will be shown using convergence residual plots for a number of initial test cases (for instance see Fig. 8 for the simple flap deployment) that the convergence of unsteady problems is improved with respect to that of a purely CD spacetime formulation.

As done in the upwind case, the spacetime flux at any inclined face in the mesh may be split into space and time contributions. Therefore, Equation (14) is still applicable in the current case where time fluxes can be worked out using Equations (17) and (18). For the space fluxes, however, the usual central-difference formulation is used in this case, so for example

(37) \begin{equation}\mathbf{f}_x\left(\mathbf{W^\textrm{f}}\right) = \mathbf{f}_x\left(\dfrac{\mathbf{W}^{+}+\mathbf{W}^{-}}{2}\right).\end{equation}

4.3 Discretisation of viscous and turbulence terms

In order to ensure numerical stability of the solution (55), a second order central-difference scheme was used for the viscous terms $\mathbf{F}_\textrm{vis}$ , regardless of the stencil chosen for the convective part of the Navier-Stokes equations. For the Spalart-Allmaras turbulence model, a second order central-difference scheme was also used for the diffusive terms, but with a first-order upwind scheme for the convective terms.