Extension of anisotropic metric-based mesh adaptation to moving mesh simulations 

Get Complete Project Material File(s) Now! »

Connectivity-change moving mesh strategy

When dealing with a unique body-fitted moving mesh, two main approaches are used.
The first one is simply to move the mesh for as long as possible, with a constant topology, and remesh when the quality becomes critical. For each remeshing, the simulation has to be stopped, a new mesh must be generated (a whole new mesh or parts of it), new data structures must be created, 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. Between two remeshings the simulation is fully ALE and free from interpolation errors due to con-nectivity changes. However, for larger displacements, the number of remeshings increases to prevent invalid elements from appearing, especially if shearing appears. This can be costly in terms of CPU, but also result in poor accuracy due to the numerous solution transfer steps. This strategy is for instance used in [Murman 2003] with hexaedral Cartesian meshes, or [L¨ohner 1990a, Hassan 2000, Jayaraman 2000] on unstructured meshes.
The second approach tries to preserve a good mesh quality throughout the simulation by performing local remeshing operations such as vertex insertion, vertex collapse, connectivity changes and vertex displacements. It is used for example in [Dobrzynski 2008, Comp`ere 2010]. The advantage of this method is that it maintains an acceptable mesh quality without need-ing 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 tend to lead to a loss of accuracy due to interpolation errors.
In our case, we want to address the case of large and complex displacements of the bound-aries, with potentially highly non-linear trajectories as may occur in FSI simulations. This requires that the moving mesh strategy fulfills a number of properties:
• It must above all preserve the validity of the mesh, and never create elements with a negative volume.
• It must preserve the quality of the elements of the mesh all along the simulation, to ensure a good solution accuracy and reasonable solver time steps.
• It must be efficient in terms of CPU since the mesh must be moved at each solver iteration. It must not take for more than a small fraction of the total simulation time.
• It must be robust, and be able to deal with very small elements, and large shearing con-straints.
• It must be able to preserve the anisotropic adaptation if any.
For efficiency purposes, we add that we want to remesh as little as possible, and that most cases should be run without remeshing at all.
Remark 1. What we call remeshing entails: stopping the simulation, calling meshing software to generate either a whole new mesh or parts of it, and interpolate the solution on the new mesh. This is costly because (a) it usually requires the use of an external code, with its own data structures, requiring CPU-expensive data transfers, (b) the meshing step can be long (and might not even complete) in the case of complex geometries, (c) an interpolation step which produces not too much error is long too, especially if the whole mesh or large portions of it are considered. On the other hand, our mesh optimizations (a) are performed on the go, within the numerical solver code, using the same semi-dynamic data-structures, (b) concern only a limited set of elements, and are performed element by element, only if the quality of the element requires them, (c) no global interpolation procedure is required since the number of vertices is constant.
The previous strategies do not match these requirements, either in terms of preserving the mesh quality or in terms of efficiency for large displacements. That is why we adopt a strategy based on only vertex displacements and connectivity changes. Improvements in the way to compute the displacement of the vertices makes that strategy very efficient, while the connectivity changes makes it robust, as will be demonstrated in Section 1.3. Besides its CPU efficiency, this strategy presents some very interesting features considering numerical simulations:
• Since no vertex insertion or deletion occurs, and no remeshing is performed, a source of large spoiling interpolation errors is removed, and the accuracy of the numerical simulations should be improved. The use of an ALE treatment for connectivity changes would allow the simulation to be fully ALE.
• It is especially efficient with a vertex-centered Finite-Volume approach, since the number of vertices is preserved, the number of dual cells is preserved as well.
• Our unsteady adaptation algorithm, that will be detailed in Part II, Chapter 5, requires a set number of vertices within each adaptation sub-interval to control the error. Our strategy keeps constant the number of vertices and is thus compatible with our adaptation algorithm.
One might wonder why perform connectivity changes (edge/face swaps) in the first place, if they are a priori not ALE compliant. First, performing local mesh modifications within the solver is far simpler and more efficient than remeshing globally. Moreover, connectivity changes can be relatively simply interpreted in terms of evanescent cells for purposes of 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).
In the present chapter, I will first detail the moving mesh algorithm that was used during this thesis (Section 1.1), and present a variant of the algorithm developed during the thesis (Section 1.2). Then I will present various examples of purely moving mesh simulations to study the efficiency of the algorithm and compare its two variants. The behavior of the algorithm with viscous-like boundary layers is studied in Section 1.4. In these sections, we only consider displacements with no volume variation, or a small one. The case of large volume deformation is tackled at the end of the chapter in Section 1.5.
The works presented in this Chapter were published in [Barral 2014d].

Our moving mesh algorithm

 A two step process

The connectivity-change moving mesh algorithm (MMA) is a two-step process:
• A mesh deformation step, in which a trajectory is assigned to each inner vertex, depending on the displacement of the boundaries.
• A mesh optimization step, in which the trajectories of the vertices are modified, and con-nectivity changes are performed, to preserve the quality of the mesh.

Mesh deformation step

For body-fitted moving mesh simulations, the whole mesh must be deformed to follow the moving boundaries. During the mesh deformation step, a displacement field is computed for the whole computational domain, given the displacement of its boundaries. Trajectories, or in other words, positions at a future solver time step, can thus be assigned to inner vertices.
Several techniques can be found to compute this displacement field:
• implicit or direct interpolation [de Boer 2007, Luke 2012], in which the displacement of the inner vertices is a weighted average of the displacement of the boundary vertices,
• solving PDEs – the most common of which being Laplacian smoothing [L¨ohner 1998], Winslow equation [Masters 2015], a spring analogy [Degand 2002] and a linear elasticity analogy [Baker 1999].
It is this last method that we selected at first, due to its robustness in 3D [Yang 2007]. In Section 1.2, a direct interpolation procedure will also be considered.
In the case of the elasticity analogy, the computational domain is assimilated to a soft elastic material, which is deformed by the displacement of its boundaries. More precisely, the inner vertices movement is obtained by solving an elasticity-like equation with a P1 Finite Element Method (FEM):
div( ( E )) = 0 , with E = rd + T rd , (1.1) where and E are respectively the Cauchy stress and strain tensors, and d is the Lagrangian displacement of the vertices. The Cauchy stress tensor follows Hooke’s law for isotropic homo-geneous medium, where ⌫ is the Poisson ratio, E the Young modulus of the material and , µ are the Lam´e coefficients:
(E) = trace(E) Id + 2 µ E or 1 + ⌫ ⌫
E( )= trace( ) Id .
E E
In our context, ⌫ is typically chosen of the order of 0.48, which corresponds to a very soft, nearly incompressible material 1. 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 then solved with a P1 Finite Element Method (FEM) by a Conju-gate Gradient algorithm coupled with an LU-SGS pre-conditioner (See [Olivier 2011a] for more details).

Basic elasticity-based algorithm

The classic algorithm to move the mesh using the elasticity-based deformation method is simply to solve the elasticity problem at every time step and then move the mesh using the solution of the problem. This strategy is described in Algorithm 1. In this algorithm, H stands for meshes, S for solutions and v for vertex speed.
Algorithm 1 Basic Moving Mesh Algorithm
Input: H0, S0
While (t < T end)
1. t = Get solver time step Hk, Sk, v, CF L
2. Solve mesh deformation: compute vertices trajectories
(a) dbody = Compute body vertex displacement from current translation speed vbody, rotation speed ✓body and acceleration abody for [t, t + t] body
(b) d = Solve elasticity system d|@⌦h , t
(c) v = Deduce inner vertices speed from displacement d
3. Sk+1 = Solve Equation of State Hk, Sk, t, v
4. Hk+1 = Move the mesh Hk, t, v
5. t = t + t
EndWhile
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. Performing it at every solver time step makes it all the more so. In order to reduce significantly the CPU time of that part, this step needs to be enhanced.

Improvements to the mesh deformation steps

Several improvements to the mesh deformation step are brought compared to the basic approach.
Local material properties alteration. First, the linear elasticity analogy allows us to modify the local material properties according to the distortion and e↵orts born by each element. In particular, small elements are made sti↵er than others, which has the advantage of limiting their deformation, and thus avoiding the creation of false elements. To do so, the way the Jacobian of the transformation from the reference element to the current element is accounted for in the FEM matrix assembly is modified as in [Stein 2003]. The classic P1 FEM formulation of the linear elasticity matrix leads to the evaluation of quantities of the form:
+ 2µ and |K| is the volume of tetrahedron K. The above quantity where > 0 is the sti↵ening power and K is the reference element. This technique locally multiplies and µ by a factor proportional to |K|χ . determines how sti↵er than large elements small elements are. In this work, we chose = 1.
Rigidify regions. Some regions of the mesh can also be rigidified, i.e. their vertices are all moved with the same displacement. Typically, this is done for a few layers of tiny elements around moving rigid bodies, that are moved together with the body. Rigidifying elements close to complex features of the geometries reduces the risk of creating bad elements in those regions, and also removes sti↵ elements from the FEM matrix. We can also chose to enclose moving bodies with a complex geometry into boxes or spheres, and rigidify the whole region inside these boxes and spheres: the objects moved are then much simpler.
Elasticity on reduced regions. The elasticity can also be solved on a reduced region: if the domain is very large compared to the size of the moving bodies and the length of the displacement, it is sometimes possible to solve the elasticity problem for a few layers around the moving bodies only. The idea is the same as the previous one, but now it is the far field that is rigidified. It has to be noted that the elastic region has to be large enough to let the moving bodies evolve in it.
Elasticity dedicated mesh. A way to reduce significantly the CPU time of the mesh defor-mation computation is to consider two meshes: one on which the physical equations are solved, and one on which the elasticity problem is solved. The elasticity-dedicated mesh can be much coarser than the computational one, since a precise elasticity solution is not needed – especially when using mesh optimizations as described below. This strategy could also be mandatory when considering adapted meshes with high aspect-ratio elements, on which the FEM elasticity system could be practically impossible to solve.
More precisely, we consider two meshes with the same boundary mesh, the elasticity-dedicated mesh being coarser than the computational one. The displacement of the boundary vertices is the same for both meshes (whether it is imposed to the bodies, or whether a FSI motion is considered: in this case, the displacement of the bodies is computed on the computational mesh, then transferred to the elasticity mesh). The elasticity FEM system is then assembled on the coarse mesh, and solved. The solution of the elasticity problem is then transferred to the computational mesh using a P2-Lagrangian scheme, which is sufficiently accurate considering the intrinsic smoothness of the solution of the above elasticity equation. Then both meshes are moved with the same time step. If mesh optimizations are performed, they are performed on both meshes, so that the elasticity mesh remains valid too. The cost of moving the coarse mesh is shown to be very small compared to the gain in solving the elasticity problem.
Reduced number of solutions. The main improvement to the basic algorithm 1 is that we do not solve an elasticity problem at each FSI solver time step t. Instead, we consider a larger time step Δt, and the elasticity problem is solved for this time step, and the trajectory of the vertices is kept constant on this whole larger time step. While there is a risk of a less good mesh displacement solution, this strategy is worthy if mesh optimizations are used to preserve the mesh quality, as explained in Section 1.1.3.
High order trajectories. Since the trajectories are computed for a large time step, it is worthwhile to consider vertex paths that are more complex than mere straight lines (linear trajectories). This is especially true if the moving bodies undergo rotation or accelerated move-ments. The paths of inner vertices is improved providing a constant acceleration a to each vertex in addition to its speed, which results in an accelerated and curved trajectory. During time frame [t, t + Δt], the position x and the velocity v of a vertex are updated as follows: x(t + t) = x(t) + t v(t) + t2
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 con-sidered. 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:
Δt v(t) = 3 x(t)
Δt2 a = 2 x(t)
+ 4 x(t + Δt/2)
4 x(t + Δt/2)
x(t + Δt)
+ 2x(t + Δt) .
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 computing by the sign of the volume of the elements all along their path [Alauzet 2014a].

Mesh optimization step

Mesh optimizations are performed regularly to preserve the quality of the mesh. The optimiza-tion procedure consists in one or several passes of vertex smoothing and one or several passes of generalized swapping, and is described below.
For 3D adapted meshes, an element’s quality is measured in terms of the element’s shape by the quality function: QM(K) = P 2 [1, +1] , (1.2)
where `M(e) and |K|M are the edge lengths and element volume in metric M. Metric M is a 3 ⇥3 symmetric positive definite tensor prescribing element sizes, anisotropy and orientation to the mesh generator (see Chapter 4). QM(K) = 1 corresponds to a perfectly regular element and QM(K) < 2 correspond to excellent quality elements, while a high value of QM(K) indicates a nearly degenerated element. For non-adapted (uniform) meshes, the identity matrix I3 is chosen as metric tensor.
Mesh smoothing. 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 Pi as their vertex. For each tetrahedron Kj of the ball of Pi, a new optimal position Pjopt for Pi can be proposed to form a regular tetrahedron: = Gj + r 2 n Pjopt j , 3 `(nj)
where Fj is the face of Kj opposite vertex Pi, Gj is the center of gravity of Fj, nj is the inward normal to Fj and `(nj) the length of nj. The final optimal position Piopt is computed as a weighted average of all these optimal positions {Pjopt}KjP i , the weight coefficients being the quality of Kj. 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.
Edge/face swapping. The second mesh optimization tool to improve mesh quality is gener-alized swapping/local-reconnection. Let ↵ and be the two tetrahedra vertices opposite the common face P1P2P3. Face swapping consists of suppressing this face and creating the edge e = ↵ . 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 P1P2P3. This swap is called 3 ! 2.
A generalization of this operation exists and acts on shells of tetrahedra [Alauzet 2014a, Frey 2008]. For an internal edge e = ↵ , the shell of e is the set of tetrahedra having e as common edge. From a shell of size n, a non-planar pseudo-polygon formed by n vertices P1…Pn can be defined. Performing a three-dimensional swap of edge ↵ requires (i) suppressing edge ↵ and all tetrahedra of the shell, (ii) defining a triangulation of the pseudo-polygon P1…Pn and (iii) finally creating new tetrahedra by joining each triangle of the pseudo-polygon with vertices ↵ and . The number of possible triangulations depends on n the number of vertices of the pseudo-polygon. The di↵erent swaps are generally denoted n ! m where 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.

Table of contents :

Introduction 
I 3D FSI moving mesh simulations
Introduction
1 Connectivity-change moving mesh strategy 
1.1 Our moving mesh algorithm
1.1.1 A two step process
1.1.2 Mesh deformation step
1.1.3 Mesh optimization step
1.1.4 Optimization procedure
1.1.5 Handling of boundaries
1.1.6 Algorithm
1.1.7 Choice of the di↵erent parameters for a robust algorithm
1.2 Mesh deformation with Inverse Distance Weighted method
1.2.1 IDW method
1.2.2 Implementation
1.3 Examples
1.3.1 Rigid-body examples
1.3.2 Deformable-body examples
1.4 Boundary layers for deformable geometries
1.4.1 One boundary layer
1.4.2 Several boundary layers
1.5 Large volume variations
1.5.1 Engine piston
1.6 Conclusion
2 ALE solver 
2.1 Euler equations in the ALE framework
2.2 Spatial discretization
2.2.1 Edge-based Finite-Volume solver
2.2.2 HLLC numerical flux
2.2.3 High-order scheme
2.2.4 Limiter
2.2.5 Boundary conditions
2.3 Time discretization
2.3.1 The Geometric Conservation Law
2.3.2 Discrete GCL enforcement
2.3.3 RK schemes
2.3.4 Application to the SSPRK(4,3) scheme
2.3.5 Practical computation of the volumes swept
2.3.6 MUSCL approach and RK schemes
2.3.7 Computation of the time step
2.3.8 Handling the swaps
2.4 FSI coupling
2.4.1 Movement of the geometries
2.4.2 Discretization
2.4.3 Explicit coupling
2.5 Implementation
2.5.1 Non-dimensionalization
2.5.2 Parallelization
2.6 Conclusion
3 Numerical Examples 
3.1 Validation of the solver
3.1.1 Static vortex in a rotating mesh
3.1.2 Flat plate in free fall
3.1.3 Piston
3.2 Numerical examples
3.2.1 AGARD test cases
3.2.2 F117 nosing up
3.2.3 Two F117 aircraft crossing flight paths
3.2.4 Ejected cabin door
3.3 Parallel performance
3.4 Conclusion
Conclusion
II Extension of anisotropic metric-based mesh adaptation to moving mesh simulations 
Introduction
4 Basics of metric based mesh adaptation 
4.1 State of the art
4.1.1 Meshing status
4.1.2 History of metric-based mesh adaptation
4.1.3 Other mesh adaptation approaches
4.2 Principle of metric-based adaptation
4.2.1 Euclidian and Riemannian metric spaces
4.2.2 Unit mesh
4.2.3 Operations on metrics
4.2.4 The non-linear adaptation loop
4.3 The continuous mesh framework
4.3.1 Duality discrete-continous: a new formalism
4.3.2 Continuous linear interpolation
4.3.3 Summary
4.4 Multiscale mesh adaptation
4.4.1 Optimal control of the interpolation error and optimal meshes
4.4.2 Control of the error in Lp norm
4.5 Conclusion
5 Unsteady mesh adaptation 
5.1 Error estimate
5.1.1 Error model
5.1.2 Spatial minimization for a fixed t
5.1.3 Temporal minimization
5.1.4 Error analysis for time sub-intervals
5.1.5 Global fixed-point mesh adaptation algorithm
5.2 From theory to practice
5.2.1 Computation of the mean Hessian-metric
5.2.2 Choice of the optimal continuous mesh
5.2.3 Matrix-free P1-exact conservative solution transfer
5.2.4 The remeshing step
5.2.5 Software used
5.3 Choice of the mean Hessian-metric
5.4 Numerical examples
5.4.1 2D shock-bubble interaction
5.4.2 3D circular blast
5.4.3 3D shock-bubble interaction
5.4.4 3D blast on the London Tower Bridge
5.5 Parallelization of the mesh adaptation loop
5.5.1 Choices of implementation
5.5.2 Analysis of parallel timings
5.5.3 Conclusion
5.6 Conclusion
6 Extension of unsteady adaptation to moving meshes 
6.1 ALE metric
6.2 Analytic examples
6.2.1 Procedure
6.2.2 Functions considered
6.2.3 Results
6.3 Update of the adaptation algorithm
vi Contents
6.3.1 Error analysis
6.3.2 Algorithm
6.3.3 Update of the metric for optimizations
6.3.4 Handling of the surface
6.4 3D numerical examples
6.4.1 Shock tube in expansion
6.4.2 Moving ball in a shock tube
6.4.3 Nosing-up f117
6.4.4 Two F117 aircraft flight paths crossing
6.5 Conclusion
Conclusion
Conclusion and perspectives 
Acknowledgments
Appendices

GET THE COMPLETE PROJECT