Three-dimensional CFD simulations with large displacement of the geometries using a connectivity-change moving mesh approach

This paper deals with three-dimensional (3D) numerical simulations involving 3D moving geometries with large displacements on unstructured meshes. Such simulations are of great value to industry, but remain very time-consuming. A robust moving mesh algorithm coupling an elasticity-like mesh deformation solution and mesh optimizations was proposed in previous works, which removes the need for global remeshing when performing large displacements. The optimizations, and in particular generalized edge/face swapping, preserve the initial quality of the mesh throughout the simulation. We propose to integrate an Arbitrary Lagrangian Eulerian compressible flow solver into this process to demonstrate its capabilities in a full CFD computation context. This solver relies on a local enforcement of the discrete geometric conservation law to preserve the order of accuracy of the time integration. The displacement of the geometries is either imposed, or driven by fluid–structure interaction (FSI). In the latter case, the six degrees of freedom approach for rigid bodies is considered. Finally, several 3D imposed-motion and FSI examples are given to validate the proposed approach, both in academic and industrial configurations.


Introduction
Fluid-structure interaction (FSI) simulations are required for a wide variety of subjects, from the simulation of jellyfish [23] to the releasing of a missile [47]. The recent development of computing capacities has made it possible to run increasingly complex simulations where moving bodies interact with an ambient fluid in an unsteady way. However, engineers are still far from performing such simulations on a daily basis, largely due to the difficulty of handling the moving meshes induced by the moving geometries.
When the displacement of the geometry is small enough, slightly deforming the original mesh [9,19] can generally be acceptable. But when large deformation of the boundaries is considered, the mesh quickly becomes distorted and the numerical error due to this distortion quickly becomes too great, until the elements of the mesh finally become invalid, and the simulation has to be stopped. Specific strategies need to be developed to deal with large displacement moving boundary problems. In the case of FSI problems, another difficulty arises from the fact that the displacement of the boundaries is by definition an unknown, and the deformation of the mesh cannot be imposed a priori. Addressing the issue of the mesh movement cannot be separated from addressing the issue of the solver that will compute on these moving meshes. Depending on the strategy employed to deal with the movement, specific numerical methods must be designed to take into account the displacement of the mesh. The question this paper focuses on is: how can we efficiently move the mesh for large displacement 3D FSI simulations and what numerical schemes need to be associated to such a strategy?
Three main approaches to address the mesh movement problem can be found in the literature. The first approach consists in having a single body-fitted mesh [11,29], and moving it along with the moving boundaries. The mesh may thus undergo large deformation. The second approach is the Chimera (or overset) method [12], in which each moving body has its own body-fitted sub-mesh, and the sub-meshes move rigidly together with their body and can overlap one another. Finally the embedded boundary approach [14,39] uses meshes that are not body-fitted at all: the bodies are embedded in a fixed grid, and techniques such as level-sets are used to recover their moving boundaries. All three approaches have their own strengths and weaknesses. This paper is aligned with our previous works on anisotropic mesh adaptation, with the ultimate goal to use moving geometries in our adaptation framework [8]. For this reason, we focus on body-fitted approaches with one single mesh.
The first strategy to handle body-fitted moving meshes is simply to move the mesh for as long as possible, and remesh (i.e., generate a whole new mesh or a part of it) when the mesh quality becomes critical [11,29]. For each remeshing, the simulation has to be stopped, a new mesh must be generated, and the solution must be transferred to the new mesh. This approach can be efficient, especially when small displacements are considered and very few remeshings are necessary, because the solver and the meshing aspects are decoupled, and between two remeshings the simulation is fully ALE and free from interpolation errors. However, for larger displacements, the number of remeshings increases to prevent invalid elements from appearing, and this can both be costly and result in poor accuracy due to the solution transfer step. Hence a second strategy has been developed [16,20], based on the use of local remeshing operations, such as vertex insertion, vertex collapse, connectivity changes and vertex displacements, to preserve a good mesh quality throughout the simulation. The advantage of this method is that it maintains an acceptable mesh quality without needing to stop, remesh and resume the simulation. However, it requires fully dynamic mesh data structures that are permanently updated, which can lead to a loss of CPU efficiency, and the numerous mesh modifications can lead to a loss of accuracy.
Our approach tries to overcome to these drawbacks and is described in detail in [1]. It aims at moving meshes with large displacements of the geometry without ever having to remesh. [By remeshing, we mean stopping the simulation, generating a new mesh (the entire mesh or part of it) and interpolating the solution on the new mesh.] A limited set of mesh modifications are used to preserve the mesh quality throughout the simulation: only connectivity changes (edge swaps) and vertex displacement are performed. This is for several reasons. Notably, performing local mesh modifications within the solver is far simpler than remeshing globally, and connectivity changes can be relatively simply interpreted in terms of evanescent cells for purposes of Arbitrary Lagrangian Eulerian (ALE) numerical schemes. In many body-fitted moving mesh strategies, a lot of CPU time is dedicated to computing the displacement of the mesh, in order to make it follow the moving boundaries. Thanks to frequent mesh optimizations, the cost of this step is reduced by computing the mesh deformation for a large number of solver time steps (i.e., we do it only a few times during the simulation). It is important to note that this approach works best if vertex-centered solvers are considered, because the connectivity changes preserve the number of degrees of freedom.
Some studies try to impose a mesh motion that is directly adapted to the physical phenomena in question, using for instance either so-called Moving Mesh PDEs [32] or a Monge-Ampère equation [15]. However, interesting these approaches may be, they still seem to be time-consuming, especially in 3D, due to the solution of a non-linear equation, and it is unsure whether they can handle complex 3D geometries. Therefore we prefer to prescribe arbitrary movements to the mesh, and use our mesh adaptation framework [42] when necessary.
Regarding numerical solvers, we consider a classic framework for moving meshes: the Arbitrary Lagrangian Eulerian (ALE) framework, which is based on a formulation of the equations that takes into account an arbitrary movement of the vertices. This technique was introduced in the 1970s in [21,31,33]. Since then, so many developments have been made in that field that a complete list of them would not fit in this paper. However, one may in particular refer to [11,24,25,29,30,46,49], which mainly focus on improving temporal schemes for ALE simulations.
To our knowledge, very few examples of ALE solvers coupled with connectivity-change moving mesh techniques can be found in literature. In [36] a conservative interpolation is proposed to handle the swaps. In [28,51] an ALE formulation of the swap operator is built. However, these studies are limited to 2D. Driven by the requirements of industry, we are interested in designing a method that works in 3D. In this paper, a linear interpolation is carried out after each swap instead of using a specific ALE formulation, and we will evaluate the numerical error due to these swaps.
The goal of the present paper is to demonstrate that three-dimensional FSI simulations can be run efficiently by coupling an ALE solver to our connectivity-change moving mesh strategy. The first part of this paper focuses on recalling important aspects of the moving mesh algorithm. The second, third and fourth parts describe in detail the solver used. In the fifth part, some validation test cases are presented, and finally some examples of complex 3D ALE simulations are given and analyzed. In this paper, we only focus on rigid movements that are involved in rigid-body FSI.

Mesh-connectivity-change moving mesh strategy
To handle moving boundaries, we adopt a body-fitted approach, with a single mesh: the inner vertices of the mesh are moved following the moving boundaries to preserve the validity of mesh (i.e., to prevent the mesh from getting tangled). Our strategy involves two main parts: • Computing the mesh deformation Inner vertices are assigned a trajectory depending on the displacement of the boundaries, and thus a position for future time steps. • Optimizing the mesh The trajectories computed in the mesh deformation phase are corrected, and the connectivity of the mesh is modified to preserve the quality of the mesh.
This strategy, detailed below, has proven to be very powerful in 3D [1], since large displacement of complex geometries can be performed while preserving a good mesh quality without any global remeshing (i.e., without ever generating a whole new mesh).

Linear elasticity mesh deformation method
During the mesh deformation step, a displacement field is computed for the whole computational domain, given the displacement of its boundaries. Trajectories can thus be assigned to inner vertices, or in other words, positions at a future solver time step. Several techniques can be found to compute this displacement field: implicit or direct interpolation [13,44], or solving PDEs-the most common of which being Laplacian smoothing [40], a spring analogy [19] and a linear elasticity analogy [6]. It is this last method that we selected, due to its robustness in 3D [58]. The computational domain is assimilated to a soft elastic material, which is deformed by the displacement of its boundaries.
The inner vertices movement is obtained by solving an elasticity-like equation with a ℙ 1 finite element method (FEM): where and  are, respectively, the Cauchy stress and strain tensors, and is the Lagrangian displacement of the vertices. The Cauchy stress tensor follows Hooke's law for an isotropic homogeneous medium. Dirichlet boundary conditions are used and the displacement of vertices located on the domain boundary is strongly enforced in the linear system. The linear system is solved by a conjugate gradient algorithm coupled with an LU-SGS pre-conditioner. An advantage of elasticity-like methods is the opportunity they offer to adapt the local material properties of the mesh, especially its stiffness, according to the distortion and efforts borne by each element. In particular, the stiffness of the elements is increased for small elements, in order to limit their distortion. More details can be found in [1].

Improving mesh deformation algorithm efficiency
The computation of the mesh deformation-here the solution of a linear elasticity problem-is known to be an expensive part of dynamic mesh simulations, and the fact that it is usually performed at every solver time step makes it all the more so. We propose to combine several techniques to improve the time efficiency of this step. Some regions are rigidified, more specifically a few layers around tiny complex details of the moving bodies, with very small elements. They are moved with exactly the same rigid displacement as the corresponding body, thus avoiding very stiff elements in the elasticity matrix. On the other hand, the elasticity can be solved only on a reduced region, if the domain is big compared to the displacement. A coarse mesh can also be used to solve the elasticity problem, the displacement of the vertices then being interpolated on the computational mesh.
The major improvement we proposed is to reduce the number of mesh deformation computations: the elasticity problem is solved for a large time frame of length Δt instead of doing it at each solver time step t . While there is a risk of a less effective mesh displacement solution, it is a worthwhile strategy if our methodology is able to handle large displacements while preserving the mesh quality. Solving the previously described mesh deformation problem once for large time frame could be problematic in the case of: (1) curved trajectories of the boundary vertices and (2) accelerating bodies. To enhance the mesh deformation prescription, accelerated-velocity curved, i.e., high-order, vertex trajectories are computed.
The paths of inner vertices can be improved if a constant acceleration is provided to each vertex in addition to its speed, which results in an accelerated and curved trajectory.

3
During time frame [t, t + Δt] , the position and the velocity of a vertex are updated as follows: Prescribing a velocity vector and an acceleration vector to each vertex requires solving two elasticity systems. For both systems, the same matrix, thus the same pre-conditioner, is considered. Only boundary conditions change. If inner vertex displacement is sought for time frame [t, t + Δt], boundary conditions are imposed by the location of the body at time t + Δt∕2 and t + Δt. These locations are computed using body velocity and acceleration. Note that solving the second linear system is cheaper than solving the first one as a good prediction of the expected solution can be obtained from the solution of the first linear system. Now, to define the trajectory of each vertex, the velocity and acceleration are deduced from evaluated middle and final positions: In this context, it is mandatory to make sure that the mesh remains valid for the whole time frame [t, t + Δt], which is done by computing the sign of the volume of the elements all along their path [1].

Local mesh optimization
In order to preserve the mesh quality between two mesh deformation computations, it has been proposed [1] to couple mesh deformation with local mesh optimization using smoothing and generalized swapping to efficiently achieve large displacement in moving mesh applications. Connectivity changes are really effective in handling shear and removing highly skewed elements. Here, we briefly recall the mesh optimization procedure.
For 3D meshes, the quality of an element is measured in terms of the element shape by the quality function: where ( ) and |K| are edge length and element volume. Q(K) = 1 corresponds to a perfectly regular element and Q(K) < 2 corresponds to excellent quality elements, while a high value of Q(K) indicates a nearly degenerated element. (2) The first mesh optimization tool is vertex smoothing which consists in relocating each vertex inside its ball of elements, i.e., the set of elements having P i as their vertex. For each tetrahedron K j of the ball of P i , a new optimal position P opt j for P i can be proposed to form a regular tetrahedron: where F j is the face of K j opposite vertex P i , G j is the center of gravity of F j , j is the inward normal to F j and ( j ) the length of j . The final optimal position P opt i is computed as a weighted average of all these optimal positions {P opt j } K j ⊃P i , the weight coefficients being the quality of K j . This way, an element of the ball is all the more dominant if its quality in the original mesh is bad. Finally, the new position is analyzed: if it improves the worst quality of the ball, the vertex is directly moved to its new position.
The second mesh optimization tool to improve mesh quality is generalized swapping/local-reconnection (Fig. 1). Let and be the two tetrahedra vertices opposite the common face P 1 P 2 P 3 . Face swapping consists of suppressing this face and creating the edge = . In this case, the two original tetrahedra are deleted and three new tetrahedra are created. This swap is called 2 → 3 . The reverse operator can also be defined by deleting three tetrahedra sharing such a common edge and creating two new tetrahedra sharing face P 1 P 2 P 3 . This swap is called 3 → 2.
A generalization of this operation exists and acts on shells of tetrahedra [1,26]. For an internal edge = , the shell of is the set of tetrahedra having as common edge. The different edge swaps are generally denoted n → m, where n is the size of the shell and m is the number of new tetrahedra. In this work, edge swaps 3 → 2 , 4 → 4 , 5 → 6 , 6 → 8 and 7 → 10 have been implemented. In our algorithm, swaps are only performed if they improve the quality of the mesh. These operations are well-known in the field of mesh generation [26,56], but are not necessarily efficient in the context of this work. Notably, performing too many of them results in slow codes, whereas the use of bad quality functions results in poor quality meshes. The interest of the method used in this paper lies how and when optimizations are performed. The mesh optimizations are performed element by element, and only when they are needed. Smoothing is performed for every vertex, provided it increases the quality of the corresponding ball. Swaps are only performed, i.e., when the quality of an element decreases and passes a certain threshold. Tetrahedra are treated in quality order from the worst to the best one. The operation is performed only if it verifies quality criteria on the current position of the mesh and on the final position given by the mesh deformation. A key to performing efficient swaps in the moving mesh context is to allow a slight quality degradation in the future. Details on this optimization step can be found in [1].

Handling of boundaries
The mesh of the boundaries is moved rigidly, and the vertices are not usually moved on the surface (no displacement in the tangential directions). However, in some cases, such as when a body is moving very close to the bounding box of the domain, it can be useful to move the vertices of the bounding box as well. In this case, we can allow tangential displacement on the boundary. The risk of deforming a curved surface being too great, we only do this for planes aligned with the Cartesian frame. To do so, the displacements along the tangential axes are simply considered as new degrees of freedom. For instance, for a plane (x, y), the displacements along the x-axis and the y-axis are considered as degrees of freedom and are added to the elasticity system. The displacement along the z-axis is still set to 0, and thus is not added to the system.

Moving mesh algorithm
The overall connectivity-change moving mesh algorithm is described in Algorithm 1, where the different phases described above are put together. When coupled with a flow solver (see Sect. 3), the flow solver is called after the optimization phase. In this algorithm,  stands for meshes,  for solutions, Q for quality (see Relation (2)), ⏐ Ω h for the displacement on the boundary, and and for speed and acceleration. Δt and t are time steps whose meaning is detailed below.

CFL CFL
In Algorithm 1, three time steps appear: a large one Δt for the mesh deformation computation, a smaller one t opt corresponding to the steps where the mesh is optimized, and the solver time step t solver . Δt , is currently set manually at the beginning of the computation. After each mesh deformation solution, the quality of the mesh in the future is analyzed: if the quality is too low, the mesh deformation is problem is solved again with a smaller Δt (Algorithm 1 step 2(d)). Moreover, if the mesh quality degrades, a new mesh deformation solution is computed (Algorithm 1 step 3(g)). t opt is computed automatically, using the CFL geom parameter as described below. Determining t solver will be discussed in Sect. 3. If the solver time step t solver is greater than the optimization time step, then the solver time step is truncated to follow the optimizations. If t solver is smaller than the optimization time step-which is almost always the case-several iterations of the flow solver are performed between two optimization steps.

Moving mesh time steps
A good restriction to be imposed on the mesh movement to limit the apparition of flat or inverted elements is that vertices cannot cross too many elements on a single move between two mesh optimizations. Therefore, a geometric parameter CFL geom is introduced to control the number of stages used to perform the mesh displacement between t and t + Δt . If CFL geom is greater than one, the mesh is authorized to cross more than one element in a single move. In practice, CFL geom is usually set between 1 and 8. The moving geometric time step is given by: where h( i ) is the smallest height of all the elements in the ball of vertex P i . In practice, when coupled with a flow solver, the actual time step is the minimum between the flow solver time step and the geometric one.

Arbitrary Lagrangian Eulerian flow solver
An Arbitrary Lagrangian Eulerian (ALE) flow solver has been coupled to the moving mesh process described in Algorithm 1. In this section, we discuss in detail the implemented solver, and all the choices that were made from the numerous possibilities available in the literature.

Euler equations in the ALE framework
We consider the compressible Euler equations for a Newtonian fluid in their ALE formulation. The ALE formulation allows the equations to take arbitrary motion of the mesh into account. Assuming that the gas is perfect, inviscid and that there is no thermal diffusion, the ALE formulation of the equations is written, for any arbitrary closed volume C(t) of boundary C(t) moved with mesh velocity w: where and we have noted the density of the fluid, p the pressure, = (u x , u y , u z ) its Eulerian velocity, = ⋅ , q = ‖ ‖ , the internal energy per unit mass, e = 1∕2 q 2 + the total energy per unit mass, h = e + p∕ the enthalpy per unit mass of the flow, ext the resultant of the volumic external forces applied on the particle and the outward normal to interface

Spatial discretization
As regards spatial discretization of the solver, we use an edge-based finite-volume approach, with an HLLC Riemann approximate solver and second-order MUSCL gradient reconstruction. The main difference when translating these schemes from the standard formulation to the ALE formulation is the addition of the mesh velocities in the wave speeds of the Riemann problem.

Edge-based finite volume solver
The domain Ω is discretized by a tetrahedral unstructured mesh  . The vertex-centered finite volume formulation consists in associating a control volume denoted C i (t) with each vertex P i of the mesh and at each time t. The dual finite volume cell mesh is built by the rule of medians. The common boundary C ij (t) = C i (t) ∩ C j (t) between two neighboring cells C i (t) and C j (t) is decomposed into several triangular interface facets. The normal flux ij (t) along each cell interface is taken to be constant (not in time but in space), just like the solution ij on the interface. Rewriting System (4) for C(t) = C i (t) , we get the following semi-discretization at P i : is the mean value of state in cell C i at time t -V i is the set of all neighboring vertices of P i , i.e., the mesh vertices connected connected to P i by an edge The computation of the convective fluxes is performed monodimensionally in the direction normal to each finite volume cell interface. Consequently, the numerical evaluation of the flux function ij at interface C ij can be achieved by solving, at each time step, a one-dimensional Riemann problem in direction ij = with initial values L = i on the left of the interface and R = j on the right. The normal speed to the interface is temporarily noted for clarity.

HLLC numerical flux
The methodology provided by Batten [10] can be extended to the Euler equations in their ALE formulation. The HLLC flux is then described by three waves phase velocities: and two approximate states: where = ⋅ is the interface normal velocity and . are Roe average variables [52]. The HLLC flux through the interface is finally given by: The HLLC approximate Riemann solver has the following properties. It automatically: (1) satisfies the entropy inequality; (2) resolves isolated contacts exactly; (3) resolves isolated shocks exactly and (4) preserves positivity.

High-order scheme
The previous formulation reaches at best a first-order spatial accuracy. A MUSCL type reconstruction method has been designed to increase the order of accuracy of the scheme. The idea is to use extrapolated values ij and ji of at interface C ij to evaluate the flux, where = ( , , p) is the vector of physical variables. The following approximation is performed: ij = ( ij , ji , ij , ij ) with ij and ji linearly interpolated state values on each side of the interface: In contrast to the original MUSCL approach, the approximate "slopes" (∇ ) ij and (∇ ) ji are defined for each edge using a combination of centered, upwind and nodal gradients.
The centered gradient related to edge i j , is defined implicitly along edge i j via relation: Upwind and downwind gradients, which are also related to edge i j , are computed using the upstream and downstream tetrahedra associated with this edge. These tetrahedra are, respectively, denoted K ij and K ji . K ij (resp. K ji ) is the unique tetrahedron of the ball of P i (resp. P j ) whose opposite face is crossed by the straight line prolongating edge i j ; see Downstream K ij and upstream K ji tetrahedra associated with edge i j where is the ℙ 1 -Galerkin gradient on element K and P is the basis function associated with P. Parametrized nodal gradients are built by introducing the -scheme: where ∈ [0, 1] is a parameter controlling the amount of upwinding. For instance, the scheme is centered for = 0 and fully upwind for = 1.
The most accurate -scheme is obtained for = 1∕3 , also called the V4-scheme. This scheme is third-order for the two-dimensional linear advection problem on structured triangular meshes. In our case, for the non-linear Euler equations on unstructured meshes, a second-order scheme with a fourth-order numerical dissipation is obtained [18]. The high-order gradients are given by:

Limiter
The previous MUSCL schemes are neither monotone nor positive. Therefore, limiting functions must be coupled to the previous high-order gradient evaluations to preserve the monotonicity and positivity of the scheme. To this end, the gradient of Relation (6) is replaced by a limited gradient denoted (∇ lim ) ij . Here, the three-entry limiter introduced by Dervieux [17], which is a generalization of the SuperBee limiter, will be used : with The operator L defined above is applied component by component.

Boundary conditions
Boundary conditions are computed vertexwise. Several conditions are used, but only one, the slipping condition, is applied to moving bodies. In the context of ALE simulation, this condition has to take into account the displacement of the body. Consequently, we impose weakly 1 where i is the DGCL compatible unitary boundary face normal and i is the boundary face velocity. The standard ALE slipping boundary flux of vertex P i reduces to: �K j � is the mean outward normal of the boundary interface and . However, when high-order numerical schemes are considered, such a boundary condition creates oscillations in the density and the pressure when shock waves impact normally the boundary. We thus prefer considering a mirror state and apply an approximate Riemann state to diminish these oscillations.
We thus have to evaluate the flux between the state on the boundary and the ALE mirror state : as the mirror state verifies To evaluate the boundary flux, we consider the HLLC numerical flux between the state and the mirror state: Note that by definition we have Thus, if Condition (9) is satisfied, then i = i and the flux slip i , i , i simplifies to the form in Relation (10). In general, this condition is not satisfied, so we use Relation (11).

Time discretization
Temporal discretization is a more complex matter. In this work, we chose an explicit time discretization, which is simpler than implicit discretizations. Our time discretization is compliant with the discrete geometric conservation law, which can be used to rigorously determine when the geometric parameters that appear in the fluxes should be computed.

The geometric conservation law
We need to make sure that the movement of the mesh is not responsible for any artificial alteration of the physical phenomena involved, or at least, to make our best from a numerical point of view for the mesh movement to introduce an error of the same order as the one introduced by the numerical scheme. If System (4) is written for a constant state, assuming ext = 0 , we get, for any arbitrary closed volume C = C(t): As the constant state is a solution of the Euler equations, if boundaries transmit the flux towards the outside as it comes, we find a purely geometrical relation inherent to the continuous problem. For any arbitrary closed volume C = C(t) of boundary C(t) , Relation (12) is integrated into: which is usually known as the geometric conservation law (GCL). From a geometric point of view, this relation states that the algebraic variation of the volume of C between two instants equals the algebraic volume swept by its boundary. The role of the GCL in ALE simulations has been analyzed in [22]. It has been shown that the GCL is neither a necessary nor a sufficient condition to preserve time accuracy; however, violating it can lead to numerical oscillation [46]. In [24] the authors show that compliance with the GCL guarantees an accuracy of at least the first order in some conditions. Therefore, most would agree that the GCL should be enforced at the discrete level for a large majority of cases.

Discrete GCL enforcement
A new approach to enforcing the discrete GCL was proposed in [46,57,58], in which the authors proposed a framework to build ALE high-order temporal schemes that reach approximately the design order of accuracy. The originality of this approach consists in precisely defining which ALE parameters are true degrees of freedom and which are not. In contrast to other approaches [35,38,48], they consider that the times and configurations at which the fluxes are evaluated do not constitute a new degree of freedom to be set thanks to the ALE scheme. To maintain the design accuracy of the fixed-mesh temporal integration, the moment at which the geometric parameters, such as the cells' interfaces' normals or the upwind/downwind tetrahedra must be computed, is entirely determined by the intermediate configurations involved in the chosen temporal scheme. The only degree of freedom to be set by enforcing the GCL at the discrete level is . Incidentally, it is implicitly stated that w is never involved alone but only hidden in the term ‖ ‖ which represents the instantaneous algebraic volume swept.
In practical terms, the interfaces normal speeds are found by simply rewriting the scheme for a constant discrete solution, which leads to a small linear system that is easily invertible by hand. This procedure is detailed in the next section for one Runge-Kutta scheme. Any fixed-mesh explicit RK scheme can be extended to the case of moving meshes thanks to this methodology, and the resulting RK scheme is naturally DGCL. Even if this has not been proven theoretically, the expected temporal order of convergence has also been observed numerically for several schemes designed using this method [57].

RK schemes
Runge-Kutta (RK) methods are famous multi-stage methods to integrate ODEs. In the numerical solution of hyperbolic PDEs, notably the Euler equations, the favorite schemes among the huge family of Runge-Kutta schemes are those satisfying the strong stability preserving (SSP) property [53,55]. In what follows, we denote by SSPRK(S,P) the S-stage RK scheme of order P. We adopt the following notations: with Superscript notation X s indicates that the quantity considered is the X obtained at stage s of the Runge-Kutta process. For instance, C s i is the cell associated with vertex P i when the mesh has been moved to its sth Runge-Kutta configuration. Coefficients c s 0≤s≤S indicate the relative position in time of the current Runge-Kutta configuration: t s = t n + c s with = t n+1 − t n . Finally, we denote by A s ij the volume swept by the interface around edge e ij of cell C i between the initial Runge-Kutta configuration and the sth configuration.

Application to the SSPRK(4,3) scheme
This approach was used, for example, to build the thirdorder 4-step Runge-Kutta scheme [51], whose Butcher and Shu-Osher representations are given in Table 1. where A s ij is the volume swept by the interface around edge e ij of cell C i between the initial and the s th Runge-Kutta configuration. A natural necessary condition for the above relations to be satisfied is to have, for each interface around edge e ij of each finite volume cell C i : Therefore, the normal speed of the interface around edge e ij of cell C i must be updated in the Runge-Kutta process as follows : and the s ij and A s ij are computed on the mesh once it has been moved to the sth Runge-Kutta configuration.

Practical computation of the volumes swept
The interface of a finite volume cell is made up of several triangles, connecting the middle of an edge to the center of gravity of a face and the center of gravity of that tetrahedron; see Fig. 3. The two triangles of the interface sharing one edge within a tetrahedron are coplanar (i.e., the middle of an edge, the center of gravity of the two faces neighboring this edge and the center of gravity of tetrahedron are coplanar). The union of these two triangles is called the facet associated to the edge and the tetrahedron.
At configuration t s = t 0 + c s , (with t 0 = t n ), the outward non-normalized normal s ij and the volume swept A s ij are com-  Butcher representation Shu-Osher representation t s = t n + c s puted as described in [49]. As the cell interface is made up of several facets, the total swept volume is the sum of the volumes swept by each facet. Let us assume that the facet considered is associated with edge ij = P i P j and belongs to tetrahedron K = (P 0 , P 1 , P 2 , P 3 ) . In what follows, i ≠ j ≠ l ≠ m ∈ [[0, 3]] , G denotes the center of gravity of tetrahedron K, M m denotes the gravity center of face F m = (P i , P j , P l ) of tetrahedron K and M l the center of gravity of face F l = (P i , P j , P m ) . The outward non-normalized normal of the facet is given by: where The volume swept by the facet is: with the mean velocity written as: with Finally, the total volume swept by the interface around edge e ij of cell C i is obtained by summing over the shell of the edge, i.e., all the tetrahedra sharing the edge: It is important to understand that normals ̃ s ij,K are pseudonormals which are used only to compute the volumes swept by the facets. They must not be mistaken for normals to facets taken at t s , s ij,K , which are used for the computation of ALE fluxes.

Volumes swept by boundary interfaces
The pseudo-normals and swept volumes of boundary faces are computed in a similar way. Let K = (P 0 , P 1 , P 2 ) be a boundary triangle, as in Fig. 3. Let M 0 , M 1 and M 2 be the middles of the edges and G the center of gravity of K. The triangle is made up of three quadrangular finite volume interfaces: (P 0 , M 2 , G , M 1 ) associated with cell C 0 , (P 1 , M 0 , G , M 2 ) with C 1 and (P 2 , M 1 , G , M 0 ) with C 2 . Each boundary interface I i is made of two sub-triangles, noted T ij and T ik with j, k ≠ i.
The volume swept by interface I i between the initial and current Runge-Kutta configurations is the sum of the volumes swept by its two sub-triangles: where w s G ij is the velocity of the center of gravity G ij of triangle T ij and ̃ s ij is the pseudo-normal associated with T ij , computed between the initial and current Runge-Kutta configurations.
The six triangles T ij are coplanar, so their pseudo-normals have the same direction. Moreover, as median cells are used, their pseudo-normals also have the same norm, which is equal to one sixth of the norm of the pseudo-normal to triangle K. This common pseudo-normal is therefore equal to: Finally, the total volume swept by the boundary interface is:

MUSCL approach and RK schemes
Regarding spatial accuracy, we have seen that the order of accuracy can be enhanced using the MUSCL-type reconstruction with upwinding. However, in the ALE context, one must determine how and when upwind/downwind elements should be evaluated to compute upwind/downwind gradients which are necessary for the -schemes. This question is neither answered in the literature and generally approximations are carried out. For instance, some papers propose to use upwind and downwind elements at t n for the whole Runge-Kutta process. However, this choice should be consistent with the considered time integration scheme. Following the framework of [57], it is clear that preserving the expected order of accuracy in time imposes that the upwind/ downwind elements and the gradients are computed on the current Runge-Kutta configuration, i.e., on the mesh at t s . Therefore, similarly to geometric parameters, the upwind/ downwind elements and the gradients should be re-evaluated at each step of the Runge-Kutta stage.

Computation of the time step
The maximal allowable time step for the numerical scheme is: where h(P i ) is the smallest height in the ball of vertex P i , c i is the speed of sound, i is the Eulerian speed (the speed of the fluid, computed by the solver) and w i is the speed of the mesh vertex. The global time step is then given by = CFL min P i ( (P i )).

Handling the swaps
An ALE formulation of the swap operator, satisfying the DGCL, was proposed in 2D in [51], but its extension to 3D is very delicate because it requires to handle 4D geometry, and it has not been carried out yet. Instead of considering an ALE scheme for the swap operator, our choice in this work is to perform the swaps between two solver iterations, i.e., at a fixed time t n . This consequently means that during the swap phase, the mesh vertices do not move, and thus the swaps do not impact the ALE parameters and w , unlike [51] where the swaps are performed during the solver iteration, i.e., between t n and t n+1 . After each swap, the solution should be updated on the new configuration. Two interpolation methods are considered here. The first, and simplest, one is to perform a linear interpolation to recover the solution. As only the connectivity changes and not the vertices positions, the solution at the vertices does not change, i.e., nothing has to be done. This interpolation is DGCL compliant, since the constant state is preserved (in fact, any linear state is preserved), but it does not conserve the mass ( i.e., it does not conserve the integral of the conservative variable) which is problematic for conservative equations when discontinuities are involved in the flow.
The second method is the P 1 -exact conservative interpolation following [2,5]. It is a simplified version of the latter because the cavity of the swap configuration is fixed. The mass conservation property of the interpolation operator is achieved by element-element intersections. The idea is to find, for each element of the new configuration, its geometric intersection with all the elements of the previous configuration it overlaps and to mesh this geometric intersection with simplices. We are then able to use a Gauss quadrature formula to exactly compute the mass which has been locally transferred. High-order accuracy is obtained through the reconstruction of the gradient of the solution from the discrete data and the use of some Taylor formulae. Unfortunately, this high-order interpolation can lead to a loss of monotonicity. The maximum principle is recovered by correcting the interpolated solution in a conservative manner, using a limiter strategy very similar to the one used for finite volume solvers. Finally, the solution values at vertices are reconstructed from this piecewise linear by element discontinuous representation of the solution. The algorithm is summarized in Algorithm 2, where m K stands for the integral of any conservative quantities (density, momentum and energy) on the considered element. This method is also compliant with the DGCL.
Moreover, after each swap, the data of the finite volume cells (volume and interface normals) are updated, together with the topology of the mesh (edges and tetrahedra). This requires to have a flow solver with dynamic data. In Sects. 6.1 and 6.2, we will quantify the error due to the use of swaps in numerical simulations.

FSI coupling
The moving boundaries can have an imposed motion, or be driven by fluid-structure interaction. A simple solid mechanics solver is coupled to the flow solver described previously. The chosen approach is the 6-DOF (6 Degrees of Freedom) approach for rigid bodies.

Movement of the geometries
In this work, the bodies are assumed to be rigid, of constant mass and homogeneous, i.e., their mass is uniformly distributed in their volume. The bodies we consider will never break into different parts. Each rigid body B is fully described by: Physical quantities Its boundary B and its associated inward normal , its mass m assumed to be constant, its d × d matrix of inertia  G computed at G which is symmetric and depends only on the shape and physical nature of the solid object. 2 Kinematic quantities The position of its center of gravity = (x(t), y(t), z(t)) , its angular displacement vector = (t) and its angular speed vector d ∕dt.
ext denotes the resultant vector of the external forces applied on B, G ext the kinetic moment of the external forces applied on B computed at G and the gravity vector. We assume that the bodies are only submitted to forces of gravity and fluid pressure. The equations for solid dynamics in an inertial frame then read:

Discretization
The equation governing the position of the center of gravity of the body is easy to solve since it is linear. Its discretization is straightforward. However, the discretization of the second equation, which controls the orientation of the body, is more delicate. Since the matrix of inertia  G depends on , it is a non-linear second-order ODE. The chosen discretization is extensively detailed in [50] and is based on rewriting of the equations in the frame of the moving body.

Explicit coupling
As the geometry must be moved in accordance with the fluid computation, the same time integration scheme has been taken to integrate the fluid and the solid equations. Therefore, time-advancing of the rigid bodies ODE System is performed using the same RKSSP scheme as the one used to advance the fluid numerical solution. The coupling is loose and explicit as the external forces and moments acting on rigid objects are computed on the current configuration.

Implementation and parallelization
Most parts of the code are parallelized with p-threads using an automatic p-thread parallelization library [45] coupled with a space filling curve renumbering strategy [4]. The Hilbert space filling curve based renumbering strategy is used to map mesh geometric entities, such as vertices, edges, triangles and tetrahedra, into a one-dimensional interval. In numerical applications, it can be viewed as a mapping from the computational domain onto the memory of a computer. The local property of the Hilbert space filling curve implies that entities which are neighbors on the memory 1D interval are also neighbors in the computational domain. Therefore, such a renumbering strategy has a significant impact on the efficiency of a code. We can state the following benefits: it reduces the number of cache misses during indirect addressing loops, and it reduces cache-line overwrites during indirect addressing loops, which is fundamental for shared memory parallelism. The automatic parallelization library cuts the data into small chunks that are compact in terms of memory (because of the renumbering), then uses a dynamic scheduler to allocate the chunks to the threads to limit concurrent memory accesses. In the case of a loop performing the same operation on each entry of a table, the loop is split into many sub-loops. Each sub-loop will perform the same operation (symmetric parallelism) on equallysized portions of the main table and will be concurrently executed. It is the scheduler's job to make sure that two threads do not write on the same memory location simultaneously. When indirect memory accesses occur in loops, memory write conflicts can still arise. To deal with this issue, an asynchronous parallelization is considered rather than of a classic gather/scatter technique, i.e., each thread writes in its own working array and then the data are merged.
Our moving mesh strategy performs frequent mesh optimizations that impact the ALE speeds of the vertices, so it is much more efficient to have the whole moving mesh part included in the solver code. To handle the moving mesh, semi-dynamic data structures are necessary. Since we do not do vertex insertion or deletion, only some of the data structures have to be dynamic (edge and tetrahedra lists), whereas some remain static (the list of vertices, as well as all the data associated to vertices). In the case of dynamic data, the appropriate organization of the memory which reduces concurrent accesses can be lost, and thus imposes re-sorting the data according to the Hilbert numbering from time to time in order to maintain a good speedup.
Some parts of the code have not yet been parallelized, in particular the elasticity solution and the mesh optimization step. The optimizations are difficult to parallelize because they are greatly dependent on the order in which they are performed. Although changing this order leads to the same overall result, it is not exactly the same, and thus, the dynamic handling of the library makes each run non reproducible. To have a parallel and reproducible mesh optimization step, the operators should be rewritten with different algorithms to make them independent of the process order. Such algorithms are generally less efficient in terms of CPU time.

Verification of the solver
In the following section, we present academic test cases in order to assess the behavior of the different parts of the solver: the fluid ALE solver, the rigid body dynamic solver and the FSI coupling.

Static vortex in a rotating mesh
The first test case is used to validate the ALE solver, and study the impact of the moving mesh on the accuracy of the solution. We study the conservation of a static vortex (similar to [34]) on a rotating mesh.
This case is an extension of a 2D case: we consider an extruded cylinder, and the initial solution is constant along the z axis. The domain is a cylinder of radius 5 and of height 1. The initial solution is defined as follows. Let r 0 = 1 be a characteristic radius and the reference state be: Let us define the following quantities: For a point of cylindrical coordinates (r, , z) , the enthalpy and speed of the vortex are: And finally the initial solution is: This initial configuration results in a cylindrical vortex rotating around the z axis, that remains in a constant state over time (density, pressure and speed are preserved). To avoid errors due to boundary conditions, the exact solution is enforced for every vertex farther than 90% of the cylinder radius (typically just a few layers of elements). The mesh is rotated rigidly around the z axis. Several rotation speeds were considered, all of which led to the same conclusions. Here, we only present the results for an angular speed of 0.34° per time unit, which corresponds approximately to the speed of the vortex at a radius r = 5 . We address both space and time convergence.

Space convergence
For the space convergence study, we consider a set of uniform meshes of sizes varying from 10,000 vertices to 1.8 million vertices; see Table 2. The simulation is run until (26) ∞ = 1, u x∞ = 0, u y∞ = 0, u z∞ = 0, p ∞ = 1 and = 1.4.
time t = 540 , which corresponds to half a revolution of the cylinder. Since an exact solution is known, it is easy to compute an error with respect to this exact solution on each mesh. Two cases are compared: in one case the mesh is fixed, and in the other case the mesh is rotated as described previously. Since no mesh optimizations are performed, it allows us to make sure that the extra ALE terms are well computed. The results are shown in Fig. 4.
On the left graph, we plot the error varying with time for the different meshes, fixed (red) and rotating (blue). We can see that, as expected, the error diminishes as the size of the mesh grows, but more significantly the curves of the error for the moving and static meshes are almost superimposed. This confirms the accuracy of the ALE scheme.
To study the spatial order of convergence, the error is plotted with respect to the mesh size at different times on the right-hand graph. We can see that the order of convergence reached is slightly greater than 2, which is coherent with the space scheme used.

Time convergence
This study was carried out in 2D for reasons of efficiency. Our goal was to study the impact of the Runge-Kutta schemes used and the CFL parameter. We consider a disc of radius 5, with the same initial solution as previously (Eq. 27) and rotating with the same speed. Three time discretizations are used: a standard first-order explicit scheme (SSPRK(1,1)), a fivestep second-order Runge-Kutta scheme (SSPRK(5,2)) and a four-step third-order Runge-Kutta scheme (SSPRK(4,3)) described in Sect. 3.3.4. Four CFL numbers are set for each case: CFL max , 3∕4 × CFL max , 1∕2 × CFL max and 1∕10 × CFL max . CFL max is, respectively, 1, 4 and 2 for the schemes SSPRK(1,1), SSPRK(5,2) and SSPRK (4,3). The results are gathered in Fig. 5. On the first graph, the error over time is plotted for the different temporal schemes and CFL numbers, while on the second, we plot the error varying with the CFL number at different time steps. On both graphs, we can see that the error decreases with the order of the Runge-Kutta scheme and with the CFL number. Whereas the standard explicit scheme is very sensitive to the time step, the error does not vary much for the SSPRK(4,3) and SSPRK (5,2). What is interesting to notice is that we have to use a CFL number of 0.1 CFL max for a standard explicit scheme to reach the same level of accuracy as the SSPRK(4,3) scheme with a CFL number at its maximum or the the SSPRK(5,2) scheme with a CFL number of 0.75 times its maximum. To reach a non-dimensional time equal to 3 with the standard explicit scheme with a CFL of 0.1, 30 time steps are required, whereas only 5 are required for the SSPRK(5,2) scheme with a CFL of 3, and 6 for the SSPRK(4,3) with a CFL of 2. Thus we can go 6 times faster with the high-order Runge-Kutta schemes and still obtain the same solution accuracy, which is obviously interesting.
Influence of swaps It is important to study the impact of edge/face swaps on the solution accuracy in 3D. However, it is difficult to set up a meaningful test case with swaps. The number of swaps should be constant (in proportion to the size of the mesh),  ) 48  211  691  942  1330  1965  3122  5287  10,218   # time steps  4673  7280  15,608 14,321 14,321 17,275 29,996 32,245 44,680   0  100  200  300  400  500  600  and they should also be uniformly distributed in both time and space.
On all the simulations run with the connectivity-change moving mesh strategy, we noticed that, on average, less than one swap per 10,000 tetrahedra and per time step is performed. To fit this observation, in this example, we swap all the bad elements to preserve the mesh quality every 15 solver time steps. Then, if not enough swaps were performed, we randomly swap elements to reach a number of 1 swap per 666 tetrahedra. Not all random swaps are actually performed, because some of them affect the quality of the mesh too drastically, so the number of swaps is not perfectly controlled but this method allows us to have an average number of swaps close to one per 10,000 tetrahedra and per time step for all the meshes, evenly distributed in the domain, as stated in Table 3.
We again consider 3D meshes, rotating with the same speed of 0.34° per time unit, and the simulations are run up to 270s. In the absence of discontinuity in the solution, a linear interpolation is performed after the swaps. The results are gathered in Fig. 6. The error is slightly higher with swaps (in green) than without (in blue and red), however the error curves remain very close to those without swaps. The discrepancy can mainly be explained by the non conservative characteristic of the linear interpolation. Even without a conservative formulation, the error introduced in this example is not so large, and the same second-order convergence rate is observed. In the next example, we analyze the impact of a conservative formulation in the case of discontinuities in the solution.
To sum up, this static vortex test case allowed us to validate our ALE solver. We asymptotically reach an order of convergence of 2 for the spatial error, and we made sure there was no additional error introduced by the ALE terms. As regards temporal convergence, the importance of a highorder Runge-Kutta scheme was established. Finally, we showed that edge/face swapping-artificially created in this case, but mandatory to preserve the quality of the mesh (and thus the accuracy of the solution) when complex geometric displacements are involved-only creates a small error, even if used without a specific conservative treatment.

Sod's shock tube with a rotating circle
To show the impact of the choice of the interpolation when swaps are performed in the presence of discontinuities, we consider the well-known Sod's shock tube problem [54].  The solution is computed until non-dimensioned time t = 0.25 , and a shock wave and a contact discontinuity propagate on the right-hand side of the tube while a rarefaction wave propagates on the left-hand side. To analyze the ALE scheme and the impact of the swaps, we define a circular region of center (0.75, 0.1) and radius 0.05 which performs two full rotations within the simulation time frame; see Fig. 7. As the circular region is rotating, the mesh is sheared and swaps are performed around the circular region to maintain a good mesh quality. As the shock wave and the contact discontinuity move through the rotating region during the simulation time frame, we can analyze the impact of the interpolation stage when swaps are performed. For comparison, we also run the simulation on the same fixed mesh. Figure 8 shows the results on the fixed mesh (top), and on the moving mesh using the linear interpolation when swaps are performed (middle) and using the P 1 -exact conservative interpolation (bottom). Extraction of the density profile along the line y = 0.0483256 is given in Fig. 9. We can observe that some defects in the solution appear in the contact discontinuity when the linear interpolation is considered due to the rotation of a part of the domain. These defects disappear when the conservative interpolation is used. The defects are pointed out on the density profile when we zoom on the contact discontinuity; see Fig. 9 (right). But, they remain minor when we observe the overall profile.

Piston
The third validation test case is an FSI piston case from [37]. It is originally a 1D problem, which is extended to 3D. As shown in Fig. 10, we consider a gas contained in a 3D cylindrical chamber, closed at one end by a wall, and at the other end by a moving piston of mass m p and of cross-section A p . Besides the forces of pressure, the piston is submitted to a restoring force modeled by a spring of rigidity k p . A natural frequency of the piston can be defined by f p = 1 2 √ k p m p . Let L 0 be the length of the chamber at rest, and u(t) be the displacement of the piston with respect to the position at rest. The gas is initially at rest, and the piston is held in position u 0 . We consider that all the variables are uniform on each cross-section of the cylinder (1D assumption).
It is shown in [37] that the piston is submitted to a resultant force applied to its center of gravity:     8 Sod's shock tube problem with a rotating region. From top to bottom, density iso-values and iso-lines at non-dimensioned time t = 0.25 for the fixed mesh, the moving mesh with linear interpolation, and the moving mesh with the P 1 -exact conservative interpolation The quantities of interest are the displacement and the pressure of the piston. For the pressure, we consider the pressure at its center of gravity. These two quantities are plotted in Fig. 11. Regarding the piston movement, it is periodic, as expected, with a period T p = 1∕f p (see Table 4), the amplitude of the oscillations being damped due to FSI effects. These curves are similar to the results of [37] for the cases m p = 10 kg, 100 kg and 1000 kg. To highlight the FSI effect, we added a case with a smaller mass m p = 4 kg, where the reduction of the amplitude of the oscillations is even clearer. As concerns the piston pressure, the results also correspond to [37]. We ran the case with several mesh sizes (from 200,000 to one million vertices) and several temporal schemes, and all results are consistent. The pressure fields at a time around 2/5 of the final time are shown in Fig. 12 for the four masses considered. We can see that the 1D structure of the solution is well preserved.

Some application examples
Finally, we analyze the behavior of our strategy on two more complex, industrial-like, examples of simulations in three dimensions. These examples are challenging due both to the size of the meshes considered, and to the large displacement of the geometries. Note that a linear interpolation was used after edge/face swaps as no shocks are present in the flow-field, and it does not alter the conclusions of these examples.

Two F117 aircraft crossing flight paths
The first example models two notional F117 aircraft with crossing flight paths. This problem illustrates well the efficiency of the connectivity-change moving mesh algorithm in handling large displacements of complex geometries without global remeshing. When both aircraft cross each other, the mesh deformation encounters a large shearing due to the opposite flight directions. The connectivity-change mesh deformation algorithm easily handles this complex displacement thanks to the mesh local reconnection. Therefore, the mesh quality remains very good during the whole displacement, without any remeshing step.
As concerns the fluid simulation, the aircraft are moved with an imposed motion of translation and rotation at a speed of Mach 0.4, in an initially inert uniform fluid: at t = 0 the speed of the air is null everywhere. Transmitting boundary conditions are used on the sides of the surrounding box, while slipping conditions are imposed on the two F117 bodies. After a short phase of initialization, the flow is established when the two F117s pass each other, and the density fields around the aircraft and in their wake interact. Acoustic waves are created in front of the F117s due to the instantaneous setting in motion of the aircraft. This is not realistic, but our aim was to validate our moving mesh approach rather than run a physically realistic simulation. In Fig. 13, a zoom on the geometry of the two aircraft is shown. In Fig. 14, we show both the moving mesh aspect of the simulation and the flow solution at different time steps.
The mesh is composed of 585,000 vertices and (initially) 3.5 million tetrahedra. The whole simulation requires 22 elasticity solutions and 1600 optimization steps for a total of 2,500,000 swaps. The final mesh average quality of 1.4 is excellent and we notice that 99.8% of the elements have a quality smaller than 2 and only 52 elements have a quality higher than 5.
This simulation was run in a reasonable time: 18 h were necessary to do 39,000 time steps on a machine with 20 hyperthreaded i7 cores at 2.5 GHz. Very few elasticity solutions are requires compared to the number of solver time steps, and the total time required for the solution of the elasticity problems is only 25 min, which represents only 2.3% of the total time. The good quality of the mesh ensures an acceptable solution accuracy throughout the simulation. The optimization steps (swapping and smoothing) only take 25 min. As for the impact of the swaps, it is difficult to evaluate, since the simulation cannot be run without them. One can notice that on average, only 66 swaps per solver time step were performed,

Ejected cabin door
The last example is the FSI simulation of the ejection of the door of an over-pressurized aircraft cabin. This test case has been proposed by aircraft designers, and the aim is to evaluate when the door hinge will yield under cabin pressure. Generally, such experiments are done in a hangar and numerical simulations must be able to predict if the door will impact the observation area. From a purely moving mesh point of view, the difficulty is that the geometry is anisotropic and rolls over the elements while progressing inside a uniform mesh composed of 965,000 vertices. Another difficulty lies at the beginning of the movement when the door is ejected from its frame. The gap between the door and the frame is very small, and the mesh is sheared in the interstice; see Fig. 15. However, no remeshing is needed. The connectivity-change moving mesh strategy is able to get rid of these skewed elements quite rapidly to finally achieve an excellent quality throughout the door displacement.
From the FSI point of view, we present a simplified version of the case-and thus probably not very representative of the actual physics involved. However, we make sure that the solution is physically coherent with the simplified initial conditions. At time t = 0 , the volume is divided into two parts: outside and inside the cabin; see Snapshots of the solution at different time steps are shown in Fig. 16. We can see features similar to mach diamonds at the rear of the door, due to the high speed of the door. As concerns the door movement, the door is blown away from the cabin as expected. However, it acquires a slow rotation movement. For the time frame considered ( T end = 0.2s ), gravity has a negligible impact on the trajectory. The mesh is composed of 965,000 vertices and 5.6 millions of tetrahedra. Due to the difficult part where the door is moving in its frame, the whole simulation requires 29 elasticity solutions and 540 optimization steps for a total of 875,000 swaps. The final mesh average quality of 1.33 is excellent and we notice that 99.95% of the elements have a quality smaller than 2 and only 1 element has a quality higher than 5.
The total time of the simulation is 2 h on a machine with 20 hyperthreaded i7 cores at 2.5 GHz, for 1122 solver time steps, including 43 min of elasticity solution. The number of elasticity solutions is still small compared to the number of time steps (29 elasticity solutions and 1122 solver time steps), and elasticity solutions do not account  for more than a third of the total time. Note that the a priori unpredictable trajectory of FSI problems generally results in a slightly higher number of elasticity steps than with analytic imposed motion. The optimization step only takes 3 min. Once again, only an average of 0.0001% swaps per tetrahedra and solver time step are performed, which is not a lot.

Conclusion
A strategy to run complex three-dimensional simulations with moving geometries has been presented in this paper. This strategy lies first on a robust connectivity-change moving mesh algorithm, which couples an elasticitylike mesh deformation method and mesh optimizations. In particular, a reduced number of elasticity solutions improves the efficiency of the algorithm, while edge/face swapping makes it possible to deal with the strong shearing of the mesh that occurs when the geometries undergo large displacement.
To run CFD simulations on such moving meshes, a finite volumes flow solver for the compressible Euler equations has been extended to the ALE framework. As regards the temporal accuracy, the SSPRK schemes considered are based on the strict application of the discrete geometric conservation law (DGCL), which is supposed to preserve the accuracy of the schemes. The displacement of the geometries can be either imposed a priori, or governed by a 6-DOF fluid-structure coupling.
The goal of the paper was to demonstrate that this strategy allows us to run complex three-dimensional simulations. Challenging examples of such simulations, were presented and analyzed, with imposed motion or fully FSI. Despite the large displacement of the geometries, no global remeshing was necessary to perform the simulations, while preserving a good mesh quality. The unsteady solutions of these examples are physically coherent, and the total CPU time of the simulations is reasonable. The use of high-order Runge-Kutta schemes was shown to help reduce the cost of the simulation. Whereas the handling of the moving mesh usually takes up a large proportion of the simulation total CPU time, our connectivity-change moving mesh algorithm only accounts for a small fraction of that CPU time, most of the time being spent on actually solving the equations.
Some issues still need to be addressed. A thorough analysis of the impact of swaps on the solution accuracy needs to be carried out. We have shown that a conservative interpolation method improves significantly the accuracy of the solution in the presence of discontinuities. A comparison of this approach with the one described in [51] would probably be revealing. The elasticity time step is currently constant and fixed a priori for the whole simulation. The efficiency of the moving mesh part would be improved if this time step could be automatically adapted to the movement of the geometries. Finding such an optimal time step is harder than it seems. The case of non rigid bodies also has to be addressed: the moving mesh algorithm can handle deformable bodies [1,7], but the FSI coupling would be much more complex.
Finally, this entire strategy was devised to fit within our metric-based mesh adaptation framework [42]. An unsteady version of the adaptation algorithm already exists [3]. It was extended to the case of moving meshes in 3D recently [8], and a goal-oriented version of the process is developed.

Software used
The solver described in this paper and used to run all the simulations is Wolf. The meshes were generated using GHS3d [27] and Feflo.a [43]. The visualization of the meshes and solutions was done with Vizir [41]. Fig. 15 Test case of the door ejection. Surface mesh and cut plane just after the beginning of the simulation. Note that the gap between the door and its frame is very small, and meshed with only one layer of elements, which increases the difficulty of the moving mesh problem Fig. 16 Test case of the door ejection. Snapshots of the moving geometry and the mesh (left) and pressure (right)