A spacetime formulation for unsteady aerodynamics with geometry and topology changes

Abstract

A spacetime formulation is presented to solve unsteady aerodynamic problems involving large deformation or topological change such as store separation, slat and flap deployment or spoiler deflection. This technique avoids complex CFD meshing methods, such as Chimera, by the use of a finite-volume approach both in space and time, and permits a locally varying real timestep. The use of a central-difference scheme in the time direction can yield non-physical transient solutions as a consequence of information travelling backwards in time. Therefore, an upwind formulation is provided and validated against one-dimensional and two-dimensional test cases. A hybrid formulation (central in space, upwind in time) is also given and unsteady cases are computed for a spoiler and spoiler/flap deployment, with all three formulations compared, demonstrating that the use of an upwind time stencil yields more representative physical solutions and improves the rate of convergence.

Nomenclature

Symbol

Description

a

speed of sound

$\alpha$

angle-of-attack

b

aerofoil’s semi-chord

E

total energy of fluid

$\mathbf{e}_i$

unit vector along i direction

$\gamma$

heat capacity ratio of fluid

$\mathbf{F}$

matrix of fluxes

$\mathbf{f}_i$

column vector of fluxes in i

k

reduced frequency

$\varphi_i$

slope limiter along i direction

$M_n$

normal Mach number

$\mathbf{n}$

normal vector

$n_i$

normal vector component in i

$n_\textrm{ne}$

number of neighbouring cells

$\nabla$

divergence operator

p

pressure of fluid

$\mathbf{R}$

column vector of residuals

$r_j$

monitor of slope limiter at cell j

$\rho$

density of fluid

S

surface area

$S_i$

projection of surface area normal to i

t

time coordinate

$t^*$

pseudo-time coordinate

$\tau,t_1,t_2$

tangential directions

$\theta$

angle of deflection

$\omega$

angular velocity, angular frequency

JST

Jameson-Schmidt-Turkel

CD

central-difference

CSUT

hybrid stencil (central in space, upwind in time)

VL

upwind stencil (Van Leer)

$\mathbf{U}$

column vector of conserved variables

$u_n$

normal velocity

u, v

fluid velocity in Cartesian coordinates

V

control volume

$\mathbf{W}$

column vector of primitive variables

x, y

spatial coordinates

Subscripts and superscripts

Symbol

Description

$+$

upstream

−

downstream

$\widetilde{\phi}$

generic variable $\phi$ in spacetime

c

cell centre

f

face centre

1.0 Introduction

Applications of computational fluid dynamics can be found in such diverse fields as medicine, aerospace, marine engineering or the oil and gas industry. Although steady state fluids problems can sometimes be computationally expensive, there exist mature techniques and numerical methods to solve them accurately. In comparison, unsteady problems with complex motions or topology changes can be very intricate and are still an active area of research. The study of interaction between helicopter rotor-blades and a fuselage constitutes a clear example of the complexity of unsteady aerodynamics. Other common complex problems include store separation, flap and spoiler deployment, or the transient process that takes place within an internal combustion engine when valves open and close. Finding accurate and efficient solutions to these problems using the most common CFD methods remains limited by the capability of existing mesh generation/deformation techniques and interpolation algorithms. A different meshing technology needs to be used depending on the problem under consideration which, inevitably, limits the ability to automate simulations and slows down the design cycle for industrial application.

Integration across the four-dimensional space-time domain is required to obtain unsteady solutions of the three-dimensional Navier-Stokes equations. Conventional methods have traditionally decoupled this process into two consecutive and different steps: a finite-volume (or finite-element) integration in a three-dimensional space and a finite-difference integration in time. The novelty introduced by the spacetime method implemented here is that the integrations in space and time are treated similarly through the use of spacetime finite-volumes, effectively solving unsteady problems of dimension N as steady problems of dimension $N+1$ . This implies modifying the fluid equations of motion through the divergence theorem to remove temporal derivatives, which are replaced by temporal fluxes. The fully conservative nature in space and time facilitates the solution of problems with geometry and topology changes, where cells can appear and disappear between consecutive points in time. A further strong advantage is that since the cells in spacetime can vary in size, the local real timestep becomes easily variable, permitting high and lower temporal resolution in regions (or at moments) of interest.

However, coupling time and space, and effectively solving both as one, leads to non-physical behaviour when using a central-difference scheme since information is propagated backwards in time as a consequence of the temporal stencil being used. The solution at a certain time level may be affected by the solution at later times and this is not physically correct (note that slight influence may still exist where an upwind scheme uses a gradient found from a central stencil). This may appear as shown in Figs. 9 and 10 of Flamarique et al. [1] where there is an oscillation in pressure before instantaneous motion of a flap, or as in Fig. 8 of Rendall et al. [2] where there is a small phase shift in the unsteady pressures for a pitching aerofoil. Although for a periodic problem the phase shift may be acceptable, for instantaneous motion it is not physical.

Understanding the direction in which characteristics of the hyperbolic problem transport disturbances across the spacetime domain is essential to obtain accurate and meaningful solutions, so this work makes a novel spacetime comparison between (i) a dissipation-based central scheme (ii) an upwind method and (iii) a hybrid dissipation-based central scheme in space, but using an upwind stencil in time. This comparison is particularly important in the time direction where waves are unidirectional. An upwind method is seen to resolve this issue, even if gradients are still found using a full stencil in all directions. Inviscid solutions are also presented showing the novel capability of a spacetime approach specifically for aerospace applications, including landing and spoiler deployment cases, as well as application to a viscous pitching aerofoil.

The goal of the work presented here is therefore to expand the range of cases demonstrated for spacetime solution, including those with topological change and viscous effects, and to explore the implications of different upwind and central schemes.

2.0 Background

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 [4] 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 [5, 6]. 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 [6], and can introduce numerical errors unless special care is taken to minimise them. They still cannot deal with arbitrary motions such as aeroelastic problems [7], 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 [8, 9]. A method for the study of helicopter rotor-fuselage interaction is proposed by Steijl et al. [9] 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 [8, 10].

Immersed boundary methods [11, 12, 13] or Cartesian cut-cell grids [13, 14] 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 [15].

Meshless methods [16, 17, 18, 16, 19], in an attempt to circumvent the bottleneck of automated mesh generation for complex geometries with sharp edges [16], replace the traditionally used grid by a dense cloud of points based on which conservation laws can be discretised [20]. 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 [18], which can be a time consuming task.

As outlined before, re-meshing can deal with arbitrary motion, preserving good-quality meshes at all times [21, 22, 23]. 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 [24] 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. [25] 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. [26] 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. [27] and Ray [28] use a DG formulation to give their own interpretation of the spacetime method. In particular, Thompson et al. [27] 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. [29] successfully apply the spacetime method developed earlier by Chang [30, 31, 32] 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 [33, 34] and Golubev et al. [35, 36]. 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. [37] 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 [38] 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. [2, 39, 40, 41] 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. [42] develop a high-order discontinuous Galerkin spacetime formulation for fully unstructured meshes. In fact, like Thompson et al. [27], 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 [43], using a generalised conjugate gradient method for the linear solve, and extended this to viscous cases [44] 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. [45] and indicates the feasibility of such an approach for three-dimensional unsteady problems.

3.0 Formulation

The numerical solution of the two-dimensional Navier-Stokes equations for viscous compressible flows requires integration in both space and time, yielding

(1) \begin{equation} \int_{t_{0}}^{t_F}\int_{V}\left(\frac{\partial\textbf{U}}{\partial t} + \nabla \cdot \textbf{F}_{\textrm{inv}} - \nabla \cdot \textbf{F}_{\textrm{vis}} \right) d\Omega dt = \textbf{0,}\end{equation}

where $d\Omega$ is the differential volume, V the integration volume and $t_{0}$ and $t_{F}$ the start and finish times and where the vector of conserved variables $\mathbf{U}$ is

(2) \begin{equation} \mathbf{U}= \begin{Bmatrix} \rho \\ \rho u \\ \rho v \\ \rho E \end{Bmatrix},\end{equation}

where $\rho$ is density, u,v are the components of velocity and E is specific energy. The matrix of inviscid fluxes $\mathbf{F}_{\textrm{inv}}$ is

(3) \begin{equation}\mathbf{F}_{\textrm{inv}} =\begin{bmatrix}\rho u\;\;\; & \rho v \\\rho u^2 + p\;\;\; & \rho u v \\\rho u v\;\;\; & \rho v^2 + p \\\left(\rho E + p\right)u\;\;\; & \left(\rho E + p\right)v\end{bmatrix},\end{equation}

and the matrix of viscous fluxes $\mathbf{F}_{\textrm{vis}}$ is

(4) \begin{equation} \mathbf{F}_{\textrm{vis}} = \begin{bmatrix} 0 & 0 \\ \sigma_{xx} & \sigma_{xy} \\ \sigma_{xy} & \sigma_{yy} \\ u\sigma_{xx} + v\sigma_{xy} + q_x\;\;\;\; & u\sigma_{xy} + v\sigma_{yy} + q_y \end{bmatrix},\end{equation}

where $\sigma_{...}$ denotes Favre-averaged shear stress (including turbulence effects) on xx, xy and yy planes and q is the heat flux. The shear stresses are

(5) \begin{equation}\sigma_{ij} = \left(\mu + \mu_t\right) \left( 2S_{ij} - \dfrac{2}{3}\dfrac{\partial u_i}{\partial x_j}\delta_{ij} \right),\end{equation}

where $\mu$ is the dynamic viscosity, $\mu_t$ is the eddy viscosity (computed through the negative Spalart-Allmaras turbulence model), $\delta_{ij}$ is the Kronecker-Delta and $S_{ij}$ is the strain rate tensor defined as

(6) \begin{equation}S_{ij} = \dfrac{1}{2}\left(\dfrac{\partial u_j}{\partial x_i}+\dfrac{\partial u_i}{\partial x_j}\right).\end{equation}

The heat flux is given by

(7) \begin{equation}q_j = \left(\frac{\mu c_p}{Pr} + \frac{\mu_t c_p}{Pr_t}\right)\dfrac{\partial T}{\partial x_j},\end{equation}

where $c_p$ is the specific heat capacity at constant pressure, Pr is the Prandtl number and T is the temperature. The negative Spalart-Allmaras turbulence model [46, 47] is solved along with Equation (1) to compute the eddy viscosity $\mu_t$ .

The key aspect of the spacetime method is the treatment of the time integration identically to the space integration through the use of four-dimensional spacetime finite-volumes. Within this framework any unsteady problem of dimension N can be effectively solved as another steady problem of dimension $N+1$ . In the current two-dimensional problem, a three-dimensional divergence operator can be defined as (the tilde $ \widetilde{} $ over variables and operators denotes spacetime)

(8) \begin{equation} \widetilde{\nabla} = \frac{\partial}{\partial t}\mathbf{e}_t + \frac{\partial}{\partial x}\mathbf{e}_x + \frac{\partial}{\partial y}\mathbf{e}_y, \end{equation}

which leads to the spacetime formulation of the Euler Equations (1) as

(9) \begin{equation} \int_{\widetilde{V}}{\widetilde{\nabla} \cdot \widetilde{\mathbf{F}} \, d\widetilde{\Omega}} = \mathbf{0,}\end{equation}

where $\widetilde{V}$ and $d\widetilde{\Omega}$ are now spacetime volumes.

One can express this as a closed surface integral via the divergence theorem as

(10) \begin{equation} \widetilde{\mathbf{R}} := \oint_{\partial \widetilde{V}}{\widetilde{\mathbf{F}} \widetilde{\mathbf{n}} \, d\widetilde{S}} = \mathbf{0,}\end{equation}

where $\mathbf{R}$ are the residuals. The matrix of spacetime fluxes $\widetilde{\mathbf{F}}$ is

(11) \begin{equation} \widetilde{\mathbf{F}} = \begin{bmatrix} \left. \begin{matrix} \\ \\ \mathbf{U} \\ \\ \\ \end{matrix} \right\vert \mathbf{F}_{\textrm{inv}} - \mathbf{F}_{\textrm{vis}} \end{bmatrix},\end{equation}

and the spacetime unit normal vector $\widetilde{\mathbf{n}}$ (which has a component now in time as well as the spatial directions) is

(12) \begin{equation} \widetilde{\mathbf{n}} = \begin{Bmatrix} n_t \\ n_x \\ n_y \end{Bmatrix}.\end{equation}

Similarly, the negative Spalart-Allmaras turbulence model has been re-written in a spacetime formulation (see [48] for more details) and integrated along with Equation (10).

Equation (10) can be regarded as the integration of the two-dimensional unsteady Navier-Stokes equations across a theoretical three-dimensional space, and can be solved by iterating until residuals converge to zero, as

(13) \begin{equation} \frac{\partial \widetilde{V} \mathbf{U}}{\partial t^*} + \widetilde{\mathbf{R}} = \mathbf{0,}\end{equation}

where $t^{*}$ is pseudo-time.

Note that implicit schemes can also be used [43]; in general, the common families of solution methods may all be applied. Unlike more frequently applied methods, the use of a finite-volume approach for the discretisation in time, as well as in space, automatically ensures conservation of mass, momentum and energy and, more importantly, allows the use of a variable real timestep across the spatial domain without causing a non-physical behaviour of the solution. Notice the potential gain in efficiency over conventional time-stepping techniques due to the fact that a bigger timestep can be used in areas of freestream flow, far away from the perturbations, while still retaining sufficiently small timesteps in areas where rapid changes occur. In terms of solution accuracy the spacetime method also brings the possibility of incorporating higher-order schemes often used for spatial discretisations into the temporal dimension.

The spacetime framework works well with any arbitrary motion, from big boundary displacements, like a helicopter rotor blade, through to geometric topological changes such as a store separation or a slotted flap deflection. There is no need for further modifications to the solver in any of the former cases, mainly as a consequence of the finite-volume approach in time. All the information related to boundary motions is implicitly given by the spacetime mesh. An example of this is depicted in Fig. 1 where the pitching movement of a NACA 0012 aerofoil is given by a twisted wing in which the span-wise direction represents the time. Moreover, no additional connectivity relationship is required between cells at different time levels* since this is implicitly accomplished by the spacetime mesh itself, therefore allowing for appearing/disappearing cells without the need for interpolation.

Figure 1. Example of spacetime mesh: pitching NACA 0012 aerofoil.

Along with the development of a spacetime framework there is a need for new grid generation techniques. Being able to generate unstructured grids in the time direction brings the possibility to refine the timestep size in some areas of the domain while keeping a coarse one in others where temporal resolution is not required. Currently, it is possible to use available 3D grid generators to create 2D unsteady meshes well-suited for the spacetime framework. However, there is no available technology to automatically generate a truly unstructured four-dimensional grid to be used in 3D+t problems. Although not yet mature, Behr [49] introduces a simple meshing technique that allows unstructured grids to be created, not only for 2D+t problems but also 3D+t. Likewise, Ungor et al. [50] and, some time later, Abedi et al. [51] have put some efforts towards the development of 2D+t grids by a mesh-marching technique. Mesh adaptation in four dimensions has been been demonstrated by Caplan et al. [45], indicating the feasibility of the approach.

4.0 Numerical implementation

Using a dual time-stepping method [52], backward differences are implicitly being used for the time derivatives. This means that the solution at each time level depends only upon the solution calculated at previous time levels. Within the spacetime framework, however, there exists the possibility to use a central-difference scheme [2, 39, 41, 40] resulting in information being propagated backwards in time [2]. Namely, whatever happens in the future affects the solution in the past, which violates a principle of causality. A spacetime version of the upwind flux-splitting method proposed by Van Leer [53] is implemented in this work to address this, with the two-dimensional formulation given below.

The methods compared in this work therefore comprise:

1. 1. Central scheme in space and time, with JST dissipation. This is referred to as CD — central difference

2. 2. Central scheme in space, but with an upwind stencil in time and still applying JST dissipation. This is referred to as CSUT — cental in space, upwind in time

3. 3. An upwind in space and time method

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 [53] 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 [53] 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 [54] 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 [53], 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.

5.0 Results

5.1 Periodic semi-infinite piston

As a starting point in the validation of the upwind formulation of the spacetime method a periodic semi-infinite piston was tested. This is a simple one-dimensional test case and the fact that there exists an analytical solution, Equation (40), makes it a very suitable correlation. A periodicity condition was applied in time, i.e. left and right vertical boundaries on the spacetime mesh in Fig. 2 were connected. The top boundary was modelled as a moving solid wall and for the bottom one, non-reflecting boundary conditions were used. The motion is sinusoidal about the position $x_{0}$

(38) \begin{equation}x\left(t\right) - x_0 = \dfrac{\Delta L}{2} \cos\left( \dfrac{2\pi}{T} t \right),\end{equation}

and the reduced frequency

(39) \begin{equation}k=\dfrac{\pi \Delta L}{T a_\infty},\end{equation}

was $0.016$ . In Equations (38) and (39), $\Delta L$ and T are the amplitude and the period of the piston’s motion, respectively, and $a_\infty$ is the speed of sound at initial conditions. This setup allows the results to be compared with piston theory at the moving wall [56]

(40) \begin{equation}\frac{p}{p_\infty} = \left( 1+\frac{\gamma-1}{2}\frac{u_w}{a_\infty} \right)^{\frac{2\gamma}{\gamma-1}},\end{equation}

where $\gamma$ is the ratio of specific heats and $u_{w}$ is the velocity of the wall.

Figure 2. One-dimensional periodic semi-infinite piston: spacetime mesh (top) and pressure contours (bottom). Note that the amplitude of the upper wall motion is small.

Pressure contours for both the central-difference and upwind schemes are depicted in Fig. 2. Also, the pressure at the moving wall is compared against theoretical results over one whole oscillation in Fig. 3. These non-dimensional results correspond to a motion of amplitude $\Delta L=10.41\;\textrm{cm}$ at $1,000\;\textrm{rpm}$ with sea-level ISA atmosphere conditions, i.e. $\rho_\infty = 1.225\;\textrm{kg}\cdot\textrm{m}^{-3}$ and $p_\infty = 101325\;\textrm{Pa}$ . The value of the heat capacity ratio used was $\gamma=1.403$ and the maximum piston velocity at each cycle $V_\textrm{max}=5.45\;\textrm{m}\cdot\textrm{s}^{-1}$ . Results for both central-difference and upwind are in good agreement with piston theory. No noticeable differences appear between the central-difference and upwind stencils, which can be explained by the periodicity of the problem.

Although information can only travel forwards in time, in periodic problems information may seem to go also backwards in time thus justifying the use of a central-difference time stencil. Bearing in mind that this is just an illusion, the explanation relies on the fact that, at each cycle, the problem is influenced by any previous temporal state, hence later stages of the previous cycle (ahead in physical time domain) determine the solution at the earliest stages of the current cycle (behind in physical time domain). Of course, the real physical influence in time is only ever forwards, and this is merely a consequence of using a periodic condition for expediency.

Figure 3. Comparison of central-difference and upwind results to piston theory.

5.2 Piston with sharp change of direction

To validate the upwind formulation of the spacetime method in a non-periodic case, in which the solution at any time level can only be influenced by the solution at previous time levels, the one-dimensional piston given by the spacetime mesh in Fig. 4 was computed. Initially, up to $t=0.4$ , the piston travels downwards at a constant speed, compressing the gas inside the chamber. In contrast with the problem defined in section 5.1, the bottom boundary was set to a solid wall, leading to some wave reflections. At time $t=0.4$ the piston inverts its velocity and from this point onwards it moves upwards at a constant speed, expanding the gas inside the chamber. The aim of this configuration is to analyse whether the upwind stencil used solves the issue of pressure waves propagating backwards in time. As shown in Fig. 5, the upwind formulation improves considerably the prediction of sudden and fast movements in time when compared to a central-difference formulation. Upwind schemes commonly avoid spatial oscillations, but the important point is that they can also avoid them in time for these types of problem; the typical oscillatory behaviour near shocks for CD solvers is observed here in the time direction when a sudden change in the movement of a boundary occurs. The upwind formulation avoids these oscillations, successfully yielding a much smoother solution. The difference in the quality of the solutions is again explained by the fact that pressure waves always travel forwards in time. Figure 4 depicts pressure contours for this configuration.

Figure 4. Spacetime mesh and $C_p$ contours (for upwind in space and time case) for the one-dimensional finite piston with sharp movement.

Figure 5. Comparison between central-difference and upwind in space and time results for piston with sharp movement.

5.3 Periodic Pitching aerofoil

The first two-dimensional problem presented in this paper is a pitching NACA 0012 aerofoil simulation. The spacetime mesh was constructed from a two-dimensional structured mesh by stacking up planes in the time direction, as shown above in Fig. 1. The first and last planes were connected to achieve the periodic boundary condition and an oncoming flow set with a freestream velocity $M_\infty=0.85$ . The aerofoil oscillates sinusoidally at a reduced frequency of $k=\dfrac{\omega c}{2U_\infty}=0.0814$ about its quarter chord point according to

(41) \begin{equation}\alpha=\alpha_0+\Delta \alpha \sin \omega t,\end{equation}

where $\alpha_0=0\;\textrm{deg}$ and $\Delta \alpha = 2.51\;\textrm{deg}$ . $C_p$ distribution plots at $\omega t=0$ , $\omega t=\dfrac{2\pi}{3}$ and $\omega t=\dfrac{4\pi}{3}$ have been depicted in Fig. 6. The CSUT (Central-difference in Space, Upwind in Time) and CD (central-difference) formulations yield very similar $C_p$ distributions, likely down to the periodicity of the problem (note that the only difference between CSUT and CD comes from the time stencil). On the other hand, the upwind formulation gives a slightly different solution, probably due to upwinding in space. Given the periodicity of the problem all three methods converge at comparable rates, as depicted in Fig. 7, although the upwind method still produces a slightly improved convergence history.

Figure 6. $C_p$ distributions for pitching NACA 0012 at $M_\infty=0.85$ for $\omega t=0,\dfrac{2\pi}{3},\dfrac{4\pi}{3}$ .

Figure 7. Convergence residuals for pitching NACA 0012 at $M_\infty=0.85.$

5.4 Simple flap

Using radial basis functions to deform the two-dimensional mesh and stacking up planes in the time direction, as done in the case of the pitching aerofoil (section 5.3), a spacetime mesh was created to simulate the deflection of a simple flap on a NACA-0012 (Fig. 8). Initially the aerofoil is flying with an angle-of-attack $\alpha=0\;\textrm{deg}$ at $M_\infty=0.7$ . After some time the flap deflects an angle $\Delta \theta = 13.5\;\textrm{deg}$ at a reduced frequency $k=\dfrac{\omega b}{U_\infty}=0.375$ , where $\omega$ is the angular velocity, b is the semi-chord and $U_\infty$ is the freestream velocity. An unsteady solution for the transient process was sought and all three formulations were used to solve it. $C_p$ contour plots and distributions along the chord at various time levels/slices are depicted in Figs 9 and 10, respectively. The hybrid CSUT formulation avoids the non-physical behaviour of CD in the transient part (time slices 1, 2 and 3) by the use of a more realistic time stencil and matches the solution of the upwind scheme better. As concluded from slice 4 in Fig. 10, the CSUT steady state solution resembles the central-difference scheme more closely since time fluxes have a negligible impact on it and only space fluxes, which are worked out using central-differences in both cases, have an effect on the steady state solution. Moreover, both the upwind and CSUT solutions converge much faster than CD due to no information moving backwards in time, as demonstrated by convergence residuals in Fig. 8. This represents a significant improvement and demonstrates the applicability of the method to transient problems.

Figure 8. Simple flap deflection on a NACA 0012 aerofoil flying at $M_\infty=0.7$ : spacetime mesh (left) and convergence residuals (right).

Figure 9. $C_p$ contours for simple flap deflection at $M_\infty=0.7$ : (a) central-difference (left), (b) upwind (centre), and (c) central-difference in space, upwind in time (right).

Figure 10. $C_p$ distributions at four different time slices corresponding to those in Fig. 9 for simple flap deflection at $M_\infty=0.7.$

5.5 Spoiler

The mesh generation of two-dimensional problems like the pitching aerofoil (section 5.3), or even the simple flap deflection (section 5.4), in spacetime may be done using a three-dimensional structured mesh generator or simply by stacking up two-dimensional meshes. However, in the case of more complex and arbitrary boundary motions the use of unstructured meshes is crucial. Figure 11 shows the unstructured spacetime mesh used in the solution of the current spoiler deployment case. Initially a NACA 0012 aerofoil at an angle of incidence of $\alpha=0\;\textrm{deg}$ is immersed in a flow at $M_\infty=0.25$ . At some point, a spoiler deflects up to an angle of $\theta=45\;\textrm{deg}$ at a reduced frequency $k=\dfrac{\omega b}{U_\infty}=0.262$ , and the whole transient process is successfully captured. As in previous sections, pressure contour plots are depicted in Fig. 12 for several time levels and convergence residuals are given in Fig. 13. The upwind simulation converges much faster than the central-difference counterpart, as would be expected from the fact that it uses a more realistic time stencil. The convergence history of the CSUT formulation lies in between, improving the convergence of the central-difference but, at the same time, producing a more physically relevant solution, closer to that of the upwind formulation.

Figure 11. Spacetime mesh for spoiler deployment.

Figure 12. $C_p$ contours for spoiler deployment at $M_\infty=0.25$ : (a) central-difference (left), (b) upwind (centre) and (c) central-difference in space, upwind in time (right).

Figure 13. Convergence residuals: spoiler deployment at $M_\infty=0.25$ (top) and landing case at $M_\infty=0.15$ (bottom).

5.6 Simplified landing case

One of the main benefits of the spacetime formulation is its versatility and the fact that it is capable of handling very complex boundary motions with relative ease. Problems like slat and flap deployment can be solved without the need for further modifications to the solver nor the implementation of intricate interpolation methods that connect cells between different time levels. With conventional time-marching techniques, defining flow properties at positions where there was no fluid domain at a previous time level can sometimes be problematic. Moreover, this relies heavily on the accuracy of the interpolation method used. The use of a finite-volume method, conservative by nature, in space and time simplifies the problem considerably.

A simplified version of all the motions that a wing undergoes during approach and landing was modelled, i.e. a slat and flap deployment on an aerofoil flying at $M_\infty=0.15$ followed by an increase of its angle of incidence and a spoiler deployment which, in turn, decreases the incidence, all of which happens while approaching to the ground. Figure 14 shows the spacetime mesh used to represent the geometry for this problem and define the motions involved in it. As mentioned above, the solver can be left unchanged speeding up the overall simulation process, from meshing through to running the CFD code. Pressure contour plots with streamlines have been depicted in Figs 15 and 16 for all three formulations (upwind, central-difference and CSUT) and the history of convergence residuals is plotted in Fig. 13. In order to understand what the streamlines represent it is important to realise that the reference frame chosen for this simulation is not fixed to the aerofoil or the ground; it moves with the aerofoil in the horizontal direction but remains fixed in the vertical direction, i.e. null vertical velocity.

Figure 14. Spacetime mesh for landing case, i.e. aerofoil with a slat and flap deployment followed by an increase of its angle of incidence and a spoiler deployment which, in turn, decreases the incidence, all of which happens while approaching to the ground.

Figure 15. $C_p$ contours and streamlines for landing case at $M_\infty=0.15$ : (a) central-difference (left), (b) upwind (centre) and (c) central-difference in space, upwind in time (right). Continues in Figure 16.

Figure 16. Continuation from Figure 15. $C_p$ contours and streamlines for landing case at $M_\infty=0.15$ : (a) central-difference (left), (b) upwind (centre) and (c) central-difference in space, upwind in time (right).

Unexpectedly, the rate of convergence of the hybrid solver is not much better (or at least the difference is negligible) than that of the central-difference solver, as can be observed in Fig. 13. A possible explanation for this is that the gradients of the fluid properties are worked out, as usual, using a central-difference approach, despite the use of an upwind stencil in the time direction. Time fluxes are therefore slightly affected by future events. Although this is a minor contribution it may become more noticeable when the solution changes rapidly.

5.7 Viscous simulation of AGARD R-702(3E3) Case 3

Unlike previous inviscid cases, here the full Navier-Stokes model is used in spacetime, by including the viscous terms in addition to the Euler equations. The mesh was constructed from a two-dimensional structured mesh by stacking up grid planes in the time direction, as shown in Fig. 17. An O-grid of size $201 \times 60$ has been used to generate the spacetime meshes for the inviscid problems with 150 physical time-steps. Similarly, a C-grid of identical size, $201 \times 60$ , has been used with 100 physical time-steps in the case of viscous problems. Thus, the size of the physical time-step is computed as $\Delta t = \dfrac{T}{150}$ in the former case and $\Delta t = \dfrac{T}{100}$ in the latter, where $T=\dfrac{\pi c}{kU_\infty}$ is the period. In order to ensure a proper resolution of the boundary layer the first grid line normal to the wall is at a distance $\sim 10^{-5}$ , where the chord of the aerofoil is $c=1$ . This ensures $y^+ \sim \mathcal{O}\left(1\right)$ .

Figure 17. Example of spacetime geometry for pitching NACA 0012. Meshes are constructed by stacking up two-dimensional grid planes in the t direction. An O-grid is used in the inviscid case (left) and a C-grid is used in the viscous case (right).

The aerofoil follows a pitching motion about its quarter chord described by equation (41) as in the example of section 5.3 above. In this case, however, the aerofoil is submerged in a freestream flow at Mach $M_\infty = 0.6$ and $Re_\infty=4.8\times 10^6$ . The motion is defined by a mean angle of incidence $\alpha_0 = 2.44\text{ deg}$ , an amplitude $\Delta \alpha = 4.89\text{ deg}$ and a reduced frequency $k=0.0810$ . The value of $\Delta \alpha$ , which effectively describes the motion along with the value of the circular frequency $\omega$ , is implicit in the definition of the spacetime geometry whereas the value of $\alpha_0$ is just the mean angle of incidence which can be (and has been) given as a parameter to the solver directly.

Radial basis functions (RBFs) are used to deform the two-dimensional mesh at each $t=$ constant plane after a geometry transformation. The problem considered here is periodic, hence the first and last planes are connected to achieve the periodic boundary condition. Since the spacetime framework is conservative, both in space and time, periodic problems like this are particularly well suited because the solution can be said to have converged to the final solution once the $\ell^2$ -norms of the residuals have dropped beyond a certain threshold (notice that the residuals in spacetime represent the change in the solution throughout the whole period), provided that the numerical scheme is stable and convergent.

Both Euler and RANS solutions are computed. A central-difference scheme and an upwind method based on Van Leer fluxes are used in both cases. Moreover, a Roe-based spacetime solver is also used for comparison in the inviscid case. The Spalart-Allmaras model is solved with a first-order upwind discretisation and pseudo-time marching, alongside the remaining flow equations.

$C_p$ distributions at different phase angles $\omega t$ are depicted and compared with experimental data extracted from ref. [57] in Figs 18 and 19. Results via the spacetime method correlate well, especially in the inviscid cases. The upwind solutions match the experimental data more closely when compared to the CD formulation, with the exception of two phase angles, namely $\omega t = 59.85\text{ deg}$ and $\omega t = 264.81\text{ deg}$ . At all other times the CD solver seems to over-predict pressure slightly, which is particularly noticeable at the leading edge. These results can probably be explained by the observed phase lead of the central-difference scheme due to information propagating backwards in time, as mentioned earlier. Viscous spacetime solutions seem to under-predict the pressure coefficient in this case too, especially at high angles of attack where, perhaps, the turbulent boundary layer of the Spalart-Allmaras model delays, or even avoids, separation. Moreover, it is important to bear in mind that a larger physical time-step than in the Euler solutions has been used (due to availability of computational resources), hence a lower temporal-accuracy is obtained; this was achieved by using larger cells in the mesh in the time direction. It is interesting to notice the oscillatory solution of the viscous CD solution at $\omega t = 135.51\text{ deg}$ , typical of central-difference solvers around shock waves. In this case, however, this is a transient effect coming from the integration in the time direction. This behaviour is similar to that observed in the non-periodic simple flap problem, presented and explained in example 5.4. In Fig. 20 plots of the locus of the pitching moment coefficient, $C_m$ , and normal forces, $C_N$ , are also compared with experimental data. Viscous solutions yield a better prediction of normal forces in this case. Finally, the moment coefficient $C_m$ is predicted well by the RANS-SA CD solution whereas the viscous Van Leer and all three inviscid solutions deviate from the experimental data. The periodic nature of this test case favours the otherwise non-physical central time stencil of the CD solver in this case.

Figure 18. $C_p$ distribution plots for pitching NACA 0012 with amplitude $\Delta \alpha = 2.44$ deg at $M_\infty=0.6$ , $\alpha_0 = 4.86\text{ deg}$ and $k=\frac{\omega c}{2U_\infty}=0.0810$ . Comparison of CFD results via spacetime solver against experimental data from AGARD R-702(3E3) Case 3. Continues in Figure 19.

Figure 19. Continuation from Figure 18. $C_p$ distribution plots for pitching NACA 0012 with amplitude $\Delta \alpha = 2.44$ deg at $M_\infty=0.6$ , $\alpha_0 = 4.86\text{ deg}$ and $k=\frac{\omega c}{2U_\infty}=0.0810$ . Comparison of CFD results via spacetime solver against experimental data from AGARD R-702(3E3) Case 3.

Figure 20. $C_N$ (left) and $C_m$ (right) coefficients for pitching NACA 0012 with amplitude $\Delta \alpha = 2.44$ deg at $M_\infty=0.6$ , $\alpha_0 = 4.86\text{ deg}$ and $k=\frac{\omega c}{2U_\infty}=0.0810$ . Comparison of CFD results via spacetime solver against experimental data from AGARD R-702(3E3) Case 3.

6.0 Conclusions

The spacetime methodology for solving fluid equations for aerodynamic problems with large boundary motions and/or topology transformation has been presented and extended to more elaborate motions, as well as incorporating upwind methods in space and time. The applicability and versatility of the spacetime framework in unsteady aerodynamics has been successfully demonstrated. It allows the solution of an unsteady problem of dimension N as a steady problem of dimension $N + 1$ . Due to its conservative nature in space and, particularly, in time, spacetime is very well-suited for periodic problems where initial and final states are connected. Moreover, the use of a coupled finite-volume formulation in space and time allows different timestep sizes at different locations (through the use of higher dimensional unstructured meshes). Therefore to improve efficiency of simulations a bigger timestep can be used where changes are more gradual, and smaller timesteps where changes in the fluid properties happen faster.

Industrial applications could benefit substantially from the use of the spacetime framework due to its potential for highly automated CFD simulations which would, in turn, speed up the design cycle. It has been shown that the solver can be left unchanged throughout the whole range of problems described in this work. There is no need to implement intricate interpolation methods, which can lead to losses in accuracy, in order to connect cells between different time levels. Spacetime solvers cope well with appearing/disappearing cells, which brings great flexibility to the method, and shared- and/or distributed-memory parallelisation can be applied to the spacetime method. In particular, a shared-memory parallelisation based on OpenMP has been implemented and used throughout the simulations in this work. A small change is that coupling with structural solvers would require definition of an interaction time between the fluid and structure, over which conventional strong coupling cycles could be used. Conventional fluid-structure coupling usually takes place over one timestep, but since a spacetime method has no single timestep, the coupling would need to be at predefined increments along the time axis. Memory requirements for a spacetime method are in general above an ALE-based counterpart if it updates the entire time domain at each pseudo-time iteration; future developments could bring this to equivalence with existing methods by operating on separate time regions in the mesh in turn.

It has been demonstrated that the upwind formulation of the spacetime framework for unsteady problems yields more representative solutions than the central-difference counterpart. This is likely due to the use of a more realistic time stencil where information is not allowed to travel backwards in time. The central-difference and hybrid formulations have been found consistently more sensitive to the freestream velocity and rate of change of boundary movements than the upwind one. This is also noticeable in the rate of convergence of the simulations where the upwind solver outperforms the hybrid one, which, in general, converges faster than the central-difference counterpart.