Computational engineering analysis with the new-generation space–time methods

This is an overview of the new directions we have taken the space–time (ST) methods in bringing solution and analysis to different classes of computationally challenging engineering problems. The classes of problems we have focused on include bio-inspired flapping-wing aerodynamics, wind-turbine aerodynamics, and cardiovascular fluid mechanics. The new directions for the ST methods include the variational multiscale version of the Deforming-Spatial-Domain/Stabilized ST method, using NURBS basis functions in temporal representation of the unknown variables and motion of the solid surfaces and fluid meshes, ST techniques with continuous representation in time, and ST interface-tracking with topology change. We describe the new directions and present examples of the challenging problems solved.


Introduction
In computational engineering analysis, one frequently faces the challenges involved in solving flow problems with moving boundaries and interfaces (MBI).This wide class of problems include fluid-structure interaction (FSI), fluid-object interaction (FOI), fluid-particle interaction (FPI), free-surface and multi-fluid flows, and flows with solid surfaces in fast, linear or rotational relative motion.The computational challenges still being addressed include accurately representing the small-scale flow patterns, which requires a reliable multiscale method.They also include contact or near contact between moving solid surfaces and other cases of topology change (TC) or near TC, such as those in flappingwing aerodynamics, wind-turbine aerodynamics, and cardiovascular fluid mechanics.These three specific classes of problems played a key role in motivating the development of the computational-analysis methods discussed in this article.
A method for flows with MBI can be viewed as an interface-tracking (moving-mesh) technique or an interfacecapturing (nonmoving-mesh) technique, or possibly a combination of the two.In interface-tracking methods, as the interface moves, the mesh moves to follow (i.e."track") the interface.The Arbitrary Lagrangian-Eulerian (ALE) finite element formulation [1] is the most widely used moving-mesh technique.It has been used for many flow problems with MBI, including FSI (see, for example, ).The Deforming-Spatial-Domain/Stabilized Space-Time (DSD/SST) method [30,[36][37][38][39][40][41][42], introduced in 1992, is also a moving-mesh method.The costs and benefits of moving the fluid mechanics mesh to track a fluid-solid interface were articulated in many papers (see, for example, [39,43]).
Moving-mesh methods require mesh update methods.Mesh update typically consists of moving the mesh for as long as possible and remeshing as needed.With the key objectives being to maintain the element quality near solid surfaces and to minimize frequency of remeshing, a number of advanced mesh update methods [40,[44][45][46][47][48] were developed in conjunction with the DSD/SST method, including those that minimize the deformation of the layers of small elements placed near solid surfaces.
Perceived challenges in mesh update are quite often given as reasons for avoiding interface-tracking methods in classes of problem that can be, and actually have been, solved with such methods.A robust moving-mesh method with effective mesh update can handle FSI or other MBI problems even when the solid surfaces undergo large displacements (see, for example, FPI [46,49] with the number of particles reaching 1,000 [46], parachute FSI [40,43,[50][51][52][53][54][55][56], flapping-wing aerodynamics [57][58][59], and wind-turbine rotor and tower aerodynamics [60].It can handle FSI or other MBI problems also even when the solid surfaces are in near contact or create near TC, if the "nearness" is sufficiently "near" for the purpose of solving the problem.Examples of such problems are FPI with collision between the particles [46,49], parachute-cluster FSI with contact between the parachutes of the cluster [52][53][54]56], flapping-wing aerodynamics with the forewings and hindwings crossing each other very close [57][58][59], and wind-turbine rotor and tower aerodynamics with the blades passing the tower close [60]. No method is a panacea for all classes of MBI problems.As mentioned in [30], certain classes of interfaces, such as those in splashing, might be too complex to deal with an interface-tracking technique, requiring an interfacecapturing technique.The Mixed Interface-Tracking/Interface-Capturing Technique (MITICT) [46] was introduced for computations that involve both fluid-solid interfaces that can be accurately tracked with a moving-mesh method and fluidfluid interfaces that are too complex or unsteady to be tracked.Those fluid-fluid interfaces are captured over the mesh tracking the fluid-solid interfaces.The MITICT was successfully tested in 2D computations with solid circles and free surfaces [61,62] and in 3D computation of ship hydrodynamics [21].
In some MBI problems with contact between the solid surfaces, the "nearness" that can be modeled with a movingmesh method without actually bringing the surfaces into contact might not be "near" enough for the purpose of solving the problem.Cardiovascular FSI with heart valves, where the flow has to be completely blocked at contact, is an example.The Fluid-Solid Interface-Tracking/Interface-Capturing Technique (FSITICT) [51] was motivated by such FSI problems.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.A specific application of the FSITICT was presented in [63], where the ALE method is used for interface tracking, and a fully Eulerian approach for interface capturing, with some 2D benchmark problems as test computations.This specific application was extended in [63] to 2D FSI models with flapping and contact, where the fully Eulerian interface-capturing is complemented with mesh adaptivity.
The original DSD/SST method is based on the SUPG/ PSPG stabilization, where "SUPG" and "PSPG" stand for the Streamline-Upwind/Petrov-Galerkin [139] and Pressure-Stabilizing/Petrov-Galerkin [36,140] methods.Starting in its very early years, the DSD/SST method also included the "LSIC" (least-squares on incompressibility constraint) stabilization.The ST-VMS method [30,41,42] is the variational multiscale version of the DSD/SST method.It was called "DSD/SST-VMST" (i.e. the version with the VMS turbulence model) when it was first introduced in [41].The VMS components are from the residual-based VMS method given in [141][142][143][144].We demonstrated the increased accuracy brought by the ST-VMS method the first time with the computations reported in [41,83,145].We have been using the ST-VMS method in most of our computations since then.The original DSD/SST method was named "DSD/SST-SUPS" in [41] (i.e. the version with the SUPG/PSPG stabilization), which was also called "ST-SUPS" in [30].
The ST methods give us the the option of using higherorder basis functions in time, including the NURBS, which have been used very effectively as spatial basis functions (see [4,8,146,147]).We started using that option with the methods and concepts introduced in [41].This of course increases the order of accuracy in the computations [41,42,129], and the desired accuracy can be attained with larger time steps, but there are positive consequences beyond that.The ST context provides us better accuracy and efficiency in temporal representation of the motion and deformation of the moving interfaces and volume meshes, and better efficiency in remeshing.This has been demonstrated in a number of 3D computations, specifically, flapping-wing aerodynamics [57][58][59]79,80], separation aerodynamics of spacecraft [55], and wind-turbine aerodynamics [60].The mesh update method based on using NURBS basis functions in mesh motion and remeshing was named "ST/NURBS Mesh Update Method (STNMUM)" in [60].
There are some advantages in using a discontinuous temporal representation in ST computations.For a given order of temporal representation, we can reach a higher order accuracy than one would reach with a continuous representation of the same order.When we need to change the spatial discretization (i.e.remesh) between two ST slabs, the temporal discontinuity between the slabs provides a natural framework for that change.There are advantages also in continuous temporal representation.We obtain a smooth solution, NURBS-based when needed.We also can deal with the computed data in a more efficient way, because we can represent the data with fewer temporal control points, and that reduces the computer storage cost.These advantages motivated the development of the ST computation techniques with continuous temporal representation (ST-C) [148].
There are different types of nonmoving-mesh methods that can compute MBI problems involving an actual contact between solid surfaces or other cases of TC.Some of those methods give up on the accurate representation of the interface, and most give up on the consistent representation of the interface motion.The DSD/SST formulation does not need to give up on either, even where we have an actual contact or some other TC, provided that we can update the mesh even there.Using an ST mesh that is unstructured both in space and time, as proposed for contact problems in [46], would give us such a mesh update option.However, that would require a fully unstructured 4D mesh generation, and that is not easy in computing real-world problems.An ST interface-tracking method that can deal with TC was introduced in [138], and it is called ST-TC.It is a practical alternative to using unstructured ST meshes, but without giving up on the accurate representation of the interface or the consistent representation of the interface motion, even where there is an actual contact between solid surfaces or other TC.The ST-TC method is based on special mesh generation and update, and a master-slave system that maintains the connectivity of the "parent" mesh when there is a TC.
In this article we provide an overview of these four new directions we have taken the ST methods and show how that brought solution to the three specific classes of problems mentioned in the first paragraph.In Sect.2, we briefly review the Navier-Stokes equations of incompressible flows.The ST-SUPS and ST-VMS methods are described in Sects.3 and 4. Methods based on temporal representation with NURBS basis functions, including the STNMUM, are given in Sect. 5.The ST-C and ST-TC are described in Sects.6 and 7.In the three sections following that, we present examples of the challenging problems solved.In Sect.8, we present aerodynamic analysis of flapping wings of an actual locust and an MAV, in Sect.9 aerodynamic analysis of wind turbines, and in Sect. 10 a proof-of-concept computation with two pairs of symmetrically-flapping surfaces with coordinated opening/closing actions.The concluding remarks are given in Sect.11.

Governing equations
Let Ω t ⊂ R n sd be the spatial domain with boundary Γ t at time t ∈ (0, T ).The subscript t indicates the time-dependence of the domain.The Navier-Stokes equations of incompressible flows are written on Ω t and ∀t ∈ (0, T ) as where ρ, u and f are the density, velocity and the external force, respectively.The stress tensor σ σ σ is defined as Here p is the pressure, I is the identity tensor, μ = ρν is the viscosity, ν is the kinematic viscosity, and ε ε ε(u) is the strainrate tensor.The essential and natural boundary conditions for Eq. ( 1) are represented as u = g on (Γ t ) g and n • σ σ σ = h on (Γ t ) h , where (Γ t ) g and (Γ t ) h are complementary subsets of the boundary Γ t , n is the unit normal vector, and g and h are given functions.A divergence-free velocity field u 0 (x) is specified as the initial condition.

ST-SUPS (DSD/SST-SUPS) method
In the DSD/SST method (see, e.g., [36][37][38][39][40][41][42]), the finite element formulation is written over a sequence of N ST slabs Q n , where Q n is the slice of the ST domain between the time levels t n and t n+1 (see Fig. 1).At each time step, the integrations are performed over Q n .The ST finite element interpolation functions are continuous within a ST slab, but discontinuous from one ST slab to another.The notation (•) − n and (•) + n will denote the function values at t n as approached from below and above.Each Q n is decomposed into elements Q e n , where e = 1, 2, . . ., (n el ) n .The subscript n used with n el is for the general case where the number of ST elements may change from one ST slab to another.The essential and natural boundary conditions are enforced over (P n ) g and (P n ) h , the complementary subsets of the lateral boundary of Fig. 1 ST slab in an abstract representation the ST slab.The finite element trial function spaces S h u for velocity and S h p for pressure, and the test function spaces V h u and V h p = S h p are defined by using, over Q n , typically first-order polynomials in space and time, but can also be of higher-order functions.
The ST-SUPS (DSD/SST-SUPS) method (from [39]) is written as follows: given where r M and r C are the the residuals of the momentum equation and incompressibility constraint (continuity equation).The ST-SUPS method has all the ingredients of the semi-discrete SUPG/PSPG finite element formulation.That includes the test functions, domain integrations, stress terms that have been integrated by parts, boundary integrations, and the SUPG, PSPG and LSIC stabilization terms with stabilization parameters τ SUPG , τ PSPG and ν LSIC .The stabilization is residual based because r M and r C appear as factors in the stabilization terms.The stabilization parameters τ SUPG , τ PSPG and ν LSIC originate from stabilized finite element methods for fluid dynamics (see, e.g., [36,39,139,[149][150][151][152]).There are various ways of defining these parameters.
In the computations included in this article, we mostly use the definitions given in [39], with some new options introduced in [40,60,79,145].For more ways of calculating the stabilization parameters in finite element computation of flow problems, see [30,106,.

Conservative form
The conservative form of the ST-VMS method is written as follows: given The stabilization parameter τ SUPS is calculated in essentially the same way as τ SUPG is calculated.The notation "SUPS," introduced in [41], indicates that there is a single stabilization parameter for the SUPG and PSPG stabilizations, instead of two separate parameters.
Remark 1 One of the main differences between the ALE and ST forms of the VMS method is that the ST form retains the fine-scale time derivative term ∂u ∂t ξ ξ ξ (ξ ξ ξ is the vector of element coordinates).Dropping this term is called the "quasistatic" assumption (see [15] for the terminology).This is the same as the "WTSE" option in the DSD/SST formulation (see Remark 2 of [40]).We believe that this makes a significant difference, especially when the polynomial orders in space or time are higher (see [41]).

Convective form
Remark 2 The 6th and 7th terms of the ST-VMS method in either form are the SUPG/PSPG and LSIC stabilization terms, respectively.If we exclude the last two terms of the convective form, the method reduces to the ST-SUPS (DSD/SST-SUPS) method under the condition τ PSPG = τ SUPG .

Temporal representation with NURBS basis functions
The concept of using NURBS basis functions, in conjunction with the ST methods, in temporal representation of the unknown variables and motion of the solid surfaces and fluid meshes was first introduced in [41].

ST basis functions
An ST basis function can be written as a product of its spatial and temporal parts: where θ ∈ [−1, 1] is the temporal element coordinate, and n en and n ent are the number of spatial and temporal element nodes.Figure 2 shows an example of temporal NURBS basis functions.Temporal NURBS basis functions can be used in an ST slab for the representation of the unknown variables and test functions as well as the spatial coordinates.
As pointed out in [30,41,42,60,79], different components (i.e.unknowns), and the corresponding test functions, can be discretized with different sets of temporal basis func-Fig.2 Temporal NURBS basis functions tions.This was shown in [30,41,42,60,79] by introducing a secondary mapping Θ ζ (θ ) ∈ [−1, 1], where Θ ζ (θ ) is a strictly increasing function, and rewriting the generalized ST basis function for the element indices (a, α) as For example, we can discretize time and position as Here Θ t (θ ) and Θ x (θ ) are the secondary mappings for time and position, and t α and x α are the time and position values corresponding to the basis function T α .

Motion of solid surfaces
As an example, Fig. 3 shows, from [79], how the motion of the forewing (FW) of a locust is represented temporally.
In the preliminary computations reported in [79], quadratic NURBS basis functions were used in the temporal representation of the wing motion.Based on those computations, using even higher-order temporal basis functions was proposed in [79], so that the acceleration is continuous.Fig. 4 shows, from [57], how the path of a point on the hindwing (HW) of the locust is represented with cubic NURBS basis functions over four temporal patches (see [57] for details).

Rotation representation with constant angular velocity
With temporal NURBS basis functions, as described in [30,41,42,60,79], we can represent a path accurately, such as representing a circular arc exactly with quadratic NURBS.We can also specify a speed along that path, such as a constant angular velocity for a circular arc [30,60,79].In this this subsection, from [60], we describe how to do that.
For the circular arc, with quadratic NURBS, n ent = 3.With the secondary-mapping concept described in Sect.5.1, the velocity can be expressed as follows: Thus, the speed along the path can be specified only by modifying the secondary mapping.For a circular arc, two methods were introduced in [79] and also described in [30]; one is modifying the secondary mapping for position and the other one is modifying both such that dt dθ is constant.We note that, in theory, the secondary-mapping selections do not make any difference as long as the relationship dΘ x dΘ t is the same.In our implementation, to keep the process general, we search for the parametric coordinate θ by using an iterative solution method [30,79].We use the latter set of the secondary mappings, having constant dt dθ .

Mesh computation and representation
Given the fluid mechanics mesh on a moving solid surface, we compute the fluid mechanics volume mesh using the mesh moving techniques [40,[44][45][46] developed in conjunction with the DSD/SST method.As proposed in [79] and also described in [30], these mesh moving techniques are applied to computing the meshes that will serve as temporalcontrol points.This allows us to do mesh computations with longer time in between, but get the mesh-related information, such as the coordinates and their time derivatives, from the temporal representation whenever we need.Obviously this also reduces the storage amount and access associated with the meshes.However, because of the longer time between the control meshes, linear interpolation of the surfaces between control points in time might be needed in computing those meshes with the mesh moving technique mentioned.
Remark 3 Getting the meshes used in the computations from the temporal representation can be done independent of which time direction was used in computing the control meshes.For example, in flapping-wing aerodynamics, it does not matter if the control meshes were computed while the wings were flapping forward or backward in time.

Remeshing
In many computations remeshing becomes unavoidable.Two choices were proposed in [79] and also described in [30].To explain those choices, let us assume that when we try to move from control mesh , we find the quality of M β+1 c to be less than desirable.In the first choice, called "trimming," we remesh going back to M , where p is the order of the NURBS basis functions.Then whenever our solution process needs a mesh, depending on the time, we use the control meshes belonging to either only the un-remeshed set or only the remeshed set (Fig. 5).In the second choice, we perform knot insertion p times in the temporal representation of the surface at the right-most knot before the maximum value of the basis function corresponding to t β+1 c , making that knot a new patch boundary.Then we do the mesh moving computation for the control meshes associated with the newly-defined basis functions, not only the one at the new patch boundary, but also going back ( p − 1) basis functions (Fig. 6).We favor the second choice, because we believe that in many cases the need for remeshing is generated by a topological change, which we can avoid going over with a large step if we use the knot insertion process.

Fig. 5 Remeshing and trimming NURBS. A set of un-remeshed meshes (top). A set of remeshed meshes (middle). Common basis functions (bottom)
Fig. 6 Remeshing with knot insertion.For the set of un-remeshed meshes, there are p newly-defined basis functions and the corresponding control points are marked "New."We carry out the mesh moving computations for those meshes

ST-C
In this section, from [148], we describe the version of ST-C used in extracting continuous temporal representation from Fig. 7 Continuous solution (top) and its basis functions (bottom), where ϑ is the parametric coordinate computed data.This is essentially a post-processing method, and can also be seen as a data compression method.For the version used in direct computation of the solution with continuous temporal representation, see [148].

Least-squares projection for full temporal domain
When we have the complete sequence of computed data, we can project that to a smooth representation, with basis functions that provide us that smooth representation, such as NURBS basis functions.As an example, Fig. 7 shows the goal continuous data φ C and its basis functions, where ϑ denotes the parametric temporal coordinate.The projection for each spatial node can be done independently from the other nodes.Consider the time-dependent, typically discontinuous computed data φ D for a node.We define the basis functions as T α C , where α = 0, 1, . .., and the coefficients to be determined in the projection as φ α .We use a standard least-squares projection: given φ D , find the solution φ C ∈ S C , such that for all test functions w C ∈ V C : where T represents time period of the computation, and S C and V C are the solution and test function spaces constructed from the basis functions.This approach requires that we store all the computed data before the projection, and that would be a significant computer storage cost when the number of time steps is large.

Successive-projection technique
In ST-C with the successive-projection technique (ST-C-SPT), we extract the continuous solution shown in Fig. 7 without storing all the computed data.We describe the technique here for the special case with quadratic B-splines.To explain the successive nature of the SPT, let us suppose that we have the continuous solution extracted up to t n = 4.0, as shown in Fig. 8.We assume that this continuous solution, which we will call φ C , has already replaced φ D up to t n = 4.0.With that, we describe how we extract the continuous solution up to t n+1 = 5.0, as shown in Fig. 9.With the newly computed data φ D between t n = 4.0 and t n+1 = 5.0, we solve the following projection equation: given φ D on t ∈ (4.0, 5.0), φ C on t ∈ [2.0, 4.0], and We note that Eq. ( 13) is essentially used for defining the coefficients φ α C , α = 4, 5, 6, which correspond to the basis functions T α C .How to deal with the initial part of the extraction, description of the ST-C-SPT for the general case (i.e.beyond quadratic B-splines), and comments on efficient implementation of the SPT can be found in [148].We consider two hypothetical cases of two bars to provide a context for TC.In the first case, shown in Fig. 10, the bars are initially coinciding, with just one hole in the fluid mechanics domain.Then the red bar starts moving upward, creating a second hole.In the second case, shown in Fig. 11, the bars are initially aligned with connected ends, again with a single hole in the domain.Then the red bar starts a flapping motion, up and down, creating a second hole in the domain, except when their ends become connected periodically during the flapping motion.When the red bar is in the upper position, the part of the domain below it is connected to the part of the domain above the blue bar.When the red bar is in the lower position, the part of the domain above it is connected to the part of the domain below the blue bar.These two cases are representatives of the typical TC challenges we expect to see in the classes of MBI problems we are targeting.Especially the first case is really not possible to treat in a consistent way without using an ST method.

Master-slave system
We propose a very simple technique in the ST context.Having a constraint between nodes in a finite element formulation is quite common.These constraints reduce the number of unknowns, but in our implementation we delay that unknown elimination until the iterative solution of the linear systems encountered at nonlinear iterations of a time step.The iterative solution of the linear systems is performed with reduced number of unknowns.The technique is easy to manage in a parallel-computing environment, especially if the precondi-tioner is simple enough.Typically we assign a master node to each slave node, and we use only the unknowns of the master nodes in iterative solution of the linear systems.We can use different master-slave relationships at different time levels.This is a practical alternative to, but less general than, using ST meshes that are unstructured in time.Still, we can use this concept to deal with the TC cases considered above, and the important point is that the connectivity of the "parent" mesh does not change.Consequently, the distribution model in the parallel-computing environment does not change during the computations with moving meshes.
With this technique, we need to implement one more functionality.We exclude certain elements from the integration of the finite element formulation.The exclusion principles are given below.
-Exclude all spatial elements with zero volume from the spatial integration.-Exclude all ST elements with zero ST volume from the ST integration.-We assume that checking if an ST element has zero ST volume is equivalent to checking if all the spatial elements associated with that ST element have zero volume.Therefore, for this purpose, we check the spatial-element volumes.-To identify the spatial elements with zero volume, which should have zero Jacobian at all the integration points, instead of evaluating the Jacobians, we make the determination for a given spatial element from the master-slave relationship of its nodes.The method is explained more in [138].

Design of the master-slave system
The data we need to provide to the solver is simple.It is just the master-slave relationship at each time level.However there are some restrictions, and here we explain the three that we want to emphasize.The first restriction is that we cannot have a node which is not part of any active (nonzero volume) spatial element.This is because the values at such nodes would no longer be in our equation system, and therefore would become undefined.If because of another TC such a node comes back to the equation system later as part of an active element, it would add an undefined component to the equation system.The second restriction is that when we construct the ST elements, we have to have matching lateral element-boundary faces between the active adjacent ST elements.This condition cannot be checked on the spatial mesh.Therefore we need to check it on the ST mesh.The third one is related to implementation.The master-slave relationship also extends to cases when we have boundary conditions on the master and slave nodes.In other words, the conditions at the master node also apply to the slave nodes.

Contraction and expansion
This is related to the first one of the two cases of TC described in Sect.7.1.Contraction and expansion are basically the same, except having different directions in time progression.Figure 12 shows a contraction example.The spatial element with nodes 1 and 2, for example, has zero volume at the first time level.However, it has nonzero volume at the second time level, and therefore the corresponding ST element has nonzero volume.

Flapping
This is related to the second one of the two cases of TC described in Sect.7.1.Figure 13 shows the red and blue bars at three instants in time as the red bar crosses the blue bar. Figure 14 shows, for the flapping motion, the ST trajectories of the neighboring ends of the blue and red bars.Figure 15 shows the ST element edges for the line separating the two sides of the domain containing the blue and red bars (shown as the vertical dashed line in Fig. 13).For each side of the domain, the spatial node motions along the ST element edges  have to be designed in a fashion that does not lead to mesh entanglement.Figure 16 shows the master-slave relationship for the blue-bar side of the domain, and Fig. 17 the redbar side.In addition, those two sides are in a master-slave relationship along the vertical dashed line in Fig. 13.

Aerodynamic analysis of flapping wings of an actual locust and an MAV
This section is from [57,58].The fluid mechanics computations are carried out with the conservative form of the ST-VMS method.More information on the computational para-  meters and conditions can be found in [57,58].The motion and deformation data for the wings is extracted from the high-speed, multi-camera video recordings of a locust in a wind tunnel at Baylor College of Medicine (BCM), Houston, Texas.The video recording is accomplished by using a set of tracking points marked on the FWs and HWs of the locust.
The tracking points can be seen in Fig. 18.The interested reader can find in [57] the details of how the wing motion and deformation data is extracted from the video data and represented in space and time with the methods described in Sect.5, including the STNMUM, and some additional techniques.Figure 19 illustrates, in the context of the HW wingtip trajectory, how the STNMUM is used in the computation.Figure 20 show how the body and wings compare for the locust and MAV models.Figure 21 shows for the locust the vorticity magnitude during the second flapping cycle.
Figure 22 shows for the MAV the vorticity magnitude during the third flapping cycle.In Figs.21 and 22, the color range from blue to red corresponds to a vorticity range from low to high, and lighter and darker shades of a color correspond to lower and higher values.Figure 23 shows the lift and thrust generated by the locust and MAV.

Aerodynamic analysis of wind turbines
This section is from [83] and [60].Computer modeling of wind-turbine aerodynamics is challenging because correct aerodynamic torque calculation requires correct separationpoint calculation, which requires an accurate flow field, which in turn requires good mesh resolution and turbulence model.We use an actual wind-turbine model, which is NREL 5MW offshore baseline wind turbine, and the geometric com- plexity is also a computational challenge.Including the tower in the model increases the computational challenge because of the mesh update requirements of the fast, rotational relative motion between the rotor and tower.
First we describe how we computed the aerodynamics of a rotor without the tower by using the ST-SUPS and ST-VMS methods.More information on the computational parameters and conditions can be found in [83].Figure 24 shows, from [15], the airfoil cross-sections of the wind-turbine blade superposed on the blade.Figure 25 shows time history of the aerodynamic torque generated by a single blade, as computed with the ST-SUPS, ST-VMS (conservative form) and ALE methods.The ALE results are from [15], which we take as the reference solution.The figure demonstrates how the ST-VMS method increases the accuracy in this particular computational analysis.
When we include the tower, we deal with the mesh update requirements with the methods described in Sect.5, including the STNMUM.Figure 26 illustrates, in the context of the blade tip trajectory, how the STNMUM is used in the computation.More information can be found in [60].Figure 27 shows the vorticity magnitude, computed with the conservative form of the ST-VMS method and the STNMUM.In that figure, the color range from blue to red corresponds to a vorticity range from low to high, and lighter and darker shades of a color correspond to lower and higher values.Figure 28 shows the torque for the individual blades.

Remark 4
The drop in the aerodynamic torque as the blade passes the tower is a well-known phenomenon, observed in experiments and in other computations of wind-turbine aerodynamics (see, for example, [35]).In [35], a slidinginterface technique [147] was used in conjunction with the ALE-VMS formulation [20] to compute wind-turbine aerodynamics, including the rotor-tower interaction.There are two pairs of symmetrically-flapping surfaces with zero thickness, positioned and functioning like the mitral and aortic valves would be.The 2D domain changes its area like the left ventricle would change its volume.When one of the pairs closes, the domain is separated into two.Figure 29 shows the computational domain.The flow enters through the flapping pair at the inlet (we call this pair "mitral"), and exits through the pair at the outlet (we call this pair "aortic").We carry out the computational analysis with the convective form of the ST-VMS method and the ST-TC described in Sect.7.
In representation of the deforming parts of the domain, we use quadratic NURBS spatially, and cubic NURBS temporally.The mesh is handled with the methods described in Sect.5, including the STNMUM.There is no remeshing in the customary sense of the word.Figure 30 shows the opening/closing stages of the mitral and aortic pairs over the period 1.0 s, prescribed in terms of the angle ratio for each pair.More information on the computational parameters and conditions can be found in [138].
Figures 31 and 32 show, from the preliminary computation reported in [138], the velocity magnitude at different instants.In those figures, the color range from blue to red corresponds to a velocity range from low to high. Figure 33 shows the velocity vectors and pressure around the mitral pair when the pair is closed.The results have all the good features expected from computations with moving-mesh methods, such as pressure jump with zero thickness and boundary-   [138], taking this computation beyond what was reported requires a different way of dealing with the flow at the inlet when we use stress conditions there.This was not attempted in this 2D test computation, because the objective was just to show that the ST-TC method could successfully deal with the TC challenges of this class of problems.

Concluding remarks
We have presented an overview of the new directions we have taken the ST methods in bringing solution and analysis to different classes of computationally challenging engineering problems.Moving in these new directions was motivated mostly by the following three classes of problems we targeted: bio-inspired flapping-wing aerodynamics, windturbine aerodynamics, and cardiovascular fluid mechanics.The new directions for the ST methods include (a) the VMS version of the DSD/SST method, which is called ST-VMS,

Fig. 3 Fig. 4
Fig.3FW control mesh and corresponding surface at three temporalcontrol points

Fig. 8 Fig. 9
Fig. 8 Continuous solution up to t n = 4.0 (top) and its basis functions (bottom)

Fig. 10 Fig. 11
Fig. 10 Hypothetical case of two bars that are initially coinciding, with one hole in the fluid mechanics domain (left).Then the red bar starts moving upward, creating a second hole in the domain (right)

Fig. 12 Fig. 13
Fig. 12 Contraction.The red nodes, 3 and 5, are on the contraction interface and are contacting.The white nodes are the slaves.They are in the same position as their masters, but for visualization purposes we slightly shift their positions in the figure.The numbers indicate the node numbers on the parent mesh.(Color figure online)

Fig. 14 Fig. 15
Fig. 14 Flapping.The ST trajectories of the neighboring ends of the blue and red bars.(Color figure online)

Fig. 16
Fig. 16 Flapping.Blue-bar side of the ST boundary between the two sides.(Color figure online)

Fig. 17
Fig. 17Flapping.Red-bar side of the ST boundary between the two sides.(Color figure online)

Fig. 18
Fig. 17Flapping.Red-bar side of the ST boundary between the two sides.(Color figure online)

Fig. 19 Fig. 20
Fig. 19 Mesh update with the STNMUM in the context of the HW wingtip trajectory, represented with cubic NURBS over four temporal patches.Control point numbering is local to each patch.A control point at the start of a patch and collocated with a control point at the end of the previous patch is in parentheses.A fluid mechanics volume mesh is generated for each patch at the middle control point, and the meshes at the other control points are computed with the mesh moving techniques developed in conjunction with the DSD/SST method.That gives us a temporal representation of the mesh

Fig. 21
Fig. 21 Locust.Vorticity for four equally-spaced points during the second flapping cycle

Fig. 27 Fig. 28
Fig. 27 Vorticity, computed with the ST-VMS method and the STN-MUM

Fig. 29 Fig. 30
Fig. 29 Computational domain.The 2D model was intended to resemble the left ventricle of human heart.The red lines represent the solid surfaces, and the rest of the domain boundaries are the inlet (near the domain center) and outlet (upper left).(Color figure online)

Fig. 31
Fig. 31 Velocity at t = 0.002 s and t = 0.126 s

Fig. 32
Fig. 32 Velocity at t = 0.290 s and t = 0.566 s

Fig. 33
Fig. 33 Two opening/closing pairs.Velocity vectors and pressure (colored) around the mitral pair when the pair is closed (t = 0.126 s).The red lines indicate the zero-thickness surfaces