Space–Time Variational Multiscale Isogeometric Analysis of a tsunami-shelter vertical-axis wind turbine

We present computational flow analysis of a vertical-axis wind turbine (VAWT) that has been proposed to also serve as a tsunami shelter. In addition to the three-blade rotor, the turbine has four support columns at the periphery. The columns support the turbine rotor and the shelter. Computational challenges encountered in flow analysis of wind turbines in general include accurate representation of the turbine geometry, multiscale unsteady flow, and moving-boundary flow associated with the rotor motion. The tsunami-shelter VAWT, because of its rather high geometric complexity, poses the additional challenge of reaching high accuracy in turbine-geometry representation and flow solution when the geometry is so complex. We address the challenges with a space–time (ST) computational method that integrates three special ST methods around the core, ST Variational Multiscale (ST-VMS) method, and mesh generation and improvement methods. The three special methods are the ST Slip Interface (ST-SI) method, ST Isogeometric Analysis (ST-IGA), and the ST/NURBS Mesh Update Method (STNMUM). The ST-discretization feature of the integrated method provides higher-order accuracy compared to standard discretization methods. The VMS feature addresses the computational challenges associated with the multiscale nature of the unsteady flow. The moving-mesh feature of the ST framework enables high-resolution computation near the blades. The ST-SI enables moving-mesh computation of the spinning rotor. The mesh covering the rotor spins with it, and the SI between the spinning mesh and the rest of the mesh accurately connects the two sides of the solution. The ST-IGA enables more accurate representation of the blade and other turbine geometries and increased accuracy in the flow solution. The STNMUM enables exact representation of the mesh rotation. A general-purpose NURBS mesh generation method makes it easier to deal with the complex turbine geometry. The quality of the mesh generated with this method is improved with a mesh relaxation method based on fiber-reinforced hyperelasticity and optimized zero-stress state. We present computations for the 2D and 3D cases. The computations show the effectiveness of our ST and mesh generation and relaxation methods in flow analysis of the tsunami-shelter VAWT.


Introduction
We address the computational challenges encountered in and present results from flow analysis of a vertical-axis wind turbine (VAWT) that has been proposed to also serve as a tsunami shelter. The VAWT model is based on "Life tower" [1]. In addition to the three-blade rotor, the turbine has four support columns at the periphery. The columns support the 3 Mechanical Engineering, Rice University, MS 321, 6100 Main Street, Houston, TX 77005, USA turbine rotor and the shelter. The challenges encountered in computational flow analysis of wind turbines in general include accurate representation of the turbine geometry, multiscale unsteady flow, and moving-boundary flow associated with the rotor motion. The tsunami-shelter VAWT, because of its rather high geometric complexity, poses the additional challenge of reaching high accuracy in turbine-geometry representation and flow solution when the geometry is so complex.
The ST-discretization feature of the integrated method provides higher-order accuracy compared to standard discretization methods. The VMS feature addresses the computational challenges associated with the multiscale nature of the unsteady flow. The moving-mesh feature of the ST framework enables high-resolution computation near the blades. The ST-SI enables moving-mesh computation of the spinning rotor. The mesh covering the rotor spins with it, and the SI between the spinning mesh and the rest of the mesh accurately connects the two sides of the solution. The ST-IGA enables more accurate representation of the blade and other turbine geometries and increased accuracy in the flow solution. The STNMUM enables exact representation of the mesh rotation. The general-purpose NURBS mesh generation method makes it easier to deal with the complex turbine geometry. The quality of the mesh generated with this method is improved with the mesh relaxation method.
The first flow computations for the tsunami-shelter VAWT were presented in [19], in the context of finite element discretization, as test computations with 2D and 3D models to show how the ST-SI method, introduced in the paper, worked. Preliminary test computations with isogeometric discretization were presented in [29] for the 2D model and in [30] for both the 2D and 3D models.

ST-VMS and ST-SUPS
The Deforming-Spatial-Domain/Stabilized ST (DSD/SST) method [31][32][33] was introduced for computation of flows with moving boundaries and interfaces (MBI), including fluid-structure interaction (FSI). In MBI computations the DSD/SST functions as a moving-mesh method. Moving the fluid mechanics mesh to follow an interface enables mesh resolution control near the interface and, consequently, highresolution boundary-layer representation near fluid-solid interfaces. Because the stabilization components of the original DSD/SST are the Streamline-Upwind/Petrov-Galerkin (SUPG) [34] and Pressure-Stabilizing/Petrov-Galerkin (PSPG) [31] stabilizations, it is now called "ST-SUPS." The ST-VMS is the VMS version of the DSD/SST. The VMS components of the ST-VMS are from the residual-based VMS (RBVMS) method [35][36][37][38]. 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 [16,17]), are desirable also in computations without MBI.
The ST-SUPS, ALE-SUPS, RBVMS, ALE-VMS and ST-VMS all have some embedded stabilization parameters that play a significant role (see [43]). These parameters involve a measure of the local length scale (also known as "element length") and other parameters such as the element Reynolds and Courant numbers. 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 [18,19,21,25,106,135,[141][142][143][144][145]. 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 [2,3,5,[7][8][9][146][147][148][149][150][151][152][153]. We will specify which ones we use here in Appendix A.1 and when we describe the computations in Sect. 4.

ST-SI
The ST-SI was introduced in [19], 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 ALE-SI [51,52]. 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 [19]. 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 [20] for the coupled incompressible-flow and thermaltransport 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 [12][13][14][15]19], thermo-fluid analysis of disk brakes [20], flow-driven string dynamics in turbomachinery [14,15,[128][129][130], flow analysis of turbocharger turbines [22,26,27,131,132], flow around tires with road contact and deformation [122,[133][134][135][136], fluid films [136,137], aerodynamic analysis of ram-air parachutes [29,30,138], and flow analysis of heart valves [78,79,118,[123][124][125] and ventricle-valve-aorta sequences [120].
In the ST-SI version presented in [19] the SI is between a thin porous structure and the fluid on its two sides. This enables dealing with the porosity in a fashion consistent with how the standard fluid-fluid SIs are dealt with and how the Dirichlet conditions are enforced weakly with fluid-solid SIs. This version also enables handling thin structures that have T-junctions. This method has been applied to incompressibleflow aerodynamic analysis of ram-air parachutes with fabric porosity [29,30,138]. The compressible-flow ST-SI methods were introduced in [139], 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 [139]. These, together with the compressibleflow ST SUPG method [154], 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 [139,140].

ST-IGA and STNMUM
The success with IGA basis functions in space [40,51, 66,155] motivated the integration of the ST methods with isogeometric discretization, which is broadly called "ST-IGA." The ST-IGA was introduced in [16]. Computations with the ST-VMS and ST-IGA were first reported in [16] 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 [16] 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 [16,17] and demonstrated in [21,23,24]. 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 STN-MUM [21,23,24], with the name coined in [25]. The STN-MUM 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" [16,17,21,43] 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 [18,20,122,[128][129][130]156]. 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 [21,23,43,109], bioinspired MAVs [24,107,108,110] and wing-clapping [111,112], separation aerodynamics of spacecraft [101], aerodynamics of horizontal-axis [11,25,107,108] and vertical-axis [12][13][14][15]19] wind-turbines, thermo-fluid analysis of ground vehicles and their tires [18,29,122], thermo-fluid analysis of disk brakes [20], flow-driven string dynamics in turbomachinery [14,15,[128][129][130], flow analysis of turbocharger turbines [22,26,27,131,132], and flow analysis of coronary arteries in motion [126].

General-purpose NURBS mesh generation method
To make the ST-IGA use, and in a wider context the IGA use, even more practical in computational flow analysis with complex geometries, NURBS volume mesh generation needs to be easier and more automated. To that end, a generalpurpose NURBS mesh generation method was introduced in [26]. The method is based on multi-block-structured mesh generation with existing techniques, projection of that mesh to a NURBS mesh made of patches that correspond to the blocks, and recovery of the original model surfaces. The recovery of the original surfaces is to the extent they are suitable for accurate and robust fluid mechanics computations. The method is expected to retain the refinement distribution and element quality of the multi-block-structured mesh that we start with. Because there are ample good techniques and software for generating multi-block-structured meshes, the method makes general-purpose mesh generation relatively easy. Mesh-quality performance studies for 2D and 3D meshes, including those for complex models, were presented in [27]. A test computation for a turbocharger turbine and exhaust manifold was also presented in [27], with a more detailed computation in [131], and with additional computational analysis in [132]. The mesh generation method was used also in the pump-flow analysis part of the flow-driven string dynamics presented in [14,15,130]. The performance studies and test computations demonstrated that the generalpurpose NURBS mesh generation method makes the IGA use in fluid mechanics computations even more practical.

ST-SI-IGA
An SI provides mesh generation flexibility in a general context by accurately connecting the two sides of the solution computed over nonmatching meshes. This type of mesh generation flexibility is especially valuable in complex-geometry flow computations with isogeometric discretization, removing the matching requirement between the NURBS patches without loss of accuracy. This feature was used in the flow analysis of heart valves [78,79,118,[123][124][125], ventriclevalve-aorta sequences [120], turbocharger turbines [22,26,27,131,132], spacecraft parachute compressible-flow analysis [140], and flow around tires with road contact and deformation [135,136]. It is used also in the VAWT flow analysis presented here.

Mesh relaxation based on fiber-reinforced hyperelasticity and optimized ZSS
isogeometric discretization but are also applicable to finite element discretization. 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 constitutive models and parameters can be defined individually for the elements or mesh regions. For more on the method, see [28].

Computations presented
We present computations for the 2D and 3D cases. In the 2D case, we present computations for three different meshes, two different time-step sizes, and two different tip-speed ratios.
In the 3D case, we compute with one of the combinations of the 2D case and compare the 2D and 3D results.

ST-VMS and ST-SI
For completeness, we include, mostly from [19,133], the ST-VMS and ST-SI methods.

ST-VMS
The ST-VMS is given as where are the residuals of the momentum equation and incompressibility constraint. Here, ρ, u, p, f, and h are the density, velocity, pressure, body force, and the stress specified at the boundary. The stress tensor is defined as σ σ σ = −pI+2με ε ε(u), where I is the identity tensor, μ = ρν is the viscosity, ν is the kinematic viscosity, and ε ε ε (u) = (∇ ∇ ∇u) + (∇ ∇ ∇u) T /2 is the strain-rate tensor. The test functions associated with the u and p are w and q. A superscript "h" indicates that the function is coming from a finite-dimensional space. The symbol Q n represents the ST slice between time levels n and n + 1, (P n ) h is the part of the slice lateral boundary associated with the boundary condition h, and n is the spatial domain at time level n. The superscript "e" is the ST element counter, and n el is the number of ST elements. The functions are discontinuous in time at each time level, and the superscripts "−" and "+" indicate the values of the functions just below and above the time level.

Remark 1
The ST-SUPS can be obtained from the ST-VMS by dropping the eighth and ninth integrations.
The stabilization parameters, τ SUPS and ν LSIC , will be given in Appendix A.1.

ST-SI
In describing the ST-SI, labels "Side A" and "Side B" will represent the two sides of the SI. The ST-SI version of the formulation given by Eq. (1) Fig. 1 Tsunami-shelter VAWT model. The disk-shaped structure just below the rotor is the shelter. The four columns at the periphery support the turbine rotor and the shelter where Here, (P n ) SI is the SI in the ST domain, n is the unit normal vector, v is the mesh velocity, γ = 1, and C is a nondimensional constant. The element length h will be defined in Appendix A.2. The four support columns are cylindrical with circular cross-section, and they provide enough strength to support the rotor, which is estimated to weigh 3 t, and the shelter. We carry out the computations at a constant free-stream velocity U ∞ and with prescribed rotor motion at constant angular velocity. The rotation is clockwise viewed from the top. The air density and kinematic viscosity are 1.205 kg/m 3 and 1.511×10 −5 m 2 /s. We extract from the computations the instantaneous power coefficient C POW , defined as

Tsunami-shelter VAWT model
where A and P are the projected area of the wind turbine and the power generated. We report the power coefficient as a function of the blade orientation as represented by the angle φ seen in Fig. 2. With that orientation, the flow speed seen by a blade can be calculated as where λ is the tip-speed ratio (TSR). The symbol T will denote the rotation cycle.

Computations
We present computations for the 2D and 3D cases. The computational-domain size in the wind direction is 62.5 times the rotor diameter, with a distance of 18.75 times the rotor diameter between the upstream boundary and the rotor center. In the cross-wind direction, the domain size is 37.5 times the rotor diameter. In the 3D case, the domain height is 10 times the rotor diameter. The mesh position is represented by quadratic NURBS in time. There are three patches that are 120 • each, and the secondary mapping introduced in [21] is used to achieve the constant angular velocity. We set U ∞ = 12.56 m/s. In the 2D case, the time-step sizes used correspond to Δφ = 1 • and 2 • . In the 3D case, we use only one timestep size, and that corresponds to Δφ = 1 • . The number of nonlinear iterations per time step is 4, and the number of GMRES iterations per nonlinear iteration is 200. The first two nonlinear iterations are based on the ST-SUPS, and the last two the ST-VMS. The stabilization parameters are those given by Eqs. (10)- (12) and (21)- (22). In the ST-SI (see Eq. (4)), we set C = 2.

2D case
We first compute with TSR = 2. The model geometry and the SI are shown in Fig. 3. The boundary conditions are U ∞ at the inflow, zero stress at the outflow, slip at the lateral boundaries, and no-slip on the rotor and support column surfaces. The prescribed velocities are evaluated at the time integration points, with the values extracted from the NURBS representation of the prescribed motion. We use three different meshes. We start with Mesh 1, and obtain the other two meshes by knot insertion. We halve the knot spacing to get Mesh 2, and halve it again to get Mesh 3. Figure 4 shows Mesh 3. The number of control points and elements are shown in Table 1. We compute for 10 rotations. Figure 5 shows, for Δφ = 2 • and 1 • , C POW , averaged over the three blades in the last three rotations. The results from  Figures 6 and 7 show, for Mesh 1 with Δφ = 1 • and Mesh 2 with Δφ = 2 • , the velocity magnitude in the wake of the support columns located at φ = 180 • and φ = 90 • . Mesh 1, with even Δφ = 1 • , is not able to capture the wake as well as Mesh 2 does even with Δφ = 2 • . This indicates that a reasonable level of mesh refinement is needed. Table 2 shows C POW averaged over the three blades in the last three rotations and the standard deviation (σ CPOW ) associated with that, both averaged over the rotation cycle. We believe the higher σ CPOW values we see at higher Courant numbers are caused by strong vortices, which we do not believe to be representative of what actually happens in 3D.    We also compute with TSR = 3. We use Mesh 3 with Δφ = 1 • . Figure 8 shows C POW , averaged over the three blades in the last three rotations, for TSR = 3 and 2. Rotationcycle-averaged value of C POW is 0.18, which is significantly larger than what we have for TSR = 2, translates to a total C POW value of 0.54.

3D case
We compute with TSR = 3. The boundary conditions are noslip on all turbine surfaces and the bottom boundary, U ∞ at the inflow, zero stress at the outflow, and slip at the lateral and top boundaries. All prescribed velocities are evaluated at the time integration points with the values extracted from the NURBS representation of the prescribed motion. Figure 9 shows the initial mesh, generated with the method in [26]. There are 1,544,460 control points and 955,477 quadratic NURBS elements. There are total nine SIs in the mesh. Figure 10 shows the SIs with an actual slip. Figures 11  and 12 show the SIs that are just for mesh generation purposes.
To improve the mesh quality, we use the mesh relaxation method introduced in [28], which is based on fiber-reinforced hyperelasticity and optimized ZSS. The hyperelasticity  50) and (51) in [28], with κ = 0. The material properties are κ B = 10 −2 , β B = 2, μ = 1, C 1 = 100, and C 2 = 10 −2 . We ramp-up to the optimized ZSS in 10 equal increments. For the numerical integration, 4×4×4 quadrature points are used. Figures 13  and 14 show the initial mesh and the mesh after the relaxation. We clearly see the improvement in the mesh-line orthogonality, without much change in the element aspect ratio or the local resolution.
We compute for two rotations and report the solution from the last rotation. Figure 15 shows the isosurfaces of the second invariant of the velocity gradient tensor near a blade, at different positions of the blade. Figure 16 shows C POW . The total C POW averaged over the rotation cycle is about 0.56. Figure 17 shows C POW from the 2D and 3D computations, averaged over the three blades. They are in reasonably good agreement, considering that the space dimensions are not the same. The 3D value is slightly higher in an average sense. Because the shelter and the upper part of the frame divert some of the flow to the rotor, the blades see more flow than what corresponds to the projected area calculated using the blade height.

Concluding remarks
We have presented computational flow analysis of a VAWT that has been proposed to also serve as a tsunami shelter. The turbine has, in addition to the three-blade rotor, four support columns at the periphery, which support the turbine rotor and the shelter. The computational challenges encountered in flow analysis of the tsunami-shelter VAWT include those encountered in flow analysis of wind turbines in gen- eral, such as accurate representation of the turbine geometry, multiscale unsteady flow, and moving-boundary flow associated with the rotor motion. Beyond those, because of its rather high geometric complexity, the tsunami-shelter VAWT poses the challenge of reaching high accuracy in turbine-geometry representation and flow solution when the geometry is so complex. We address the challenges with an ST computational method that integrates three special ST methods around a core ST method and mesh-related methods. The core method is the ST-VMS, which subsumes its precursor ST-SUPS. The three special methods are the ST-SI, ST-IGA, and the STNMUM. The mesh-related methods are a general-purpose NURBS mesh generation method and a mesh relaxation method based on fiber-reinforced hyperelasticity and optimized ZSS.
The ST-discretization feature of the integrated method provides higher-order accuracy compared to standard discretization methods. The VMS feature addresses the computational challenges associated with the multiscale nature of the unsteady flow. The moving-mesh feature of the ST framework enables high-resolution computation near the blades. The ST-SI enables moving-mesh computation of the spinning rotor. The mesh covering the rotor spins with it, and the SI between the spinning mesh and the rest of the mesh accurately connects the two sides of the solution. The ST-IGA enables more accurate representation of the blade and other turbine geometries and increased accuracy in the flow solution. The STNMUM enables exact representation of the mesh rotation. The general-purpose NURBS mesh generation method makes it easier to deal with the complex turbine geometry. The quality of the mesh generated with this method is improved with the mesh relaxation method.
We have presented computations for the 2D and 3D cases. In the 2D case, we have presented computations for three different meshes, two different time-step sizes, and two different tip-speed ratios. In the 3D case, we computed with one of the combinations of the 2D case and compared the 2D and 3D results. The computations show the effectiveness of our ST and mesh generation and relaxation methods in flow analysis of the tsunami-shelter VAWT. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.

A.1 ST-VMS
There are various ways of defining the stabilization parameters τ SUPS and ν LSIC . Here, τ SUPS is mostly from [142]: The first and second components are given as and where r is the solution-gradient direction: Here G ST and G are the ST and space-only element metric tensors: wherê The ST and space-only Jacobian tensors are and where θ and ξ ξ ξ are the temporal and spatial parametric coordinates. The transformation tensor D ST is defined as The definitions used for D θ and D play an important role, especially for higher-order isogeometric discretization [142,144] and simplex elements [143].
In this article, we set D θ = 1 and set D to its "RQD-MAX" version [144].

A.2 ST-SI
The element length used in the ST-SI is given as Remark 2 We note that the expression for h given by Eq. (24) is slightly different from its original form introduced in [131]: The modification, which has only minor effect, is for consistency with how the stabilization parameters are calculated from their components and for implementation convenience.