A linear-elasticity-based mesh moving method with no cycle-to-cycle accumulated distortion

Abstract

Good mesh moving methods are always part of what makes moving-mesh methods good in computation of flow problems with moving boundaries and interfaces, including fluid–structure interaction. Moving-mesh methods, such as the space–time (ST) and arbitrary Lagrangian–Eulerian (ALE) methods, enable mesh-resolution control near solid surfaces and thus high-resolution representation of the boundary layers. Mesh moving based on linear elasticity and mesh-Jacobian-based stiffening (MJBS) has been in use with the ST and ALE methods since 1992. In the MJBS, the objective is to stiffen the smaller elements, which are typically placed near solid surfaces, more than the larger ones, and this is accomplished by altering the way we account for the Jacobian of the transformation from the element domain to the physical domain. In computing the mesh motion between time levels $t_n$ and $t_{n+1}$ with the linear-elasticity equations, the most common option is to compute the displacement from the configuration at $t_n$. While this option works well for most problems, because the method is path-dependent, it involves cycle-to-cycle accumulated mesh distortion. The back-cycle-based mesh moving (BCBMM) method, introduced recently with two versions, can remedy that. In the BCBMM, there is no cycle-to-cycle accumulated distortion. In this article, for the first time, we present mesh moving test computations with the BCBMM. We also introduce a version we call “half-cycle-based mesh moving” (HCBMM) method, and that is for computations where the boundary or interface motion in the second half of the cycle consists of just reversing the steps in the first half and we want the mesh to behave the same way. We present detailed 2D and 3D test computations with finite element meshes, using as the test case the mesh motion associated with wing pitching. The computations show that all versions of the BCBMM perform well, with no cycle-to-cycle accumulated distortion, and with the HCBMM, as the wing in the second half of the cycle just reverses its motion steps in the first half, the mesh behaves the same way.

1 Introduction

Mesh moving based on linear elasticity and mesh-Jacobian-based stiffening (MJBS) [1,2,3] has been in use with moving-mesh methods since 1992 in computation of flow problems with moving boundaries and interfaces (MBIs), including fluid–structure interaction (FSI). In the MJBS, the objective is to stiffen the smaller elements, which are typically placed near solid surfaces, more than the larger ones, and this is accomplished by altering the way we account for the Jacobian of the transformation from the element domain to the physical domain. In computing the mesh motion between time levels $t_n$ and $t_{n+1}$ with the linear-elasticity equations, the most common option is to compute the displacement from the configuration at $t_n$. While this option works well for most problems, because the method is path-dependent, it involves cycle-to-cycle accumulated mesh distortion. Moving the mesh based on large-deformation mechanics equations, instead of the linear-elasticity equations, would be one good way of addressing the issue. This mesh moving approach has also been in use in FSI and MBI computations, such as those reported in [4,5,6,7,8]. In that direction, a low-distortion mesh moving method based on fiber-reinforced hyperelasticity and optimized zero-stress state (ZSS) was introduced very recently in [9]. As an alternative direction and also very recently, the back-cycle-based mesh moving (BCBMM) method was introduced in [10] as a linear-elasticity-based mesh moving method with no cycle-to-cycle accumulated distortion. In this article, for the first time, we present mesh moving test computations with the BCBMM. We also introduce a version we call “half-cycle-based mesh moving” (HCBMM) method, and that is for computations where the boundary or interface motion in the second half of the cycle consists of just reversing the steps in the first half and we want the mesh to behave the same way. To show how the BCBMM and HCBMM perform, we present detailed 2D and 3D test computations with finite element meshes, using as the test case the mesh motion associated with wing pitching.

The BCBMM and the method introduced in [9] offer alternative solutions to the same issue but have very similar overall objectives. Consequently, the context, history and directly related methods described in the introduction of [9] are also applicable here. Therefore, in much of the rest of this introduction, we include, for completeness, subsections from the introduction of [9].

1.1 Moving-mesh methods

1.1.1 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 high-resolution 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 represented, the boundary layer resolution will be limited by the fluid mechanics mesh resolution where the interface is.

1.1.2 Mixed moving-mesh/nonmoving-mesh methods

As mentioned in [11], 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) [12] was introduced 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 [13] 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.

1.1.3 Space–time (ST) methods

The Deforming-Spatial-Domain/Stabilized ST (DSD/SST) method [14,15,16], 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) [17] and Pressure-Stabilizing/Petrov-Galerkin (PSPG) [14] stabilizations, and therefore it is now called “ST-SUPS.” The ST Variational Multiscale (ST-VMS) method [8, 18, 19], 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 [20,21,22,23]. 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 [18, 19]), are desirable also in computations without MBI.

1.1.4 ALE-SUPS and ALE-VMS methods and classes of problems computed

In the moving-mesh category, the DSD/SST is an alternative to the arbitrary Lagrangian–Eulerian (ALE) method, which is older (see, for example, [24]) and more commonly used. The ALE-VMS method [11, 25,26,27,28,29,30] is the VMS version of the ALE. It succeeded the ST-SUPS and ALE-SUPS [31] and preceded the ST-VMS. To increase their scope and accuracy, the ALE-VMS and RBVMS are often supplemented with special methods, such as those for weakly-enforced Dirichlet boundary conditions [32,33,34] and “sliding interfaces” [35, 36]. The ALE-SUPS, RBVMS and ALE-VMS have been applied to many classes of FSI, MBI and fluid mechanics problems. The classes of problems include ram-air parachute FSI [31], wind turbine aerodynamics and FSI [37,38,39,40,41,42,43,44,45,46,47,48,49], more specifically, vertical-axis wind turbines (VAWTs) [46, 47, 50, 51], floating wind turbines [52], wind turbines in atmospheric boundary layers [45,46,47, 53,54,55], and fatigue damage in wind turbine blades [56], patient-specific cardiovascular fluid mechanics and FSI [25, 57,58,59,60,61,62], biomedical-device FSI [63,64,65,66,67,68,69,70], ship hydrodynamics with free-surface flow and fluid–object interaction [71, 72], hydrodynamics and FSI of a hydraulic arresting gear [73, 74], hydrodynamics of tidal-stream turbines with free-surface flow [75], passive-morphing FSI in turbomachinery [76], bioinspired FSI for marine propulsion [77, 78], bridge aerodynamics and fluid–object interaction [79,80,81], and mixed ALE-VMS/Immersogeometric computations [66,67,68, 82,83,84,85] in the framework of the FSITICT, including heart valve computations. Recent advances in stabilized and multiscale methods may be found for stratified incompressible flows in [86], for divergence-conforming discretizations of incompressible flows in [87], and for compressible flows with emphasis on gas-turbine modeling in [88, 89].

1.1.5 Classes of problems computed with the ST-SUPS and ST-VMS

The ST-SUPS and ST-VMS have also been applied to many classes of FSI, MBI and fluid mechanics problems (see [90] for a comprehensive summary of the work prior to July 2018). The classes of problems include spacecraft parachute analysis for the landing-stage parachutes [4, 11, 91,92,93], cover-separation parachutes [94] and the drogue parachutes [95,96,97], wind turbine aerodynamics for horizontal-axis wind turbine (HAWT) rotors [11, 37, 98, 99], full HAWTs [6, 43, 100, 101] and VAWTs [46,47,48,49, 102, 103], flapping-wing aerodynamics for an actual locust [11, 104,105,106], bioinspired MAVs [6, 101, 107, 108] and wing-clapping [5, 109], blood flow analysis of cerebral aneurysms [101, 110], stent-blocked aneurysms [110,111,112], aortas [7, 10, 69, 70, 113,114,115], heart valves [5, 6, 10, 69, 70, 114, 116,117,118,119,120] and coronary arteries in motion [121], spacecraft aerodynamics [94, 122], thermo-fluid analysis of ground vehicles and their tires [8, 54, 55, 117], thermo-fluid analysis of disk brakes [123], flow-driven string dynamics in turbomachinery [48, 49, 124,125,126], flow analysis of turbocharger turbines [127,128,129,130,131], flow around tires with road contact and deformation [117, 132,133,134,135], fluid films [135, 136], ram-air parachutes [54, 55, 137], and compressible-flow spacecraft parachute aerodynamics [138, 139].

1.2 Special ST methods

A number of special ST methods have been developed in conjunction with the ST-SUPS and ST-VMS to increase their accuracy and scope. They are the ST Slip Interface (ST-SI) method [102, 123], ST Topology Change (ST-TC) method [5, 116], ST Isogeometric Analysis (ST-IGA) [18, 104, 127], and their integration as the ST-SI-TC [132], ST-SI-IGA [127], and ST-SI-TC-IGA [118, 119].

1.2.1 ST-SI

The ST-SI was introduced in [102], in the context of incompressible-flow equations, to retain the desirable moving-mesh 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 [35, 36]. 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 [102]. 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 [123] 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 VAWTs [46,47,48,49, 102, 103], thermo-fluid analysis of disk brakes [123], flow-driven string dynamics in turbomachinery [48, 49, 124,125,126], flow analysis of turbocharger turbines [127,128,129,130,131], flow around tires with road contact and deformation [117, 132,133,134,135], fluid films [135, 136], aerodynamic analysis of ram-air parachutes [54, 55, 137], and flow analysis of heart valves [69, 70, 114, 118,119,120] and ventricle-valve-aorta sequences [10].

In the ST-SI version presented in [102] 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 incompressible-flow aerodynamic analysis of ram-air parachutes with fabric porosity [54, 55, 137]. The compressible-flow ST-SI methods were introduced in [138], 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 [138]. These, together with the compressible-flow ST SUPG method [140], 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 [138, 139].

1.2.2 ST-TC

The ST-TC [5, 116] 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 [5]. 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 [5, 6, 69, 70, 114, 116,117,118,119,120], ventricle-valve-aorta flows [10], wing clapping [5, 109], flow around a rotating tire with road contact and prescribed deformation [117, 132,133,134,135], and fluid films [135, 136].

1.2.3 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 fluid–solid 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 [69, 70, 114, 118,119,120], ventricle-valve-aorta flows [10], tire aerodynamics with road contact and deformation [117, 132,133,134,135], and fluid films [135, 136].

1.2.4 ST-IGA

The success with IGA basis functions in space [25, 35, 57, 141] motivated the integration of the ST methods with isogeometric discretization, which we broadly call “ST-IGA.” The ST-IGA was introduced in [18]. Computations with the ST-VMS and ST-IGA were first reported in [18] 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 [18] 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 [18, 19] and demonstrated in [104, 105, 107]. 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) [104, 105, 107], with the name coined in [100]. 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” [11, 18, 19, 104] 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 [8, 117, 123,124,125,126, 142]. The STNMUM 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 [11, 104,105,106], bioinspired MAVs [6, 101, 107, 108] and wing-clapping [5, 109], separation aerodynamics of spacecraft [94], aerodynamics of horizontal-axis [6, 43, 100, 101] and vertical-axis [46,47,48,49, 102, 103] wind turbines, thermo-fluid analysis of ground vehicles and their tires [8, 54, 117], thermo-fluid analysis of disk brakes [123], flow-driven string dynamics in turbomachinery [48, 49, 124,125,126], flow analysis of turbocharger turbines [127,128,129,130,131], and flow analysis of coronary arteries in motion [121].

The ST-IGA with IGA basis functions in space enables more accurate representation of the geometry and increased accuracy in the flow solution. It accomplishes that with fewer control points, and consequently with larger effective element sizes. That in turn enables using larger time-step sizes while keeping the Courant number at a desirable level for good accuracy. It has been used in ST computational flow analysis of turbocharger turbines [127,128,129,130,131], flow-driven string dynamics in turbomachinery [48, 49, 125, 126], ram-air parachutes [54, 55, 137], spacecraft parachutes [139], aortas [69, 70, 114, 115], heart valves [69, 70, 114, 118,119,120], ventricle-valve-aorta sequences [10], coronary arteries in motion [121], tires with road contact and deformation [133,134,135], fluid films [135, 136], and VAWTs [48, 49, 103]. The image-based arterial geometries used in patient-specific arterial FSI computations do not come from the ZSS of the artery. A number of methods [11, 26, 143,144,145,146,147,148,149,150,151,152] have been proposed for estimating the ZSS required in the computations. Using IGA basis functions in space is now a key part of some of the newest ZSS estimation methods [69, 150,151,152,153] and related shell analysis [154]. The IGA has also been successfully applied to the structural analysis of wind turbine blades [155,156,157,158,159].

1.2.5 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 [127,128,129,130,131], flow-driven string dynamics in turbomachinery [48, 49, 125, 126] and VAWTs [48, 49, 103] 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 spinning-mesh 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 [54, 55, 137] and spacecraft parachute compressible-flow analysis [139] 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 [69, 70, 114, 118,119,120] and ventricle-valve-aorta flow analysis [10] 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 [133,134,135], 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 [135, 136], the ST-SI-TC-IGA enables solution with a computational cost comparable to that of the Reynolds-equation model for the comparable solution quality [136]. With that, narrow gaps associated with the road roughness [135] 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 flow analysis of heart valves [69, 70, 114, 118,119,120], ventricle-valve-aorta flow analysis [10], turbocharger-turbine flow analysis [127,128,129,130,131], spacecraft parachute compressible-flow analysis [139], flow analysis around tires with road contact and deformation [134, 135], and flow analysis of a tsunami-shelter VAWT [49, 103].

1.3 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 moving-mesh 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.

1.4 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.

1.4.1 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.

1.4.2 Mesh moving based on linear elasticity and MJBS

In the mesh moving method introduced in 1992 [1,2,3], 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. When the MJBS was first introduced in [1,2,3], it consisted of simply dropping the 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” in [160]. It was also augmented in [160] to a more extensive kind by introducing a stiffening power $\chi$ that determines the degree by which the smaller elements are rendered stiffer than the larger ones. This approach, when $\chi$ = 1, would be identical to the method first introduced in [1], and most of the time $\chi$ = 1. To also clarify that the “Jacobian” in the name of the method is the mesh Jacobian, in [10] the method was renamed “Mesh-Jacobian-based stiffening.”

1.4.3 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 [1,2,3], 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 [2, 3]). Earlier examples of element layers undergoing rigid-body motion, in combination with deforming structured meshes, were reported in 1990 [14]. 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 [161]. In the SEMMT [11, 161,162,163,164], 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 [161]: 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 [11, 163,164,165], 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 [11, 163,164,165] and a 2D FSI model problem [163, 164].

1.4.4 Path-dependence and cycle-to-cycle accumulated mesh distortion

As mentioned in [9, 10], in computing the mesh motion between time levels $t_n$ and $t_{n+1}$ with the linear-elasticity equations, there are two closely related drawbacks associated with computing the displacement from the configuration at $t_n$. (i) The method is path-dependent. (ii) Once elements accumulate in some region, it is hard to undo that. The path-dependence leads to non-cyclic results in FSI computations that we expect to have cyclic or near-cyclic results and causes cycle-to-cycle accumulated mesh distortion. Because moving the mesh based on large-deformation mechanics equations, instead of the linear-elasticity equations, is an alternative to using the BCBMM, we will include here, from [9], an introductory explanation of that alternative. This introductory material will include the components of the mesh moving method based on fiber-reinforced hyperelasticity and optimized ZSS.

1.4.5 Mesh moving based on large-deformation mechanics

Moving the mesh based on large-deformation mechanics equations was used in a number of FSI and other MBI computations with the ST-SUPS and ST-VMS (see, for example, [4,5,6,7,8]). Additional good features of this method are, as pointed out in [4, 9], 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 [4], and neo-Hookean in [5,6,7,8]. As pointed out in [9], the MJBS can still be part of the mesh moving method, and it was in the computations reported in [7, 8].

1.4.6 Element-based mesh relaxation (EBMR)

The EBMR was introduced in 2013 [4]. 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 [4, 92, 93, 95, 96, 101]. When faced with a loss of mesh integrity, it was proposed in [4] 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 [4], 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 [4]. The EBMR was successfully used in FSI computation of clusters of spacecraft parachutes with modified geometric porosity (see [4]).

1.4.7 Locally-defined ZSS

Locally-defined ZSS originated as arterial ZSS estimation (see [69, 148,149,150,151,152,153]). It was formulated first as the EBZSS in the finite element discretization context [148, 149], then as the EBZSS in the isogeometric discretization context [150, 151], and then as the integration point based ZSS (IPBZSS) in the isogeometric discretization context [152, 153].

In the EBZSS, the ZSS is defined 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 the displacement is measured from that connected configuration. This way of formulating the structural mechanics problem was named “element-based total Lagrangian” (EBTL) method in [148]. The EBTL was used in the EBMR [4].

In the IPBZSS, the way the EBZSS is defined is extended 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. This way of formulating the structural mechanics problem was named “integration-point-based total Lagrangian” (IPBTL) method in [9].

1.4.8 Mesh relaxation and mesh moving based on fiber-reinforced hyperelasticity and optimized ZSS

Mesh relaxation and mesh moving methods based on fiber-reinforced hyperelasticity and optimized ZSS [9] were introduced targeting isogeometric discretization but are also applicable to finite element discretization. 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. The objective of the mesh relaxation is to improve the quality of the mesh after its initial creation and to have an equilibrium state with the optimized ZSS, boundary conditions and the constitutive law. The NURBS mesh used in the computational flow analysis of the tsunami-shelter VAWT [103] was created with the mesh relaxation method.

1.5 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 [11]). 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 [8, 99, 100, 102, 104, 134, 166,167,168,169]. Element lengths, including the direction-dependent element lengths, should have node-numbering invariance for all element types, including simplex elements. The direction-dependent element length expression introduced in [168] meets that requirement. It was proven in [170] that the local-length-scale expression introduced in [169] for isogeometric discretization is invariant with respect to element splitting. The expression was used in [89] in calculating the Courant number associated with the gas turbine computation. 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 [76, 171,172,173,174,175,176,177,178,179,180,181,182,183].

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.

1.6 BCBMM

In the BCBMM, in computing the mesh motion between time levels $t_n$ and $t_{n+1}$ with the linear-elasticity equations, in any cycle, the mesh motion is computed from the configurations in the first cycle. As pointed out in [10], in later cycles, as needed, the configuration the mesh motion is computed from can be changed from the first cycle to a higher cycle. As also pointed out in [10], that would be a good way if there are significant differences between the FSI solutions at different cycles before the solution becomes cyclic or near-cyclic.

1.7 HCBMM

In the HCBMM, for computations where the boundary or interface motion in the second half of the cycle consists of just reversing the steps in the first half and we want the mesh to behave the same way, we modify how we compute the mesh motion in the second half of the first cycle. At any time level in the second half, we compute the mesh motion from the configuration at the corresponding time level in the first half. In subsequent cycles, the mesh motion is computed as in the BCBMM.

1.8 Outline of the remaining sections

The linear-elasticity-based mesh moving method and its versions are described in Sect. 2. The test computations are presented in Sect. 3, and the concluding remark are given in Sect. 4.

2 Mesh moving method

First, for completeness, we include, from [9], the description of the MJBS. Then we describe, partly also from [9], its options related to cycle-to-cycle mesh distortion.

2.1 MJBS

In the MJBS (see Sect. 1.4.2), 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. With the notation of [11], the formulation can be written as

$$\begin{aligned} \int _{\Omega _{\tilde{t}}} \pmb {\varepsilon }(\mathbf{w ^h}): \pmb {\sigma }(\hat{\mathbf{y }}^h-\hat{\mathbf{y }}_{\tilde{t}}^h) \text {d}\Omega&= 0. \end{aligned}$$

(1)

Here $\Omega$ is the domain, with the subscript t denoting its time-dependence, $\pmb {\varepsilon }$ is the strain tensor, $\mathbf{w}$ is the test function, the superscript “h” indicates that the function is coming from a finite-dimensional space, $\pmb {\sigma }$ is the stress tensor, $\hat{\mathbf{y }}^h$ is the displacement from the reference configuration $\mathbf{X} ^h_\text {REF}$ to the current configuration $\mathbf{x} ^h = \mathbf{X} ^h_\text {REF} + \hat{\mathbf{y }}^h$, and $\hat{\mathbf{y }}_{\tilde{t}}^h$ is the displacement from $\mathbf{X} ^h_\text {REF}$ to the time $\tilde{t}$ configuration $\mathbf{x} _{\tilde{t}}^h = \mathbf{X} ^h_\text {REF} + \hat{\mathbf{y }}_{\tilde{t}}^h$ we compute the mesh motion from. The stress tensor is given as

$$\begin{aligned} \pmb {\sigma }(\mathbf{y} )&= \frac{E}{1+\nu } \left( \frac{\nu }{(1 - 2 \nu )} \text {tr}\left( \pmb {\varepsilon }(\mathbf{y} )\right) \mathbf{I} + \pmb {\varepsilon }(\mathbf{y} ) \right) , \end{aligned}$$

(2)

where E and $\nu$ are the Young’s modulus and Poisson’s ratio.

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:

$$\begin{aligned} J_\text {M}&= \det \left( \frac{\partial \mathbf{x} ^h_{\tilde{t}}}{\partial \pmb {\xi }} \right) , \end{aligned}$$

(3)

where $\pmb {\xi }$ is the element coordinate. The subscript “$\text {M}$” is to clarify that this is the mesh Jacobian. 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 [1,2,3], it consisted of simply dropping $J_{\text {M}}$ from Eq. (1). That is equivalent to modifying Eq. (1) as

$$\begin{aligned} \int _{\Omega _{\tilde{t}}} \pmb {\varepsilon }(\mathbf{w ^h}): \pmb {\sigma }(\hat{\mathbf{y }}^h-\hat{\mathbf{y }}_{\tilde{t}}^h) \left( \frac{J_\text {M} }{\left( J_\text {M}\right) _0}\right) ^{-1} \text {d}\Omega&= 0, \end{aligned}$$

(4)

where $\left( J_\text {M}\right) _0$ was originally intended to be a global scaling value and is now seen as a free parameter. The modification results in the smaller elements being stiffened more than the larger ones, and the element stiffening factor is $\frac{\left( J_\text {M}\right) _0}{J_\text {M}}$. The method was augmented in [160] to a more extensive kind by introducing a stiffening power $\chi$ that determines the degree by which the smaller elements are rendered stiffer than the larger ones:

$$\begin{aligned} \int _{\Omega _{\tilde{t}}} \pmb {\varepsilon }(\mathbf{w ^h}): \pmb {\sigma }(\hat{\mathbf{y }}^h-\hat{\mathbf{y }}_{\tilde{t}}^h) \left( \frac{J_\text {M} }{\left( J_\text {M}\right) _0}\right) ^{-\chi } \text {d}\Omega&= 0. \end{aligned}$$

(5)

This approach, when $\chi$ = 1, would be identical to the method first introduced in [1], and most of the time $\chi$ = 1.

Remark 1

A study was conducted in [3] to see the effect of varying $\nu$ on the mesh distortion. For positive values of $\nu$, the best performance was seen in the range $0 \le \nu \le \frac{1}{3}$, with the variations over that range making very little difference.

Remark 2

Studies leading to the selection of $\nu = -0.2$ in the computations reported in [100] indicated that the performance for $\nu < 0$ was comparable to what was seen in the range $0 \le \nu \le \frac{1}{3}$, and sometimes better.

Remark 3

It was shown in [9, 10] that in practical computations the element stiffening method proposed in [166], which is based on element volumes, is approximately equivalent to the MJBS method with $\chi$ = 1.

Remark 4

It was noted in [9, 10] 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. It was also noted in [9, 10] that such an expression is naturally applicable to isogeometric discretization.

2.2 Options related to cycle-to-cycle mesh distortion

The options are defined by how we select $\tilde{t}$ when we are computing the mesh at $t_{n+1}$. The cycle period is represented by the symbol T.

2.2.1 “TZ”

In this option, $\tilde{t} = t_0$, which corresponds to the initial mesh. There is no cycle-to-cycle accumulated distortion, but, depending on the amount of boundary or interface motion, the distortion within a cycle might reach unacceptable levels.

2.2.2 “TN”

In this option, $\tilde{t} = t_n$. This is the most common option. As pointed out in [10], it gives reasonably good performance for most problems and that is partly because when we compute the mesh motion from a configuration close to the current configuration, $J_\text {M}$ provides the method feedback on which elements to protect. While this option works well for most problems, because the method is path-dependent, it involves cycle-to-cycle accumulated distortion.

2.2.3 BCBMM

Two versions of the BCBMM were proposed in [10]. Here we will call those versions “BC1” and “BC2.” In BC1,

$$\begin{aligned} \tilde{t} = t_n - (k-k_\text {FORW}-1)T , \end{aligned}$$

(6)

where $k = \text {floor}(t_n/T) + 1$ is the current cycle number, and the function $\text {floor}(x)$ returns the greatest integer less than or equal to x. We start with $k_\text {FORW}$ = 0. With that, in any cycle, we compute the mesh motion based on the first cycle. In later cycles, as needed, we can set $k_\text {FORW}$ to higher values to compute the mesh motion based on the recent past cycles. As pointed out in [10], that would be a good way if there are significant differences between the FSI solutions at different cycles before the solution becomes cyclic or near-cyclic. In BC2, for $k \ge 2$,

$$\begin{aligned} \tilde{t} = t_{n+1} - (k-k_\text {FORW}-1)T. \end{aligned}$$

(7)

2.2.4 HCBMM

If the boundary or interface motion in the second half of the cycle consists of just reversing the steps in the first half, we expect the mesh to behave the same way. For the first cycle, that would be $\mathbf{x} ^h_t = \mathbf{x} ^h_{T-t}$. With the HCBMM, we can meet that expectation. In the first cycle, we select $\tilde{t}$ as

$$\begin{aligned} \tilde{t} = {\left\{ \begin{array}{ll} t_{n} &{}\text{ if } t_{n} < \frac{T}{2} \\ t_{N-(n+1)} &{}\text{ otherwise } \\ \end{array}\right. }, \end{aligned}$$

(8)

where N is the number of time steps in T. In subsequent cycles, the mesh motion is computed as in the BCBMM.

Fig. 1

Wing pitching. Problem setup

Fig. 2

2D case. Mesh. The green color indicates the inner mesh. (Color figure online)

Fig. 3

2D case. TZ. $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$. Entire mesh (top) and inner mesh (bottom)

Fig. 4

2D case. TN. $\epsilon _{L_2}^h(t, t -(k-2)T)$ for $k \ge 2$. Entire mesh (top) and inner mesh (bottom)

Fig. 5

2D case. TN. $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$. Entire mesh (top) and inner mesh (bottom)

Fig. 6

2D case. BCBMM. $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$. Entire mesh (top) and inner mesh (bottom)

Fig. 7

2D case. HCBMM. $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$. Entire mesh (top) and inner mesh (bottom)

Fig. 8

3D case. Mesh. Root plane (top) and cut mesh (bottom). The green color indicates the inner mesh. (Color figure online)

Fig. 9

3D case. TN. $\epsilon _{L_2}^h(t, t -(k-2)T)$ for $k \ge 2$. Entire mesh (top) and inner mesh (bottom)

Fig. 10

3D case. TN. Mesh displacement magnitude when $\theta = \theta _{\max }$ in the 7th cycle. Root plane (top) and tip plane (bottom)

Fig. 11

3D case. TN. $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$. Entire mesh (top) and inner mesh (bottom)

Fig. 12

3D case. BCBMM. BC2. $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$. Entire mesh (top) and inner mesh (bottom)

Fig. 13

3D case. HCBMM. BC2. $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$. Entire mesh (top) and inner mesh (bottom)

3 Wing pitching

As the mesh moving computation test case, we use the mesh motion associated with wing pitching, with the prescribed angular displacement coming from [184]. Figure 1 shows the problem setup. The wing has an airfoil cross-section, with the chord length varying linearly in the spanwise direction, with a shrinkage factor of 0.6 from the root to the tip. The spanwise length is the same as the root chord length. The pitching motion is a rigid-body rotation around the half-chord point. Working with nondimensional quantities, the angle of attack is given as

$$\begin{aligned} \theta&= \frac{\theta _{\max } + \theta _\text {min}}{2} - \frac{\theta _{\max } - \theta _\text {min}}{2}\cos (2\pi f t), \end{aligned}$$

(9)

where f is the pitching frequency, set as $f = 1$ in the test computations, and $\theta _{\max }$ and $\theta _\text {min}$ are the maximum and minimum values of the angle of attack, set as $\theta _{\max } = 30^{\circ }$ and $\theta _\text {min}= 10^{\circ }$. Only in the first cycle, however, as a special case, $\theta _{\max } = 26^{\circ }$.

On the root plane, the 3D computational domain is the same as the 2D computational domain. We use finite element meshes. The step-time size is 0.05. In the MJBS, we set $E=1$, which has no significance, $\nu =0.3$, and $\chi =1$. In the BCBMM, $k_\text {FORW} = 0$.

3.1 Mesh quality measures

To measure the cycle-to-cycle mesh distortion, we use the $L_2$ norm of the error in the mesh position relative to the mesh position in a reference cycle:

$$\begin{aligned} \epsilon _{L_2}^{h}(t, t_\text {ref}) = \left( \frac{\int _{\Omega _\text {ref}} \left\| \mathbf{x} ^h_{t} - \mathbf{x} ^h_{t_\text {ref}}\right\| ^2 \text {d}\Omega }{\int _{\Omega _\text {ref}}\text {d}\Omega } \right) ^{\frac{1}{2}}, \end{aligned}$$

(10)

where $t_\text {ref}$ is the corresponding time in the reference cycle.

The distortion is measured using the $L_2$ norm of the relative aspect ratio, represented by the symbol $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$ and calculated as

$$\begin{aligned} \Vert \mathbf{f} _\mathrm{AR}\Vert _{2}&= \left( \frac{\sum _{e=1}^{n_{\text {el}}}\left( f_{\mathrm{AR}}^{e}\right) ^2}{n_{\text {el}}}\right) ^{\frac{1}{2}}, \end{aligned}$$

(11)

$$\begin{aligned} f_\mathrm{AR}^{e}&= \left| \ln {\frac{{AR}^{e}}{AR_{0}^{e}}}\right| . \end{aligned}$$

(12)

Here e is the element counter, $n_{\text {el}}$ is the number of the elements, and the subscript 0 indicates the initial mesh.

In 2D, $AR^{e}$ is defined as

$$\begin{aligned} AR^{e} = \frac{\left( l_{\max }^{e}\right) ^2}{A^{e}}, \end{aligned}$$

(13)

where $A^{e}$ and $l^e_{\max }$ are the area and maximum edge length of the element. In 3D, it is defined as

$$\begin{aligned} AR^{e} = \frac{\left( l_{\max }^{e}\right) ^3}{V^{e}}, \end{aligned}$$

(14)

where $V^{e}$ is the volume of the element.

3.2 2D case

Figure 2 shows the mesh. It is made of 12,600 nodes and 6360 quadratic triangular elements. It has an inner mesh of four layers of refined elements.

3.2.1 TZ

When we inspect $\epsilon _{L_2}^h(t, t -(k-2)T)$ for $k \ge 2$, we see that, as expected, it is zero. There is no cycle-to-cycle accumulated mesh distortion. Figure 3 shows $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$, and we can again see that the mesh motion is cyclic. The distortion within each cycle, however, is reaching higher levels than, as we will see in Sect. 3.2.2, what we get with TN in the first four cycles.

3.2.2 TN

Figure 4 shows $\epsilon _{L_2}^h(t, t -(k-2)T)$ for $k \ge 2$. We can clearly see the cycle-to-cycle accumulated mesh distortion. Figure 5 shows $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$, and again we can see that the mesh motion is not cyclic. We can also see that, in the first four cycles, the distortion levels reached within each cycle is lower than what we get with TZ, with the difference being more pronounced for the inner mesh.

3.2.3 BCBMM

For both versions, BC1 and BC2, when we inspect $\epsilon _{L_2}^h(t, t -(k-2)T)$ for $k \ge 2$, we see that, as expected, it is zero. There is no cycle-to-cycle accumulated mesh distortion. Figure 6 shows $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$ for the two versions, and we can again see that the mesh motion is cyclic. We can also see that the distortion levels reached within each cycle is lower than what we get with TZ or TN, significantly lower for the inner mesh.

3.2.4 HCBMM

The airfoil motion in the second half of the cycle consists of just reversing the steps in the first half. The mesh behaves the same way as we can see in Fig. 7, which shows $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$ for both BC1 and BC2, with a barely detectable departure from the expected behavior for BC1.

3.3 3D case

Figure 8 shows the mesh. It is made of 261,198 nodes and 181,139 quadratic tetrahedral elements. As in 2D, it has an inner mesh of four layers of refined elements. The mesh on the root plane is almost the same as the 2D mesh.

3.3.1 TN

Figure 9 shows $\epsilon _{L_2}^h(t, t -(k-2)T)$ for $k \ge 2$. The behavior, as expected, is similar to what we had in the 2D case. Figure 10 shows the displacement magnitude when $\theta = \theta _{\max }$ in the 7th cycle. Figure 11 shows $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$, and again the behavior is similar to what we had in the 2D case.

3.3.2 BCBMM

We only compute with BC2, since that is our preferred version. When we inspect $\epsilon _{L_2}^h(t, t -(k-2)T)$ for $k \ge 2$, we see that, as expected, it is zero. There is no cycle-to-cycle accumulated mesh distortion. Figure 12 shows $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$, and we can again see that the mesh motion is cyclic.

3.3.3 HCBMM

Again, we only compute with BC2. As in 2D, the wing motion in the second half of the cycle consists of just reversing the steps in the first half. The mesh behaves the same way as we can see in Fig. 13, which shows $\Vert \mathbf{f} _\mathrm{AR}\Vert _{2}$.

4 Concluding remarks

Mesh moving based on linear elasticity and MJBS has been in wide use with moving-mesh methods in computation of flow problems with FSI and MBI. In the MJBS, the smaller elements, which are typically placed near solid surfaces, are protected more than the larger ones. In computing the mesh motion between time levels $t_n$ and $t_{n+1}$ with the linear-elasticity equations, computing the displacement from the configuration at $t_n$, which is a very common practice, leads to a path-dependent mesh motion and results in cycle-to-cycle accumulated mesh distortion. The BCBMM was introduced to remedy that and has no cycle-to-cycle accumulated distortion. In the BCBMM, in any cycle, the mesh motion is computed from the configurations in the first cycle. We have, for the first time, carried out mesh moving test computations with the BCBMM. We have also introduced the HCBMM, which is a special-case version of the BCBMM. The HCBMM is for computations where the boundary or interface motion in the second half of the cycle consists of just reversing the steps in the first half and we want the mesh to behave the same way. In the HCBMM, at any time level in the second half of the first cycle, we compute the mesh motion from the configuration at the corresponding time level in the first half. In subsequent cycles, the mesh motion is computed as in the BCBMM. We have presented detailed 2D and 3D test computations with finite element meshes, using as the test case the mesh motion associated with wing pitching. The computations show that all versions of the BCBMM perform well, with no cycle-to-cycle accumulated distortion, and with the HCBMM, as the wing in the second half of the cycle just reverses its motion steps in the first half, the mesh behaves the same way.