A low-distortion mesh moving method based on ﬁber-reinforced hyperelasticity and optimized zero-stress state

In computation of ﬂow problems with moving boundaries and interfaces, including ﬂuid–structure interaction, moving-mesh methods enable mesh-resolution control near the interface and consequently high-resolution representation of the boundary layers. Good moving-mesh methods require good mesh moving methods. We introduce a low-distortion mesh moving method based on ﬁber-reinforced hyperelasticity and optimized zero-stress state (ZSS). The method has been developed targeting isogeometric discretization but is also applicable to ﬁnite element discretization. With the large-deformation mechanics equations, we can expect to have a unique mesh associated with each step of the boundary or interface motion. With the ﬁbers placed in multiple directions, we stiffen the element in those directions for the purpose of reducing the distortion during the mesh deformation. We optimize the ZSS by seeking orthogonality of the parametric directions, by mesh relaxation, and by making the ZSS time-dependent as needed. We present 2D and 3D test computations with isogeometric discretization. The computations show that the mesh moving method introduced performs well.


Introduction
For computation of flow problems with moving boundaries and interfaces (MBI), including fluid-structure interaction (FSI), we introduce a low-distortion mesh moving method based on fiber-reinforced hyperelasticity and optimized zerostress state (ZSS). With the large-deformation mechanics equations, we can expect to have a unique mesh associated with each step of the boundary motion. The method has been developed targeting isogeometric discretization but is also applicable to finite element discretization. To show how the

Key features
In moving-mesh (interface-tracking) methods, as the interface moves and the fluid mechanics domain changes its shape, the mesh moves to adjust to the shape change and to follow (i.e. "track") the interface. Moving the fluid mechanics mesh to follow a fluid-solid interface enables us to control the mesh resolution near the interface, have highresolution representation of the boundary layers, and obtain accurate solutions in such critical flow regions. These desirable features do not come easily or do not come at all with the nonmoving-mesh (interface-capturing) methods. In these methods, the interface geometry is somehow represented over a nonmoving fluid mechanics mesh, with more accuracy in some methods than in some others, but the key point is that the fluid mechanics mesh does not move to follow the interface. Because the mesh is not following the interface, independent of how accurately the interface geometry is rep-resented, the boundary layer resolution will be limited by the fluid mechanics mesh resolution where the interface is.

Mixed moving-mesh/nonmoving-mesh methods
As mentioned in [1], one of course recognizes that certain classes of interfaces (such as free-surface and two-fluid flows with splashing) might be too complex to handle with an interface-tracking technique and, therefore, for all practical purposes, require an interface-capturing technique. The Mixed Interface-Tracking/Interface-Capturing Technique (MITICT) [2] was introduced in 2000 for computation of MBI problems that involve both fluid-solid interfaces that can be accurately tracked with a moving-mesh method and fluid-fluid interfaces that are too complex or unsteady to be tracked. Such fluid-fluid interfaces are captured over the mesh tracking the fluid-solid interfaces. Thinking similarly about MBI problems with an actual contact or other topology change, the Fluid-Solid Interface-Tracking/Interface-Capturing Technique (FSITICT) was introduced in [3] as the FSI version of the MITICT. In the FSITICT, we track the interface we can with a moving mesh, and capture over that moving mesh the interfaces we cannot track, specifically the interfaces where we need to have an actual contact between the solid surfaces.

Space-time (ST) methods
The Deforming-Spatial-Domain/Stabilized ST (DSD/SST) method [4][5][6], introduced for computation of flows with MBI, including FSI, is a moving-mesh method, possessing the associated desirable features. The stabilization components of this original DSD/SST are the Streamline-Upwind/Petrov-Galerkin (SUPG) [7] and Pressure-Stabilizing/Petrov-Galerkin (PSPG) [4] stabilizations, and therefore it is now called "ST-SUPS." The ST Variational Multiscale (ST-VMS) method [8][9][10], which subsumes its precursor ST-SUPS is the VMS version of the DSD/SST. The VMS components of the ST-VMS are from the residual-based VMS (RBVMS) method [11][12][13][14]. The ST-VMS has two more stabilization terms beyond those in the ST-SUPS, and the additional terms give the method better turbulence modeling features. The ST-SUPS and ST-VMS, because of the higher-order accuracy of the ST framework (see [8,9]), are desirable also in computations without MBI.

ST-SI
The ST-SI was introduced in [86], in the context of incompressible-flow equations, to retain the desirable movingmesh features of the ST-VMS and ST-SUPS in computations involving spinning solid surfaces, such as a turbine rotor. The mesh covering the spinning surface spins with it, retaining the high-resolution representation of the boundary layers, while the mesh on the other side of the SI remains unaffected. This is accomplished by adding to the ST-VMS formulation interface terms similar to those in the version of the ALE-VMS for computations with sliding interfaces [26,27]. The interface terms account for the compatibility conditions for the velocity and stress at the SI, accurately connecting the two sides of the solution. An ST-SI version where the SI is between fluid and solid domains was also presented in [86]. The SI in that case is a "fluid-solid SI" rather than a standard "fluid-fluid SI" and enables weak enforcement of the Dirichlet boundary conditions for the fluid. The ST-SI introduced in [110] for the coupled incompressible-flow and thermal-transport equations retains the high-resolution representation of the thermo-fluid boundary layers near spinning solid surfaces. These ST-SI methods have been applied to aerodynamic analysis of vertical-axis wind turbines [37,38,86], thermo-fluid analysis of disk brakes [110], flow-driven string dynamics in turbomachinery [111][112][113], flow analysis of turbocharger turbines [114][115][116][117][118], flow around tires with road contact and deformation [103,[119][120][121][122], fluid films [122,123], aerodynamic analysis of ram-air parachutes [124], and flow analysis of heart valves [99,101,[104][105][106] and ventricle-valve-aorta sequences [107].
In the ST-SI version presented in [86] the SI is between a thin porous structure and the fluid on its two sides. This enables dealing with the porosity in a fashion consistent with how the standard fluid-fluid SIs are dealt with and how the Dirichlet conditions are enforced weakly with fluid-solid SIs. This version also enables handling thin structures that have T-junctions. This method has been applied to incompressibleflow aerodynamic analysis of ram-air parachutes with fabric porosity [124]. The compressible-flow ST-SI methods were introduced in [125], including the version where the SI is between a thin porous structure and the fluid on its two sides. Compressible-flow porosity models were also introduced in [125]. These, together with the compressible-flow ST SUPG method [127], extended the ST computational analysis range to compressible-flow aerodynamics of parachutes with fabric and geometric porosities. That enabled ST computational flow analysis of the Orion spacecraft drogue parachute in the compressible-flow regime [125,126].

ST-TC
The ST-TC [92,102] was introduced for moving-mesh computation of flow problems with TC, such as contact between solid surfaces. Even before the ST-TC, the ST-SUPS and ST-VMS, when used with robust mesh update methods, have proven effective in flow computations where the solid surfaces are in near contact or create other near TC, if the nearness is sufficiently near for the purpose of solving the problem. Many classes of problems can be solved that way with sufficient accuracy. For examples of such computations, see the references mentioned in [92]. The ST-TC made moving-mesh computations possible even when there is an actual contact between solid surfaces or other TC. By collapsing elements as needed, without changing the connectivity of the "parent" mesh, the ST-TC can handle an actual TC while maintaining high-resolution boundary layer representation near solid surfaces. This enabled successful moving-mesh computation of heart valve flows [85,92,99,[101][102][103][104][105][106], ventricle-valve-aorta flows [107], wing clapping [92,93], flow around a rotating tire with road contact and prescribed deformation [103,[119][120][121][122], and fluid films [122,123].

ST-SI-TC
The ST-SI-TC is the integration of the ST-SI and ST-TC. A fluid-fluid SI requires elements on both sides of the SI. When part of an SI needs to coincide with a solid surface, which happens for example when the solid surfaces on two sides of an SI come into contact or when an SI reaches a solid surface, the elements between the coinciding SI part and the solid surface need to collapse with the ST-TC mechanism. The collapse switches the SI from fluid-fluid SI to fluidsolid SI. With that, an SI can be a mixture of fluid-fluid and fluid-solid SIs. With the ST-SI-TC, the elements collapse and are reborn independent of the nodes representing a solid surface. The ST-SI-TC enables high-resolution flow representation even when parts of the SI are coinciding with a solid surface. It also enables dealing with contact location change and contact sliding. This was applied to heart valve flow analysis [99,101,[104][105][106], ventricle-valve-aorta flows [107], tire aerodynamics with road contact and deformation [103,[119][120][121][122], and fluid films [122,123].

ST-IGA
The success with IGA basis functions in space [16,26,44,128] motivated the integration of the ST methods with isogeometric discretization, which we broadly call "ST-IGA." The ST-IGA was introduced in [8]. Computations with the ST-VMS and ST-IGA were first reported in [8] in a 2D context, with IGA basis functions in space for flow past an airfoil, and in both space and time for the advection equation. Using higher-order basis functions in time enables deriving full benefit from using higher-order basis functions in space. This was demonstrated with the stability and accuracy analysis given in [8] for the advection equation.
The ST-IGA with IGA basis functions in time enables a more accurate representation of the motion of the solid surfaces and a mesh motion consistent with that. This was pointed out in [8,9] and demonstrated in [87,88,90]. It also enables more efficient temporal representation of the motion and deformation of the volume meshes, and more efficient remeshing. These motivated the development of the ST/NURBS Mesh Update Method (STNMUM) [87,88,90], with the name coined in [83]. The STNMUM has a wide scope that includes spinning solid surfaces. With the spinning motion represented by quadratic NURBS in time, and with sufficient number of temporal patches for a full rotation, the circular paths are represented exactly. A "secondary mapping" [1,8,9,87] enables also specifying a constant angular velocity for invariant speeds along the circular paths. The ST framework and NURBS in time also enable, with the "ST-C" method, extracting a continuous representation from the computed data and, in large-scale computations, efficient data compression [10,103,[110][111][112][113]129]. The STN-MUM and the ST-IGA with IGA basis functions in time have been used in many 3D computations. The classes of problems solved are flapping-wing aerodynamics for an actual locust [1,[87][88][89], bioinspired MAVs [84,85,90,91] and wing-clapping [92,93], separation aerodynamics of spacecraft [77], aerodynamics of horizontal-axis [34,[83][84][85] and vertical-axis [37,38,86] wind-turbines, thermo-fluid analysis of ground vehicles and their tires [10,103], thermo-fluid analysis of disk brakes [110], flow-driven string dynamics in turbomachinery [111][112][113], flow analysis of turbocharger turbines [114][115][116][117][118], and flow analysis of coronary arteries in motion [108].

ST-SI-IGA and ST-SI-TC-IGA
The ST-SI-IGA is the integration of the ST-SI and ST-IGA, and the ST-SI-TC-IGA is the integration of the ST-SI, ST-TC and ST-IGA.
The turbocharger turbine flow [114][115][116][117][118] and flow-driven string dynamics in turbomachinery [112,113] were computed with the ST-SI-IGA. The IGA basis functions were used in the spatial discretization of the fluid mechanics equations and also in the temporal representation of the rotor and spinningmesh motion. That enabled accurate representation of the turbine geometry and rotor motion and increased accuracy in the flow solution. The IGA basis functions were used also in the spatial discretization of the string structural dynamics equations. That enabled increased accuracy in the structural dynamics solution, as well as smoothness in the string shape and fluid dynamics forces computed on the string.
The ram-air parachute analysis [124] and spacecraft parachute compressible-flow analysis [126] were conducted with the ST-SI-IGA, based on the ST-SI version that weakly enforces the Dirichlet conditions and the ST-SI version that accounts for the porosity of a thin structure. The ST-IGA with IGA basis functions in space enabled, with relatively few number of unknowns, accurate representation of the parafoil and parachute geometries and increased accuracy in the flow solution. The volume mesh needed to be generated both inside and outside the parafoil. Mesh generation inside was challenging near the trailing edge because of the narrowing space. The spacecraft parachute has a very complex geometry, including gores and gaps. Using IGA basis functions addressed those challenges and still kept the element density near the trailing edge of the parafoil and around the spacecraft parachute at a reasonable level.
The heart valve flow analyses [99,101,[104][105][106] and ventricle-valve-aorta flow analysis [107] were conducted with the ST-SI-TC-IGA. The method, beyond enabling a more accurate representation of the geometry and increased accuracy in the flow solution, kept the element density in the narrow spaces near the contact areas at a reasonable level. When solid surfaces come into contact, the elements between the surface and the SI collapse. Before the elements collapse, the boundaries could be curved and rather complex, and the narrow spaces might have high-aspect-ratio elements. With NURBS elements, it was possible to deal with such adverse conditions rather effectively.
In computational analysis of flow around tires with road contact and deformation [120][121][122], the ST-SI-TC-IGA enables a more accurate representation of the geometry and motion of the tire surfaces, a mesh motion consistent with that, and increased accuracy in the flow solution. It also keeps the element density in the tire grooves and in the narrow spaces near the contact areas at a reasonable level. In addition, we benefit from the mesh generation flexibility provided by using SIs.
In computational analysis of fluid films [122,123], the ST-SI-TC-IGA enables solution with a computational cost comparable to that of the Reynolds-equation model for the comparable solution quality [123]. With that, narrow gaps associated with the road roughness [122] can be accounted for in the flow analysis around tires.
An SI provides mesh generation flexibility in a general context by accurately connecting the two sides of the solution computed over nonmatching meshes. This type of mesh generation flexibility is especially valuable in complex-geometry flow computations with isogeometric discretization, removing the matching requirement between the NURBS patches without loss of accuracy. This feature was used in the flow analysis of heart valves [99,101,104,105], ventricle-valveaorta sequences [107], turbocharger turbines [114][115][116][117][118], and spacecraft parachute compressible-flow analysis [126].

Moving-mesh methods are worth the effort to make them work
It is quite evident from Sects. 1.1.4, 1.1.5 and 1.2 that movingmesh methods have been practical in more classes of complex FSI and MBI problems than commonly thought of, and, with the increased scope provided by the special methods, the ST methods can now do even more. In many classes of computations, the moving-mesh methods are worth the effort to make them work.

Mesh moving methods
Good moving-mesh methods require good mesh moving methods, and that has to be part of the effort to make the moving-mesh methods work.

Mesh moving and remeshing
In a more general context, moving-mesh methods require mesh update methods. Mesh update has two components: moving the mesh for as long as it is possible, which is the core component, and full or partial remeshing (i.e. generating a new set of elements, and sometimes also a new set of nodes) when the element distortion becomes too high. The key objectives of a mesh moving method should be to maintain the element quality near solid surfaces and to minimize remeshing frequency. Provided that at the boundaries and interfaces the normal velocity of the mesh matches the normal velocity of the fluid, we can move the mesh any way we find most suitable for meeting those objectives.

Mesh moving based on linear elasticity and mesh-Jacobian-based stiffening
In the mesh moving method introduced in 1992 [ [137][138][139], the motion of the internal nodes is determined by solving the equations of linear elasticity. As the boundary condition, the normal velocity of the mesh at the boundaries and interfaces is specified to match the velocity of the fluid. The mesh deformation is dealt with selectively based on the sizes of the elements. Selective treatment based on element sizes is attained by altering the way we account for the Jacobian of the transformation from the element domain to the physical domain. The objective is to stiffen the smaller elements, which are typically placed near solid surfaces, more than the larger ones. When this method was first introduced in [137][138][139], it consisted of simply dropping Jacobian from the finite element formulation of the mesh moving (elasticity) equations. This results in the smaller elements being stiffened more than the larger ones. The method was given the name "Jacobian-based stiffening" (JBS) in [140]. It was also augmented in [140] to a more extensive kind by introducing a stiffening power χ that determines the degree by which the smaller elements are rendered stiffer than the larger ones. This approach, when χ = 1, would be identical to the method first introduced in [137], and most of the time χ = 1. To also clarify that the "Jacobian" in the name of the method is the mesh Jacobian, in [107] the method was renamed "Mesh-Jacobian-based stiffening" (MJBS).

Solid-Extension Mesh Moving Technique (SEMMT)
In dealing with fluid-solid interfaces, sometimes we need to place thin layers of elements near the solid surfaces to fully control the mesh resolution there and have more accurate representation of the boundary layers. In the mesh moving method introduced in [137][138][139], such layers of elements move "glued" to the solid objects, undergoing rigid-body motion. No equations are solved for the motion of the nodes in these layers, because these nodal displacements are not governed by the equations of elasticity. This results in some cost reduction. But, more importantly, we have full control of the mesh resolution in these layers. Examples of automatic mesh moving combined with thin layers of elements undergoing rigid-body motion with solid objects were reported as early as 1993 (see [138,139]). Earlier examples of element layers undergoing rigid-body motion, in combination with deforming structured meshes, were reported in 1990 [4].
In computation of flows with fluid-solid interfaces where the solid is deforming, the motion of the fluid mesh near the interface cannot be represented by a rigid-body motion. Depending on the deformation mode of the solid, the mesh moving method described in Sect. 1.4.2 may need to be used. In such cases, the thin layers of elements near the solid surface become a challenge for the mesh moving method. The SEMMT was proposed in 2001 [141]. In the SEMMT [1,[141][142][143][144], the thin layers of elements are treated almost like an extension of the solid elements. In solving the equations of elasticity governing the motion of the fluid nodes, higher stiffness is assigned to the thin layers of elements compared to the other fluid elements. Two ways of accomplishing this were proposed in [141]: solving the elasticity equations for the nodes connected by the thin layers of elements separate from the elasticity equations for the other nodes, or together. If they are solved separately, for the thin layers of elements, as boundary conditions at the interface with the other elements, traction-free conditions would be used. The separate solution option is referred to as "SEMMT -Multiple Domain (SEMMT-MD)," and the unified solution option as "SEMMT -Single Domain (SEMMT-SD)." In [1,[143][144][145], test computations were presented to demonstrate how the SEMMT works as part of the mesh update method used with the ST-SUPS method. Both SEMMT options described above were employed. The tests included mesh deformation tests [1,[143][144][145] and a 2D FSI model problem [143,144].

Mesh moving based on large-deformation mechanics
In computing the mesh motion from time level t n to t n+1 with the linear-elasticity equations, typically the configuration at t n is where we compute the displacement from. That works for most problems, but, as mentioned in [107], there are two closely related drawbacks: i) the method is pathdependent, ii) once elements accumulate in some region, it is hard to undo that. The path-dependence leads to noncyclic results in FSI computations that we expect to have cyclic or near-cyclic results. While the "back-cycle-based mesh moving" (BCBMM) method introduced in [107] can remedy that, moving the mesh based on large-deformation mechanics equations provides a more general solution. That was the mesh moving method used in a number of recent FSI and other MBI computations with the ST-SUPS and ST-VMS (see, for example, [10,75,85,92,97]). Additional good features of this method are, as pointed out in [75], having many constitutive models to choose from and being able to define the ZSS locally in arbitrary orientations. The constitutive model was St. Venant-Kirchhoff in [75], and neo-Hookean in [10,85,92,97]. The MJBS can still be part of the mesh moving method, and it was in the computations reported in [10,97].

Element-based mesh relaxation (EBMR)
The EBMR was introduced in 2013 [75]. It restores, without resorting to remeshing, the mesh integrity lost during the mesh motion. The loss of mesh integrity, though not frequent because of the advanced mesh moving methods used with the ST-SUPS and ST-VMS, can happen in FSI computations with a high level of complexity. A recent example of such complexity is FSI computation of clusters of spacecraft parachutes with modified geometric porosity [74][75][76]78,79,84]. When faced with a loss of mesh integrity, it was proposed in [75] to use the EBMR, where the mesh is "relaxed" without altering the mesh at the fluid-structure interface and thus the mesh integrity is somewhat restored. As commented in [75], this is of course a less costly and less disruptive alternative to remeshing.
In the EBMR, the new mesh has the same number of nodes and elements as before, but some of the nodes are moved slightly to improve the quality of some of the elements. The motion is determined by using the large-deformation mechanics equations and an element-based ZSS (EBZSS). This is essentially a shape generated for each element. By design, the undeformed shape, made of "target elements," is the shape we want to obtain from solving the solid mechanics equations. There are several options for constructing the target element shapes, and those options can be found in [75]. The EBMR was successfully used in FSI computation of clusters of spacecraft parachutes with modified geometric porosity (see [75]).

EBZSS
In the EBZSS, we define the ZSS with a set of positions for each element. Positions of nodes (control points) from different elements mapping to the same node in the mesh do not have to be the same. In the reference configuration, all elements are connected by nodes, and we measure the displacement from that connected configuration. This way of formulating the structural mechanics problem was named in [130] "element-based total Lagrangian (EBTL)" method. The EBTL was used in the EBMR [75]. We explain more in Sect. 2.1.

IPBZSS
In the IPBZSS, we extended the way we define EBZSS to integration-point counterpart of that. The ZSS is represented in terms of the metric tensor. Converting the EBZSS representation to IPBZSS representation is straightforward and the conversion will be exact. Converting the IPBZSS representation to EBZSS representation will, in general, not be exact because the IPBZSS has more parameters than the EBZSS. Here we name this way of formulating the structural mechanics problem "integration-point-based total Lagrangian (IPBTL)" method. We explain more in Sect. 2.2.

Stabilization parameters and their mesh sensitivity
The ST-SUPS, ALE-SUPS, RBVMS, ALE-VMS and ST-VMS all have some embedded stabilization parameters that play a significant role (see [1]). These parameters involve a measure of the local length scale (quite often called "element length") and other parameters such as the element Reynolds and Courant numbers. The interface terms in the ST-SI also involve the local length scale, in the direction normal to the interface. There are many ways of defining the stabilization parameters. Some of the newer options for the stabilization parameters used with the SUPS and VMS can be found in [10,82,83,86,87,121,[146][147][148][149]. Some of the earlier stabilization parameters used with the SUPS and VMS were also used in computations with other SUPG-like methods, such as the computations reported in [61,[150][151][152][153][154][155][156][157][158][159][160][161].
Because the stabilization parameters depend on the local length scale, we expect them to have increased mesh sensitivity near the solid surfaces as the mesh moves, because the elements there are typically very thin in the direction normal to the surface. This is another reason to have a good mesh moving method, with the thin layers of elements near solid surfaces shielded from distortion. The SEMMT (see Sect. 1.4.3) accomplishes that to a large extent.

Mesh moving based on fiber-reinforced hyperelasticity and optimized ZSS
In the mesh moving method we introduce in this article, with the fibers placed in multiple directions, we stiffen the element in those directions for the purpose of reducing the distortion during the mesh deformation. We optimize the ZSS by seeking orthogonality of the parametric directions, by mesh relaxation, and by making the ZSS time-dependent as needed.

Test computations
We first conduct 2D test computations with a flow domain containing a beam undergoing three different modes of motion: translation, rotation and bending. Then we do 2D test computations with a rectangular domain undergoing compression and expansion. The 3D tests are with a box-shaped flow domain containing a square block undergoing torsion.

Outline of the remaining sections
In Sect. 2, we provide the large-deformation mechanics framework that we use in the mesh moving method we are introducing. We describe the fiber-reinforced hyperelasticity model in Sect. 3, and introduce the mesh moving method in Sect. 4. The test computations are presented in Sect. 5, and the concluding remarks in Sect. 6. The MJBS, which we also use in the test computations for comparison purposes, is described in Appendix A.

Large-deformation mechanics
We provide, mostly from [134], the large-deformation mechanics framework we will use in the mesh moving method we are introducing here. We note that there can be other ways to reach the same mesh moving functionality.
Let Ω 0 ⊂ R n sd be the material domain in the ZSS, where n sd is the number of space dimensions, and let Γ 0 be its boundary. Let Ω t ⊂ R n sd , t ∈ (0, T ), be the material domain in the deformed state, and let Γ t be its boundary. The structural mechanics equations based on the total Lagrangian formulation can be written as Here, y is the displacement, w is the virtual displacement, δE is the variation of the Green-Lagrange strain tensor, S is the second Piola-Kirchhoff stress tensor, ρ 0 is the mass density in the ZSS, f is the body force per unit mass, and h is the external stress vector applied on the subset (Γ t ) h of the boundary Γ t .

EBZSS
In the EBTL method the ZSS is defined with a set of positions X e 0 for each element e. Positions of nodes from different elements mapping to the same node in the mesh do not have to be the same. In the reference state, X REF ∈ Ω REF , all elements are connected by nodes, and we measure the displacement y from that connected state. The implementation of the method is simple. The deformation gradient tensor F is evaluated for each element: where x is the position in the deformed configuration. The deformation gradient tensors for different elements are on different states, but the terms in Eq. (1), including the second term, do not depend on the orientation. Therefore the rest of the process is the same as it is in the total Lagrangian formulation.

IPBZSS
In the IPBTL method, we extend the way we see X e 0 to integration-point counterpart of X e 0 . As we did with the EBZSS, we work with the reference domain. With the reference Jacobian Eq.
(1) can be rearranged as In our implementation, we use the natural coordinates, with covariant basis vectors where ξ α is the parametric coordinate, and α = 1, . . . , n pd , with n pd being the number of parametric dimensions. The contravariant basis vectors can be calculated with the metric tensor components as where With those vectors, we can express the deformation gradient tensor: and the Cauchy-Green deformation tensor: The Jacobian, J = det F, can be expressed as We can write det C as and from that we obtain We define the covariant basis vectors corresponding to X REF : and the components of the metric tensor are In the process leading to Eq. (20), we replace g α and g β with (G REF ) α and (G REF ) β and obtain an alternative to the expression given by Eq. (4): The Green-Lagrange strain tensor, where I is the identity tensor, can be expressed with the contravariant basis vectors as The second Piola-Kirchhoff tensor can be expressed with the covariant basis vectors as where S αβ can be expressed with the components of the metric tensors. Thus, the inner product δE : S, and all the other quantities in the weak form given by Eq. (5), can be evaluated without actually using the basis vectors G α . Based on that, we rearrange the formulation and obtain the "integration-point-based updated Lagrangian (IPBUL)" method: and We note that we obtained the second term in Eq. (27) from = ε ε ε(w) : g α g β S αβ (34) and dΩ t = J −1 dΩ 0 .

Hyperelasticity model
Here we consider a compressible-material hyperelasticity model, where the volumetric part of the strain energy is represented separately. The general strain-energy-density function can be written as with ϕ iso and ϕ vol representing the isochoric and volumetric parts. The isochoric part can be expressed with the invariants: making ϕ iso = ϕ iso (I 1 , I 2 ). In most models, ϕ iso = ϕ iso (I 1 ).

Remark 1
The last two terms on the right-hand-side of Eqs. (37) and (38) are to make the invariant calculations consistent with the plane-strain assumption in 2D computations.
For the volumetric part, we use the following expression: where κ B is the bulk modulus and β B is a positive number based on [162]. However, we do not see anything wrong with using negative values for β B , and actually we have been using β B = −2 for nearly incompressible materials. We also note that the expression is defined for β B → 0 and gives With β B < 0, the material is stiffer in expansion than in shrinkage, and with β B > 0, it is stiffer in shrinkage than in expansion (see Fig. 1).

Need for fiber-reinforced hyperelasticity
Consider a square domain expanding horizontally in a volume-preserving fashion (see Fig. 2). How an element deforms, and consequently how the mesh lines change, will depend on how the element is oriented. In Fig. 3, the element has an upright orientation and therefore does not undergo distortion. The mesh lines retain their orthogonality. In Fig. 4, the element has a rotated orientation and therefore undergoes distortion. The mesh lines lose their orthogonality. Clearly, the constitutive model does not see the mesh lines, and keeping the mesh-line orthogonality is not an objective of the model.

Constitutive model for fiber-reinforced hyperelasticity
The constitutive model is obtained by adding an invariant that represents the fiber. We consider two options for that invariant, which are commonly seen in the literature. In the first option, where λ a is the stretch in the fiber direction (see [163]). In the second option, where λ a is the stretch associated with C. With the index k being the fiber-direction counter, we write and

Remark 2
We note that there are other invariants [164], such as but we do not consider them here.
With n fib fiber directions, we express the strain-energydensity function with ϕ (C) = ϕ iso I 1 , I 2 + n fib k=1 ϕ ani I 1 , I 2 , I 2k+2 + ϕ vol (J ), (48) or with its slightly modified version: We use the neo-Hookean model for the isotropic part: where μ is the shear modulus. The anisotropic part, from [165], is given as where C 1 is a stress-like coefficient, and C 2 and 0 ≤ κ < 1 3 are nondimensional parameters. When κ = 1 3 , the model becomes an isotropic Fung-type model.

Mesh-moving method
Here we explain the method and options for the mesh moving. We assume that the initial mesh we begin the mesh relaxation with has reasonable resolution and aspect ratio. Most of the time we create thin layers of elements near the solid surfaces and we want the mesh moving method to retain such features.

Mesh relaxation method
The objective of the mesh relaxation is to have a better quality mesh and to have an equilibrium state with the optimized ZSS, boundary conditions and the constitutive law.
We let x 0 denote our initial mesh. We evaluate the corresponding metric tensor at each integration point from the expressions given by Eqs. (6) and (10). When the parametric directions are orthogonal, the metric tensor does not have off-diagonal terms. As the simplest choice, we just remove the off-diagonal terms, and write the ZSS as We write the volume-preserving version of that as We note that √ G αα / G ββ represents the aspect ratio between the α and β directions.
Let a = a α G α denote the fiber direction in the ZSS. We recall that the metric tensor was used in Eqs. (1), (5) and (27) to represent the ZSS, without using the basis vectors G α . Similarly, a α will be used in the constitutive model to represent the fiber direction, without using the basis vectors G α . We will use 2 n pd −1 fiber directions, and we recall that the index k is the fiber-direction counter. We propose two options for defining a α k . In the first option, We omit the normalization a k = 1. In the second option, Beginning with x = x 0 , metric tensor G αβ , reinforcement fibers, and appropriate boundary conditions, including the prescribed boundary positions, we carry out a steady-state structural mechanics computation based on Eq. (27). The constitutive models and parameters can be defined individually for the elements or mesh regions. To avoid any negative impact of the ZSS modifications given by Eqs. (52) and (53), we put the modification into effect gradually: The ramp-up parameter s is increased from 0 to 1 typically in 10 equal increments, with 3-5 nonlinear iterations at each increment to obtain the steady-state solution.
If the mesh quality is not satisfactory, we repeat the process with x 0 ← x or adjust the metric tensor. The key point here is that the mesh x should be in equilibrium state with the ZSS before we start the mesh moving computation.

Mesh moving computation
With the starting mesh x, metric tensor G αβ , reinforcement fibers, and the prescribed boundary positions coming from the time-dependent boundary motion, we solve the structural mechanics equations at each time step to compute the mesh positions. Let b n and b n+1 be the prescribed boundary positions at time levels n and n +1, and let x n and x n+1 be the rest of the mesh positions at those time levels. As we move from time level n to n + 1, sometimes the boundary motion might be too large relative to the thickness of the layers of elements close to the boundary. Simply setting the initial guess in the steady-state computation as where the superscript is the iteration counter, would very likely lead to convergence difficulties. As an option, we have been setting the initial guess as which requires additional positions x n−1 as input, and even with that it may not work when there are sudden condition changes, such as reversing from an expansion to a compression. Here we propose an iteration sequence i = 1, 2, . . . , L,  Table 2 Material properties for tests with the neo-Hookean model ("NH") and fiber-reinforced hyperelasticity ("FR") where the prescribed boundary positions are applied gradually with a ramp-up factor 0 <α 1 ≤ . . . ≤α L−1 ≤α L = 1.
We form the increment equations by evaluating the residuals and derivatives at with b 0 n+1 = b n+1 and x 0 n+1 = x n . With the increment Δx i n+1 , we update b n+1 and x n+1 by extrapolation: We note that for linear mesh moving equations, such as the linear-elasticity equations, the solution is independent ofα i .
Here we also introduce the option of changing the ZSS at each time level. For example, we can calculate the domain volume as it changes in time, and scale with that the metric tensor globally or in selected mesh regions, excluding regions near the solid surfaces. That would give us a good meshresolution control through ϕ vol . To avoid any negative impact on the solution, the ZSS change should be gradual.

Test computations
Here we compare the MJBS (see Appendix A), neo-Hookean model, and fiber-reinforced hyperelasticity. Both hyperelasticity models are based on Eqs. (48), (39), (50) and (51), with κ = 0. For both hyperelasticity models, the ZSS is based on Eq. (53). The fiber directions are based on Eq. (55). Table 1 shows the material properties for the MJBS, and Table 2 shows the material properties for the hyperelasticity models. We are working with nondimensional numbers. We note that for the hyperelasticity models, the Poisson's ratio can be seen as   Figure 5 shows the problem setup. There are three deformation patterns: translation, rotation and bending. We use the material properties (MJBS1, NH1, FR1) for two layers of ele-   Figure 6 shows the initial quadratic NURBS mesh x 0 . The number of elements and control points are 3220 and 3626. For the numerical integration, 4×4 quadrature points are used. We do not use mesh relaxation with the MJBS. Figure 7 shows the mesh after the relaxation. Because the relaxations are based on different material models, the two meshes are slightly different. In both cases the mesh-line orthogonality is improved with the relaxation. Table 3 shows the overall mesh-orthogonality statistics in the whole domain, for the initial mesh and meshes after relaxation.

2D tests with a beam: translation, rotation and bending
The translation, rotation and bending take place in 500 steps. The number of iterations per step (see Sect. 4.2) is 6, with α i = [0.01, 0.04, 0.16, 1, 1, 1]. We note that for the MJBS we have only one iteration per step. Figures 8,9 and 10 show the mesh at the end of the translation, rotation and bending, respectively. For all the deformation patterns, the MJBS and fiber-reinforced hyperelasticity result in a better mesh than the neo-Hookean model does. Comparing the  MJBS and fiber-reinforced hyperelasticity, we see that with the MJBS, distortion is more localized, close to the beam, while with the fiber-reinforced hyperelasticity, it spreads more to regions not close to the beam. Tables 4, 5 and 6 show, at the maximum translation, rotation and bending, the overall mesh-orthogonality statistics in the whole domain.  Overall mesh orthogonality at the maximum bending, with the MJBS, neo-Hookean model ("NH") and fiber-reinforced hyperelasticity ("FR"). L 2 and max norms of the deviation of the mesh-line angles from π/2  Figure 11 shows the rectangular domain. The mesh first undergoes compression, then expansion. We use the material properties (MJBS2, NH3, FR3) for all the elements. Even though the domain is rectangular, we intentionally use a mesh with non-orthogonal mesh lines to see the effect of the relaxation. Figure 12 shows the initial quadratic NURBS mesh x 0 . The number of elements and control points are 2600 and 2928. For the numerical integration, 4×4 quadrature points are used. We do not use mesh relaxation with the MJBS. Fig-Fig. 12 2D tests with a rectangular domain undergoing compression and expansion. The initial mesh x 0  Overall mesh orthogonality for the initial mesh and meshes after relaxation with the neo-Hookean model ("NH") and fiber-reinforced hyperelasticity ("FR"). L 2 and max norms of the deviation of the meshline angles from π/2 ure 13 shows the mesh after the relaxation. Again, with both models, the mesh-line orthogonality is improved. At points where three patches meet, it is impossible to have orthogonality, the angle between the mesh lines becomes close to 2π/3, and that is the best we can expect. Table 7 shows the overall mesh-orthogonality statistics in the whole domain, for the initial mesh and meshes after relaxation. The compression takes place in 250 steps, and the expansion in 500 steps beyond that. The number of iterations per step is 6, with α i = [0.01, 0.04, 0.16, 1, 1, 1]. We note that for the MJBS we have only one iteration per step. Figures 14 and 15 show the mesh at the end of the compression and expansion, respectively. In all three cases the mesh quality is reasonably good. When we look closer, we see that the neo-Hookean model and fiber-reinforced hyperelasticity give smoothness across patch boundaries, which is because of targeting mesh-line orthogonality. Tables 8 and 9 show, at  Figure 16 shows the box-shaped domain containing the block and the block in its initial and twisted states. The block, with Overall mesh orthogonality at the maximum compression, with the MJBS, neo-Hookean model ("NH") and fiber-reinforced hyperelasticity ("FR"). L 2 and max norms of the deviation of the mesh-line angles from π/2 Overall mesh orthogonality at the maximum expansion, with the MJBS, neo-Hookean model ("NH") and fiber-reinforced hyperelasticity ("FR"). L 2 and max norms of the deviation of the mesh-line angles from π/2

3D tests with a square block undergoing torsion
Fig. 16 3D tests with a square block undergoing torsion. Domain containing the block and the block in its initial (red) and twisted (green) states one end fixed centrally on one of the box surfaces, extends halfway through the box. The opposite end of the block undergoes a 90 • rotation about its own center. This motion defines the prescribed boundary positions in computing the rest of the mesh positions in the box. We use the material properties (MJBS1, NH1, FR1) for two layers of elements around the block, and the material properties (MJBS2, NH2, FR2) for the rest of the mesh. Figure 17 shows the initial quadratic NURBS mesh x 0 . The number of elements and control points are 136,000 and 167,728. For the numerical integration, 4×4×4 quadrature points are used. We do not use mesh relaxation with the MJBS. Figure 18 shows the mesh after the relaxation. Table 10 shows the overall mesh-orthogonality statistics in  Overall mesh orthogonality for the initial mesh and meshes after relaxation with the neo-Hookean model ("NH") and fiber-reinforced hyperelasticity ("FR"). L 2 and max norms of the deviation of the meshline angles from π/2 Fig. 19 3D tests with a square block undergoing torsion. The portion of the initial mesh x 0 used in visualizing the deformed states Overall mesh orthogonality at the maximum torsion, with the MJBS, neo-Hookean model ("NH") and fiber-reinforced hyperelasticity ("FR"). L 2 and max norms of the deviation of the mesh-line angles from π/2 the whole domain, for the initial mesh and meshes after relaxation.
The torsion takes place in 500 steps. The number of iterations per step is 3, with α i = [1, 1, 1]. We note that for the MJBS we have only one iteration per step. Figure 19 shows the portion of the initial mesh x 0 used in visualizing the deformed states. Because of the symmetry of the mesh, we will visualize only half of the NURBS patches in that portion. Figure 20 shows the NURBS patches at the end of the torsion. Table 11 shows, at the maximum torsion, the overall mesh-orthogonality statistics in the whole domain. Overall, the behavior tends to be the same as what we observed in the 2D cases.

Concluding remarks
We have introduced a low-distortion mesh moving method based on fiber-reinforced hyperelasticity and optimized ZSS. A good mesh moving method is essential in computation of FSI and MBI problems with moving-mesh methods, which enable mesh-resolution control near the interfaces and consequently high-resolution representation of the boundary layers. The mesh moving method, developed targeting isogeometric discretization, is also applicable to finite element discretization. With the large-deformation mechanics equations, we can expect to have a unique mesh associated with each step of the boundary or interface motion. With the fibers placed in multiple directions, the element is stiffened in those directions for the purpose of reducing the distortion during the mesh deformation. The ZSS is optimized by seeking orthogonality of the parametric directions, by mesh relaxation, and by making the ZSS time-dependent as needed. We conducted both 2D and 3D test computations with isogeometric discretization to demonstrate how the mesh moving method performs. The computations show that the method improves the quality of the initial mesh coming from the mesh generation, giving us a better starting mesh for the mesh motion, and reduces the element distortion experienced during the motion.
where V min and V max are the minimum and maximum element volumes. This can also be written as In practical computations, V max >> V min , and for most of the elements that matter, Therefore Because V max − V min is a global scaling value just like (J M ) 0 is, the stiffening factor given by Eq. (77) is equivalent to the one given by Eq. (73).

Remark 3
We note that a stiffening-factor expression in terms of the mesh Jacobian gives us the option of evaluating the factor at the integration points, without the need for any precalculation. We also note that such an expression is naturally applicable to isogeometric discretization.