Mesh deformation techniques in fluid-structure interaction: robustness, accumulated distortion and computational efficiency

An important ingredient of any moving-mesh method for fluid-structure interaction (FSI) problems is the mesh deformation technique (MDT) used to adapt the computational mesh in the moving fluid domain. An ideal technique is computationally inexpensive, can handle large mesh deformations without inverting mesh elements and can sustain an FSI simulation for extensive periods ot time without irreversibly distorting the mesh. Here we compare several commonly used techniques based on the solution of elliptic partial differential equations, including harmonic extension, bi-harmonic extension and techniques based on the equations of linear elasticity. Moreover, we propose a novel technique which utilizes ideas from continuation methods to efficiently solve the equations of nonlinear elasticity and proves to be robust even when the mesh is subject to extreme deformations. In addition to that, we study how each technique performs when combined with the Jacobian-based local stiffening. We evaluate each technique on a popular two-dimensional FSI benchmark reproduced by using an isogeometric partitioned solver with strong coupling.


Introduction
Fluid-structure interaction (FSI) constitutes a class of problems involving two-way dependence between structural objects and a fluid. As such, FSI is a vast topic with applications spanning a spectrum from aerospace and civil engineering [1] to biomechanical and cardiovascular simulations [9]. In FSI problems, the fluid exerts a force on the structure which deforms in response. As the structure moves, it changes the shape of the domain occupied by the fluid together with the fluid motion and, as a result, the force that the fluid exerts on the structure. This two-way coupling between the fluid and structure behavior as well as the necessity to accommodate the fluid domain motion in both the continuous and discrete formulations of the problem is what makes FSI so notoriously complex.
Since FSI problems rarely admit analytical solutions, computational methods are widely adopted in FSI research. Here, one can distinguish between the static-mesh and moving-mesh methods. While the former attempt to resolve the motion of the fluid domain implicitly, for example by means of a stationary background Cartesian mesh [21], the latter deal with the motion of the fluid domain by tracking its boundary and adapting the computational mesh correspondingly. In this work, we study and compare various mesh deformation techniques (MDTs) which can be used to adapt the fluid mesh if the moving-mesh methods are used. The main focus here lies on the robustness of a given technique, that is, how much mesh distortion it introduces and how much mesh deformation it can handle without entangling or inverting mesh elements. Additionally, we pay attention to computational complexity of MDTs, which can significantly increase the overall FSI simulation time if left unchecked.

arXiv:2006.14051v1 [cs.CE] 19 Jun 2020
To study performance of MDTs in their natural habitat, we employ a popular two-dimensional FSI benchmark from [33]. In the benchmark, an unstable flow of a viscous fluid leads to development of the vortex shedding phenomenon which results in oscillations of a flexible beam structure. The oscillations grow in magnitude until they reach a stable periodic regime which lends itself well to studying possible long-term effects of MDTs on the fluid mesh. One of such long-term effects is the accumulated mesh distortion where mesh elements increasingly become permanently distorted, deteriorating the simulation accuracy. In addition to the original benchmark, we employ its simplified version with no fluid mechanics involved to perform a large number of computationally inexpensive tests. In this way, we can concentrate on mesh deformation and conduct a detailed analysis of MDT behavior.
To reproduce the FSI benchmark, we use a partitioned solver with strong coupling and Aitken relaxation [17]. Although modern space-time (ST) methods are becoming increasingly common in FSI [1,32,25], we resort to classical arbitrary Lagrangian-Eulerian (ALE) methods [1,20,12] which are more straightforward in implementation. Our choice in favor of a basic partitioned ALE solver is justified since both ST and ALE methods make use of the same MDTs, so we can focus on mesh deformation in this paper.
All MDTs we consider here are based on solution of elliptic partial differential equations. These include existing techniques such as harmonic extension [8,36], biharmonic extension [19] and a widely adopted technique based on linear elasticity theory [31,14]. Moreover, we propose an efficient MDT based on the equations of nonlinear elasticity and a logarithmic neo-Hookean material law which we refer to as tangential incremental nonlinear elasticity (TINE). Although techniques based on the equations of nonlinear elasticity have been proposed before [27,24], TINE is novel in that it uses the idea of a tangential continuation method [7,22] to efficiently solve the corresponding nonlinear equations. As a result, TINE is only slightly more computationally expensive than the linear-elasticity-based techniques which are linear in nature. On the other hand, it can handle as much mesh deformation but does not suffer from the accumulated distortion effect.
Robustness of any MDT can be increased by additional augmentations. Probably the most popular one is the Jacobian-based local stiffening [30,23] which turns individual mesh elements stiffer or softer depending on their size and shape. In this work, we study how different MDTs react to the Jacobian-based local stiffening. Although not considered here, further possible MDT augmentations include solid layer extension [29] and element relaxation [27,26].
The research we present in this work has been conducted in the framework of isogeometric analysis [11,6]. Despite that, the results are applicable to classical finite elements methods or any other mesh-based method for solving partial differential equations.
The rest of this paper is structured as follows. Section 2 outlines geometry and settings of the FSI benchmark and fixes the necessary notation. In Section 3, we describe various MDTs considered in this work as well as the Jacobian-based local stiffening. In Section 4, we consider the simplified benchmark and conduct a detailed analysis of the short-term and long-term behavior of different MDTs in artificial FSI-like conditions. After that, we proceed to performing a full FSI simulation of the benchmark in Section 5. We study the performance of the MDTs and check if the choice of a particular MDT affects the simulation results. Finally, we discuss the results of the MDT analysis and FSI simulations, draw conclusion and outline further research directions in Section 6.

Benchmark description
In this section, we describe the geometry and mathematical model of the FSI benchmark from [33]. It studies the flow of an incompressible Newtonian fluid through a 2D channel as the fluid interacts with a submerged structure.  Figure 1 illustrates the setting; note that the geometry is intentionally non-symmetric.
The top and bottom walls of the channel are impermeable to the fluid. The fluid enters the channel through the left wall with a prescribed velocity and exits through the right wall freely. The presence of the submerged structure changes the flow of the fluid and, in response, the fluid exerts a certain force on the structure. Thus, this is an FSI problem. Depending on the prescribed material parameters of the fluid and beam as well as on the inflow fluid velocity, this fluid force can result in a noticeable deformation of the beam. The beam deformation alters the shape of the flow channel, the flow itself and, as a consequence, the force exerted on the structure by the fluid. Such systems with a twoway interaction between components are called coupled systems.
In what follows, when we use the word "structure", we refer to the flexible beam only. If the rigid disk is also considered, we explicitly mention it. Throughout this work, we use subscripts s and f to distinguish between objects related to the structure and fluid respectively. Moreover, we use subscript a when dealing with ALE mappings and related objects. Let Ω f (t) ⊂ R 2 and Ω s (t) ⊂ R 2 denote domains occupied by the fluid and structure respectively at time t ∈ [0, T ]. Additionally, we use Ω 0 f = Ω f (0) and Ω 0 s = Ω s (0) to denote the fluid and structure domains at time t = 0. Furthermore, we use Γ (t) to denote the FSI interface ∂Ω f (t) ∩ ∂Ω s (t) where the interaction between the fluid and structure takes place. Correspondingly, Γ 0 = Γ (0) denotes the FSI interface at time t = 0. In the rest of this section, we briefly formulate the equations that we use to describe the motion of the beam and the fluid as well as their interaction. The main goal here is to fix the notation necessary for the ensuing description of mesh deformation techniques (MDTs). We encourage the reader to consult the following references for more information on: finite element discretization for the Navier-Stokes equations for incompressible flows [13]; finite element discretization for nonlinear elasticity problems [35,2]; ALE formulation of FSI problems [19,1]; partitioned solution approach to FSI problems [17].

Structure motion
We assume that the structure behavior can be characterized as elastic and compressible. Let the structure in its initial (undeformed) configuration occupy the domain Ω 0 s . We can describe the structure motion in terms of a displacement field u s : Ω 0 s × [0, T ] → R 2 . Note that displacement u s can describe a rigid body motion with no deformation involved. Information on whether actual deformation of Ω 0 s takes place is contained in the deformation gradient F s = I + ∇u s . Two important objects derived from the deformation gradient are the Green-Lagrange strain tensor E s = (F T s F s − I)/2 and the Jacobian determinant J s = det F s .
For the material behavior, we use the St. Venant-Kirchhoff constitutive law which links the second Piola-Kirchhoff stress tensor S s to the Green-Lagrange strain tensor E s : The material law (1) includes the Lamé parameters λ s and µ s , which are constitutive parameters describing physical properties of the material. They can be computed from Young's modulus E s and Poisson's ration ν s as Second Piola-Kirchhoff stress tensor S s measures forces appearing in the deformed structure with respect to its initial configuration Ω 0 s . In FSI applications, it is important to have an ability to express these forces with respect to the deformed configuration Ω s (t). This can be achieved by means of the Cauchy stress tensor σ σ σ s , which is related to S s by In the presence a given external acceleration g : , the displacement u s should satisfy the local conservation equations of linear momentum Here, ρ s is the structure density which we assume to be constant in Ω 0 s , and P s = F s S s is the first Piola-Kirchhoff stress tensor. Note that equations (4) are formulated in the stationary domain Ω 0 s .

Fluid motion
To describe the fluid motion, we use the Navier-Stokes equations for incompressible flows. If a computational domain Ω f does not change with time, the Navier-Stokes equations have the following form: Here, v f : Ω f × [0, T ] → R 2 is a vector field describing the fluid velocity at a given point in Ω, ρ f is a constant fluid density, and σ σ σ f denotes a Cauchy stress tensor.
Behavior of an incompressible Newtonian fluid is characterized by the following constitutive law: where p f : Ω f × [0, T ] → R is a pressure field, and ν f denotes the kinematic viscosity of the fluid.
In FSI applications, one has to account for the motion of the fluid domain. In this work, we consider a common strategy based on the arbitrary Lagrangian-Eulerian (ALE) mappings. An ALE mapping describes the motion of the fluid domain in terms of an auxiliary displacement field u a : With an ALE mapping, the Navier-Stokes equations can be formulated in the moving domain Ω f (t): Since the fluid domain motion is driven by the structure deformation, the ALE displacement u a has to comply with the structure displacement u s on the FSI interface: Inside the fluid domain, u a can be chosen arbitrary, hence the term ALE. The only condition is that the ALE displacement should define an invertible deformation of the fluid domain. This means that the following condition has to hold: The scope of this work is to compare different options for defining the ALE displacement in the fluid domain provided the structure displacement on the FSI interface.

Interaction conditions
Physical interaction between the fluid and structure takes place on the FSI interface Γ (t). It is characterized by the following two coupling conditions: the kinematic continuity which assures that the fluid stays attached to the structure; and the dynamic continuity which maintains the balance of forces on the FSI interface. If a partitioned approach to FSI is used, one has to enforce the coupling conditions by exchanging information between the fluid and structure solvers.

Initial and boundary conditions
To complete the definition of the benchmark, we provide suitable initial and boundary conditions. The system is initialized with zero initial conditions for the fluid velocity, structure displacement and structure velocity: The main energy source of the system is an inflow boundary condition on the fluid velocity prescribed on the left end of the channel ∂Ω in f . The condition prescribes a parabolic velocity profile with the maximum inflow velocity v max serving as an adjustable parameter. In order to comply with the initial conditions, v par is made time-dependent by scaling it with time: The resulting time-dependent inflow boundary condition v f = v in provides a smooth warm-up phase for the simulation. On the right wall of the channel ∂Ω out f , a do-nothing condition σ σ σ f ·n = 0 is prescribed. Additionally, a no-slip boundary condition v f = 0 is set on the rest of the channel wall ∂Ω ns

Geometry parametrization
The research we present in this work has been conducted in the framework of isogeometric analysis [11,6].
In the spirit of IGA, we model computational domains as collections of tensor-product NURBS patches [18].
For the structure, we use a single quadratic NURBS patch, whereas the fluid domain is modeled with seven patches. Figure 2 illustrates the isoparametric mesh in the fluid domain after several applications of uniform hrefinement. For simplicity, we use matching parametrizations for neighboring fluid patches as well as for the structure patch and the surrounding fluid patches. This choice does not restrict the applicability of our results; however, it significantly simplifies the exchange of coupling information between the fluid and structure solvers. Let us consider one of the NURBS patches forming the initial configuration of the fluid domain. Its → Ω 0 f can be written in terms of control points c k ∈ R 2 and tensor-product NURBS basis functions N k : In IGA, we the same NURBS basis functions to approximate solutions to PDEs. In our case, we use N k to discretize the ALE displacement u a in space. Assuming that we use time moments t 0 = 0, . . . , t i , . . . , t N = T for the time discretization, we denote the ALE displacement at time t i by u i a . Its NURBS representation can be written as where d i k are the corresponding control points. Given u i a , we can easily construct a NURBS parametrization Now G i f can be used to discretize the Navier-Stokes equations (8)(9) in the deformed configuration of the fluid domain.
As a final comment, let us mention that it is possible to deform only a portion of the fluid domain while keeping the rest fixed. An obvious advantage of this approach is a considerable reduction in computational cost associated with construction of ALE mappings. In this work, we choose to deform only three patches of the fluid domain that are directly adjacent to the FSI interface. In Figure 2, they are highlighted in red. In our experience, these three patches are more than capable of absorbing deformation of the structure that appears in this benchmark.

Mesh deformation techniques
The scope of this paper is to compare various techniques to construct ALE mappings. Put simply, an ALE mapping is nothing else but a deformation of the computational mesh in the fluid domain. In this section, we introduce several commonly used mesh deformation techniques (MDTs) as well as propose certain variations to them. All considered techniques achieve the same goal: provided a displacement of the FSI interface at a given time, they extend it into the fluid domain. Note that although we present these techniques in a 2D setting, one can readily apply them in 3D. When comparing different techniques, we largely pay our attention to the maximum amount of deformation a particular technique can handle. That is, how much the mesh can be deformed before the bijectivity condition (11) is violated. A secondary measure is of course the overall computational cost associated with computing ALE mappings using a given MDT.
After MDTs, we describe the Jacobian-based local stiffening [23] which we use to augment each of the considered techniques. Finally, we discuss practical ways to check whether the bijectivity condition (11) is satisfied.

Harmonic extension (HE/IHE)
Probably the simplest way to extend displacement of the FSI interface into the fluid domain is by means of harmonic extension (HE) [19,8,36]. Given the interface displacement u i s at time t i , the ALE displacement u i a is computed by solving Laplace's equation in the initial configuration of the fluid domain Ω 0 f for every displacement component: The interface displacement u i s serves as a Dirichlet boundary condition on the FSI interface Γ 0 . At the rest of the boundary, the prescribed displacement is zero. Note that the method does not take into account information about the interface or ALE displacement from the previous time step.
The HE technique is the least computationally expensive method to construct ALE mappings. Let N be the number of inner control points in the fluid domain. Since all displacement components satisfy the same equation, they can be computed by solving a single linear system with an N × N matrix and an N × d right-hand side. Here, d is the problem dimension; in our case, d = 2. Moreover, once the matrix is assembled, it can be reused for all time steps, which drastically reduces the computational cost associated with matrix assembly in IGA. At each time step, one only has to update the right-hand side to take the current FSI interface displacement into account. This can be performed efficiently by storing a Dirichlet elimination matrix when assembling the main matrix.
Despite its computational efficiency, the HE technique has serious disadvantages. First of all, it treats displacement components as completely independent variables and does not promote bijectivity of the ALE mapping in any way. Second, solutions of Laplace's equation in the vicinity of corners behave like r π/ω , where r is the distance to the corner and ω is the corner angle. For reentrant corners, that is for ω > π, solutions do not belong to H 1 (Ω 0 f ) since their derivatives tend to infinity. As a consequence, the corresponding ALE mapping may lose its bijectivity. Due to these two problems, the HE technique usually is only able to handle rather small deformations.
We propose a slight improvement to the HE technique achieved by turning it into an incremental algorithm. Assume that the ALE displacement u i a at time t i is known. We can use it to deform the initial configuration of the fluid domain Ω 0 f and obtain the deformed (19) describes. We can then compute an ALE displacement increment δu i+1 a by solving Laplace's equation in the deformed configuration: And finally, we define the ALE displacement u i+1 a at time t i+1 as u i a + δu i+1 a . Note that the resulting incremental harmonic extension (IHE) technique is not equivalent to the HE technique since ALE increments are computed in deformed configurations of the fluid domain.
One advantage of the IHE technique is that it uses the ALE displacement from the previous time step. Therefore, IHE can be expected to perform slightly better than the HE technique, meaning that it can handle larger deformations. On the other hand, each step of the IHE technique is formulated in a different configuration of the fluid domain than the previous one. As a result, the technique requires matrix assembly at every time step, which makes it more computationally expensive than HE.
In what follows, we apply the same ideas to other MDTs and consider both their non-incremental and incremental versions, which often share the same advantages and disadvantages as the HE and IHE techniques.

Bi-harmonic extension (BE/IBE)
The HE technique is weak when it comes to large deformations. To overcome this problem, one can search for the ALE displacement as a solution to the bi-harmonic equation: Solutions of the bi-harmonic equation are known to be more regular in comparison to Laplace's equation and do not have problems at reentrant corners [19]. However, this bi-harmonic extension (BE) technique is often dismissed as too computationally expensive. Indeed, in order to solve the bi-harmonic equation, one has two options: either use C 1 -conforming elements, which in IGA requires G 1 -continuity between patches [5,3]; or use mixed elements with an auxiliary variable q to replace the bi-harmonic equation with two Laplace's equations [4]: In this work, we consider only the latter option since it is easier to implement for multi-patch geometries. In our interpretation, the BE technique has the following form: given the interface displacement u i s at time t i , the ALE displacement u i a is computed by solving the following linear system in the initial configuration of the fluid domain Ω 0 f : The BE technique shares many similarities with HE which make both techniques very efficient: it does not use information from previous time steps; the linear system has to be assembled only once; the multiple-righthand-sides approach can be used to compute all displacement components at once. However, the resulting linear system is twice the size of the HE linear system: the matrix is of size 2N × 2N , and the right-hand size is of size 2N × d. Moreover, the linear system has a saddle-point structure, so specialized linear solvers are necessary to solve it efficiently.
Just like with the HE technique, we propose an incremental variation of the bi-harmonic extension (IBE). An increment δu i+1 a is computed by solving the following system in the deformed configuration Ω i f : ∇q · n = 0 on ∂Ω i f .
After that, the ALE displacement u i+1 a at time t i+1 is defined as u i a + δu i+1 a . The IBE technique requires matrix assembly at each time step but can potentially handle larger mesh deformations than the BE technique.

Linear elasticity (LE/ILE)
The next MDT we consider is based on the linear elasticity theory. It is widely used in FSI applications and belongs to the state-of-the-art in the field [1,23,32]. The core idea is to treat the fluid domain as an elastic body and to construct ALE displacement as a solution to the equations of linear elasticity: Here, σ σ σ a is the Cauchy stress tensor related to the linearized strain tensor ε ε ε a = (∇u a + ∇u T a )/2 by the Hooke's law: σ σ σ a = λ a tr(ε ε ε a )I + 2µ a ε ε ε a .
The Lamé parameters λ a and µ a depend on Young's modulus E a and Poisson's ratio ν a . Since we do not apply volumetric or surface force to the fluid domain, Young's modulus does not affect the resulting ALE displacement. On the other hand, Poisson's ratio is important because it regulates resistance of the fluid mesh to volumetric changes. A too high value (close to 0.5) would result in an almost incompressible behavior, which could lead to excessive distortion of the mesh elements and numerical instabilities. In contrast to that, a too low value (close to 0 or even negative) can reduce resistance of the fluid mesh to bijectivity violation. Therefore, we recommend choosing a value between 0.3 and 0.45. Unlike the techniques based on harmonic and biharmonic extension, the linear elasticity MDT is best known in its incremental version. That is, given the ALE displacement u i a at time t i , an increment δu i+1 a is computed by solving the linear elasticity equations in the deformed configuration of the fluid domain Ω i f : After that, the ALE displacement u i+1 a at time t i+1 is defined as u i a + δu i+1 a . We refer to this technique as incremental linear elasticity (ILE). The ILE technique is known for its robustness and an ability to withstand large mesh deformations. However, little to no research has been conducted to explain its superior behavior.
With respect to computational cost, the ILE technique involves solving a linear system with a matrix of size dN × dN and a right-hand size of size dN × 1. The linear system has to be reassembled at each time step. Note that the size of the linear system depends on a dimension of the problem. Therefore, it scales worse from 2D to 3D than the IHE and IBE techniques.
For the sake of completeness, let us also study a non-incremental version of the ILE technique. The ALE displacement u i a at time t i is computed by solving the equations of linear elasticity in the initial configuration of the fluid domain Ω 0 f : We call this the linear elasticity (LE) technique. Although one can only expect it to perform well for small deformations, it is rather computationally inexpensive.
Similarly to the HE and BE techniques, the LE technique requires matrix assembly only once and lets one reuse the matrix for every time step.

Nonlinear elasticity (TINE)
The last MDT we present in this paper is based on equations of nonlinear elasticity. The idea is to construct the ALE displacement at each time step as an approximate solution to the local balance equations of linear momentum where P a = F a S a . To ensure bijectivity of the ALE mapping, we use a logarithmic variation of the neo-Hookean material law where C a = F T a F a is the right Cauchy-Green strain tensor. Similarly to the LE and ILE techniques, the Lamé parameters λ a and µ a can be computed from Young's modulus E a and Poisson's ratio ν a , of which only Poisson's ratio affects the solution of equations (44).
Due to the term ln J a in the neo-Hookean law (45), any solution of equations (44) satisfies the bijectivity condition (11). This fact makes the MDT based on equations (44) unique since it explicitly enforces the bijectivity condition. Unfortunately, equations (44) are nonlinear, and an attempt to fully solve them at each time step would make the MDT prohibitively expensive. However, since the ALE mapping should possess certain regularity in time, it is possible to use a solution of equations (44) at time t i to efficiently construct an approximate solution at time t i+1 [22]. We refer to this technique as tangential incremental nonlinear elasticity (TINE). The TINE technique can be seen as pseudo time-stepping or an example of the continuation methods for nonlinear problems [7].
Let us look under the hood of TINE. It is based on a Newton-like linearization of equations (44). To define it, we need to transform equations (44) into a weak form, also known as variation formulation. To that end, let us define a solution space V = (H 1 (Ω 0 f )) d and a test space We can then write the weak form of equations (44) as Here, δE a [w] = 1 2 F T a ∇w + ∇w T F a is the variation of the Green-Lagrange strain tensor. The Taylor expansion at point (u a , w) with an increment δu a yields R(u a + δu a , w) =R(u a , w)+ DR(u a , w) · δu a + o(||δu a ||), where DR(u a , w) · δu a is a directional derivative. We refer to [35,22] for details on computing DR(u a , w) · δu a . The idea of the TINE technique is to use one Newtonlike step find δu a ∈ V such that ∀w ∈ V 0 DR(u a , w) · δu a = −R(u a , w) per time step to compute an ALE increment and update the ALE displacement. Concretely, given the ALE displacement u i a at time t i , we find an ALE increment δu i+1 a as a solution of the linear problem After that, we define the ALE displacement u i+1 a at time t i+1 as u i a + δu i+1 a . It is natural to compare the TINE and ILE techniques which are very similar at first glance. Both are incremental techniques based on the elasticity theory; both require solution of a linear system with a matrix of size dN × dN and a right-hand side of size dN × 1; both require matrix assembly at each time step. In general, one can expect both techniques to be roughly equal in computational cost. Unlike ILE, however, the TINE technique explicitly enforces the bijectivity condition (11). Moreover, the TINE technique is based in the initial configuration of the fluid domain. As we show in Sections 4 and 5, this last observation results in crucial differences in behavior of the ILE and TINE techniques when it comes to the accumulated distortion effect.

Local stiffening
Most of the fluid mesh deformation happens along the FSI interface, where the structure displacement is applied as a Dirichlet boundary conditions to the ALE displacement. On the other hand, mesh elements in the vicinity of the stationary boundary of the fluid domain undergo almost no deformation. Therefore, their contribution into processing of the applied interface displacement is negligible. If the deformation could be redistributed away from the FSI interface towards the stationary boundary, the mesh could undergo larger deformations without becoming invalid. This is the idea behind local stiffening, which locally changes the way different elements react to the deformation.
Let G : [0, 1] d → Ω be a parametrization of the computational domain Ω. Imagine that we have to compute integrals corresponding to matrix entries of the discretized linear system. One of the simplest ways to implement local stiffening is to drop the Jacobian determinant det ∇G when transforming the integrals from domain Ω to parametric domain [0, 1] d : This local stiffening method was first proposed in [30]. For elasticity problems, (52) can be interpreted as a local change of Young's modulus which makes elements with small values of det ∇G stiffer and elements with large values softer. Therefore, the former elements undergo less deformation and are less likely to become invalid. A more advanced local stiffening method introduced in [23] does not simply drop the Jacobian determinant but changes the degree with which it enters the integrals: Here, χ 0 is called the stiffening degree. The higher the stiffening degree is, the more local stiffening is achieved. χ = 0 corresponds to no local stiffening, and χ = 1 corresponds to Jacobian dropping (52). Too high stiffening degrees, however, may result in excessive mesh distortion. We refer to this method as the Jacobian-based local stiffening. The Jacobian-based local stiffening acts differently depending on whether the mesh deformation method is formulated in the initial or in the deformed configuration of the fluid domain. Namely, if the integrals for matrix entries are computed in the initial configuration Ω 0 f , the local stiffening is based only on the initial parametrization G 0 f . However, in the case of the deformed configuration Ω i f , the local stiffening takes into account already applied deformation since the parametrization G i f of the deformed configuration is defined as (I + u i a ) • G 0 f , see equation (19). This effect has both advantages and disadvantages. From one point of view, if a particular mesh element becomes ill-shaped after the deformation, its value of det ∇G i f decreases. As a result, this element receives more stiffening, which prevents it from becoming even more ill-shaped or invalid. On the other hand, in case of the ILE technique, this deformation-aware local stiffening essentially makes material properties of the mesh deformation-dependent, which can cause irreversible plastic deformation. For other MDTs based in the deformed configuration, namely IHE and IBE, the effect is similar. As we show in Sections 4 and 5, this irreversible deformation accumulates over time and can significantly affect results of FSI simulations. We refer to this effect as accumulated distortion.
Regardless if the local stiffening is deformation-aware or not, the initial parametrization G 0 f of the fluid domain provides a major contribution to how much stiffening each element receives. Let us consider the deforming part of the fluid domain, see Figure 3. The top and bottom patches are perfect parallelograms, and det ∇G 0 f is constant. Therefore, elements of these patches receive no local stiffening with respect to each other. However, the right patch has a distinct tapered left side, where det ∇G 0 f becomes very small in comparison to the surrounding elements of the top and bottom patches. As a result, element on the left side of the right patch become much stiffer and maintain their shape. Consequently, angles of all three mesh patches that are adjacent to the beam right end do not change much during mesh deformation. In particular, they do not exceed π, which would lead to det ∇G i f becoming negative, which means that the bijectivity condition (11) is not violated. Note that local stiffening may be harder to achieve in IGA than in classical FEM. We have seen an example of a parallelogram with a natural uniform tensorproduct NURBS parametrization. In that case, det ∇G is constant throughout the geometry, and the Jacobianbased local stiffening would have no effect. With a finite element discretization of a parallelogram, it would be possible to place smaller elements where necessary. Since mappings to a reference element are independent for each FEM element, the Jacobian-based local stiffening would have an effect. In IGA, it is possible to locally control the element size distribution with local refinement, for example by means of THB-splines [34]. However, local refinement does not change the underlying geometry parametrization and values of the Jacobian determinant. Therefore, in order to apply the Jacobian-based local stiffening in IGA, one has two options: either use non-uniform parametrizations with carefully designed knot vectors to artificially construct domain regions with smaller values of the Jacobian determinant; or design patch geometries in a way that naturally defines such domain regions. An example in Figure 3 belongs to the latter approach.

Bijectivity check
Let us briefly discuss ways to check the bijectivity condition (11) in practice. A solution which takes the NURBS nature of the ALE displacement u a into account is to express J a as a NURBS function [10]. If all coefficients in a NURBS representation of J a are positive, then the displacement u a satisfies the bijectivity condition. Unfortunately, this condition is only sufficient and not a necessary one. Therefore, it may often lead to false detection of bijectivity violation. In practice, we resort to a less elegant solution of sampling J a at the Gaussian quadrature points associated with the NURBS basis of u a .
Note that whichever method is chosen, it introduces a certain computational overhead to construction of ALE mappings. Nevertheless, we recommend doing some bijectivity check at every time step, or at least with regular intervals. An ALE mapping that does not satisfy the bijectivity condition (11) quickly makes all ensuing computation results meaningless.

Benchmark ALE: mesh deformation
In order to test and compare all mesh deformation techniques (MDTs), we first consider a simplified FSI-like test based on the benchmark introduced in Section 2. Instead of solving a fully coupled FSI problem, we ignore the fluid component and let the flexible beam oscillate freely in the presence of external acceleration g = (0, l). We use the beam motion to drive deformation of the three fluid domain patches adjacent to the beam. By varying the parameter l, we can regulate magnitude of the mesh deformation. Although this mesh deformation test is artificial, it mimics real deformations occurring in the original benchmark well enough. Moreover, it is significantly less computationally expensive, which allows us to conduct more tests and better assess properties of each MDT. For this test, we refine the fluid domain parametrization by applying uniform h-refinement thrice. Figure  4 shows the corresponding computational mesh in the state of maximum deformation for loading levels l = 1, 2, 3. We use the following parameters for the structure motion: Young's modulus E s = 1.4 × 10 6 kg·m −1 ·s −2 , Poisson's ratio ν s = 0.4, density ρ s = 10 3 kg·m −3 . For MDTs based on the elasticity theory, we use Poisson's ratio ν a = 0.3. At each time step, we check whether the bijectivity condition (11) holds. We have implemented this test and the original FSI benchmark within G+Smo-an open-source C++ library for isogeometric analysis [15]-using the gsElasticity submodule. As a linear system solver, we use Pardiso-an efficient parallel direct linear solver [16]. For reference, all simulations have been performed on a laptop with a 7th generation Intel Core i7 CPU using eight hyper-threads.

Single period test
First, we consider mesh deformation over one period of beam oscillations. The goal is to study how much deformation each MDT can handle. Here, local stiffening is of crucial importance. Without it, most MDTs can handle only small loading levels l. Usually, one of the patch corners adjacent to the right end of the beam becomes larger than π, which violates the bijectivity conditions (11). With local stiffening, all MDTs can keep these angles below π at least for moderate loading levels l, see Figure 5. Figure 6 shows a plot of the maximum achievable loading level l max versus the stiffening degree χ for each MDT. Immediately, we can split all MDTs into two groups, with BE and IBE forming one group, and all other techniques belonging to the second group. The main difference between the two groups is that MDTs from the second group can handle almost no deforma- tion without local stiffening. However, as the stiffening degree χ grows, these MDTs can handle increasingly larger loading levels with maximum loading levels achieved with χ ∈ [2,3]. With χ > 3, we can observe some form of performance deterioration for all MDTs, which is likely caused by too much mesh distortion introduced by the local stiffening.
For the BE and IBE techniques, the behavior is radically different. Already without local stiffening, they can handle larger loading levels than some MDTs from the second group can achieve even with local stiffening. However, as we introduce local stiffening, BE and IBE show almost no response for χ < 1 and start to slowly perform worse for χ > 1.
Overall, MDTs can be ranged with respect to their capability to handle large deformations in the following way: IBE, ILE and TINE are the most powerful and can handle loading levels up of to 2.8-3; BE occupies the second place with the maximum loading level of 2.2; and HE, IHE and LE are the least powerful with maximum loading levels of 1.6-1.9. Figure 7 presents an analysis of computational complexity of each MDT. We have measured time required for linear system assembly, linear system solution and an ensuing check of the bijectivity condition (11). Note that this is real time and not CPU time. Although not a perfect measure of algorithm performance, real time still allows us to compare relative efficiency of different MDTs since we have implemented them in the same framework of G+Smo.
The time analysis shows that the HE, BE and LE techniques are significantly faster than their incremental versions. This result is not surprising because nonincremental techniques do not require matrix assembly at each time step. The BE and IBE techniques take the largest amount of time to solve the linear system due to the saddle-point structure of the system. At the same time, relative complexity of the elasticity equations makes the assembly time of the ILE and TINE techniques significantly larger than for other

Long-term behavior test
In the second test, we study the long-term behavior of different MDTs and their effect on the fluid mesh. To that end, we perform the simulation over a time period of 20s, which includes roughly 22 periods of beam oscillations. A quantity of interest is the L 2 -norm of the ALE displacement measured in the initial configuration of the fluid domain. A perfect MDT should return the fluid mesh to its initial state once the beam is not deformed. Therefore, the ALE norm ||u a (t)|| L 2 (Ω 0 f ) should be close to zero at the end of each oscillation period. Figure 8 shows behavior of the ALE norm over time for each MDT with the loading level l = 1.5. We have used Jacobian-based local stiffening with χ = 2 for all MDTs with the exception of the BE and IBE techniques. These techniques are able to handle the loading level l = 1.5 without local stiffening.
As Figure 8 shows, the ALE norm behaves periodically and returns to zero at the end of each oscillation period with the HE, BE, LE and TINE techniques. On the other hand, with the IHE, IBE and ILE techniques the ALE norm at the end of each oscillation period grows in a monotonous fashion. This effect has been previously reported in [28], and we refer to it as accumulated distortion. It appears only for MDTs which are based in the deformed configuration of the fluid domain. As a result, mesh deformation becomes pathdependent, the fluid mesh does not return to its initial state, and its quality deteriorates over time.  illustrates the state of the mesh at the end of the simulation with the ILE technique. The accumulated distortion effect becomes more prominent as the mesh deformation magnitude grows. To study it in more details, we perform the long-term behavior test for the IHE, IBE and ILE techniques with varying values of the loading level and stiffening degree. In Figure 9, we plot values of the ALE norm at the end of each oscillation period. We can observe that both parameters seem to increase the rate of accumulated distortion.

Benchmark FSI2: flow-induced vibrations
In this section, we perform the FSI simulation described in Section 2. From several simulation scenarios proposed in [33], we choose a scenario titled FSI2 since it corresponds to the largest magnitude of beam displacement and mesh deformation. Using this FSI simulation, we test and compare different mesh deformation techniques (MDTs) introduced in Section 3. We use stiffening degree χ = 2.5 for all MDTs but BE and IBE. To them, we apply no local stiffening. For the analysis, we refine the geometry parametrization five times using uniform h-refinement and perform the simulation for 15 seconds with a time step ∆t = 0.0025 s. For time integration, we use the Newmark method [35] with β = 0.5 and γ = 1 for the structure and the IMEX scheme [13] with θ = 0.5 and no stabilization for the fluid. We achieve the coupling of fluid and structure by means of the partitioned Fluid-Dirichlet-Structure-Neumann algorithm [8].
The FSI2 scenario is characterized by the following parameters: fluid density ρ f = 10 3 kg·m −3 , fluid kinetic viscosity ν f = 10 −3 m 2 ·s −1 , maximum inflow velocity v max = 1.5 m·s −1 , structure density ρ s = 10 4 kg·m −3 , structure Young's modulus E = 1.4 × 10 6 kg·m −1 ·s −2 , structure Poisson's ratio ν s = 0.4, gravitational acceleration g = (0, 0) T m·s −2 and mesh Poisson's ratio ν a = 0.3 (where applicable). The corresponding Reynolds number is Re = 100, which results in an unstable flow and development of vortex shedding. Alternating downward and upward forces exerted on the structure by the fluid lead to oscillations of the beam which grow in magnitude until they reach a fully periodic regime. Figure 10 illustrates typical fluid velocity field and beam deformation when the oscillations are fully developed.
To assess the simulation accuracy, we study the following quantities of interest when the beam oscillations are fully developed : x-and y-displacement components u x s (A) and u y s (A) of the point A located in the middle of the beam right end; and drag and lift forces F D and F N exerted on the structure by the fluid which are defined as Here, Σ(t) is the entire boundary of the submerged solid, including the rigid disk and the flexible beam, at time t. Since we expect the quantities of interest to behave periodically, we report them in terms of their mean ( * ) max + ( * ) min /2, amplitude ( * ) max − ( * ) min /2 and frequency. In Table 1, we compare the simulation results obtained with the TINE technique against the reference results from [33]. Overall, our results seem to undershoot the reference values by about 5%, which can be expected since we use much fewer degrees of freedom and a larger time step than the reference simulation. Despite this discrepancy, we are more than capable of reproducing a qualitatively correct behavior of the system and can use it to study the MDTs. Figures 11 and  12 illustrate behavior of the lift, drag and beam displacement.
Let us now focus on the fluid mesh deformation. When performing the FSI simulation, we apply each of the seven MDTs (HE, IHE, BE, IBE, LE, ILE and TINE) and study how a particular MDT handles mesh deformations occurring during the simulation. Figure  13 illustrates what portion of the simulation can be completed using different MDTs. As we can observe, simulations with the HE, IHE and IBE techniques had to be terminated before they could reach the end. All three techniques have failed to maintain bijectivity of the ALE mapping; however, different reasons have led to the failure. In the case of HE and IHE, the simulations have stopped at the ninth second-when the oscillations in the system start to develop. As the applied mesh deformation grows, the HE and IHE techniques fail due to their intrinsic inability to handle large deformations.
On the other hand, the IBE technique is able to handle mesh deformations occurring in the simulation but suffers from the accumulated deformation effect described in Section 4. As a result, the simulation fails at the 14th second. Figure 14 depicts behavior of the ALE norm for all MDTs. We can observe that the ILE technique suffers from even stronger accumulated deformation than IBE but still manages to maintain a bijective ALE mapping until the simulation end. Unfortunately,   the highly distorted mesh produced by the ILE technique (see Figure 15) affects the simulation results, see Figures 11 and 12. Instead of a stable periodic behavior in the fully developed oscillation regime, we observe signs of damping. We observe a similar effect when using the IBE technique: the accumulated mesh deformation results in spurious amplification of the oscillations.
Of the seven MDTs considered in this work, only three-BE, LE and TINE-were able to handle mesh deformations occurring in the simulation and maintain high mesh quality until the simulation end. Most importantly, the BE, LE and TINE techniques have demonstrated no signs of accumulated distortion. Using these techniques, we were able to reproduce a stable periodic behavior of the system and correct simulation results. Figure 16 depicts the fluid mesh at the end of the simulation deformed using the TINE technique. Of course, the BE, LE and TINE techniques differ a lot in terms of their computational cost. However, since construction of the ALE mapping corresponds only to a small portion of a total computational effort required to perform an FSI simulation, the choice of MDT does not affect the total computational cost too much. Table 2 compares computational cost of BE, LE and TINE against the ILE technique which is often considered a default option in the FSI community.

Discussion and conclusion
In this work, we have described and compared several mesh deformation techniques (MDTs) which can be used within moving-mesh methods for FSI problems. To evaluate each MDT, we have used a 2D FSI benchmark and its simplified version where the focus lies on mesh deformation. Based on the tests performed in Sections 4 and 5, we can make the following conclusions: -Out of seven MDTs that we have considered, two most robust are bi-harmonic extension (BE) and tangential incremental nonlinear elasticity (TINE). Both MDTs can handle large mesh deformations and do not suffer from the accumulated distortion effect. -BE is easier to implement, performs well even without the Jacobian-based local stiffening and is about two times less computationally expensive than TINE. Provided that the saddle-point structure of the linear system is accounted for, we recommend the BE technique as the first method to try in many FSI applications.    We would like to emphasize that the performance of all MDTs is dependent on the chosen parametrization of the fluid domain. However, it is unlikely that the MDT behavior will be qualitatively different if a different geometry parametrization is used. With respect to the further research directions, we see the following possibilities. One could study the effect of different geometry parametrizations on the MDT behavior. Moreover, isogeometric methods in FSI could benefit from alternative local stiffening approaches. The commonly used Jacobian-based local stiffening has no effect if a uniform geometry parametrization is used. Finally, it would be interesting to combine BE and TINE with additional augmentation techniques such as the solid layer extension and mesh element relaxation.