Space–time VMS isogeometric analysis of the Taylor–Couette flow

Abstract

The Taylor–Couette flow is a classical fluid mechanics problem that exhibits, depending on the Reynolds number, a range of flow patterns, with the interesting ones having small-scale structures, and sometimes even wavy nature. Accurate representation of these flow patterns in computational flow analysis requires methods that can, with a reasonable computational cost, represent the circular geometry accurately and provide a high-fidelity flow solution. We use the Space–Time Variational Multiscale (ST-VMS) method with ST isogeometric discretization to address these computational challenges and to evaluate how the method and discretization perform under different scenarios of computing the Taylor–Couette flow. We conduct the computational analysis with different combinations of the Reynolds numbers based on the inner and outer cylinder rotation speeds, with different choices of the reference frame, one of which leads to rotating the mesh, with the full-domain and rotational-periodicity representations of the flow field, with both the convective and conservative forms of the ST-VMS, with both the strong and weak enforcement of the prescribed velocities on the cylinder surfaces, and with different mesh refinements. The ST framework provides higher-order accuracy in general, and the VMS feature of the ST-VMS addresses the computational challenges associated with the multiscale nature of the flow. The ST isogeometric discretization enables exact representation of the circular geometry and increased accuracy in the flow solution. In computations where the mesh is rotating, the ST/NURBS Mesh Update Method, with NURBS basis functions in time, enables exact representation of the mesh rotation, in terms of both the paths of the mesh points and the velocity of the points along their paths. In computations with rotational-periodicity representation of the flow field, the periodicity is enforced with the ST Slip Interface method. With the combinations of the Reynolds numbers used in the computations, we cover the cases leading to the Taylor vortex flow and the wavy vortex flow, where the waves are in motion. Our work shows that all these ST methods, integrated together, offer a high-fidelity computational analysis platform for the Taylor–Couette flow and for other classes of flow problems with similar features.

1 Introduction

We use the Space–Time Variational Multiscale (ST-VMS) method [1, 2] with ST isogeometric discretization [1, 3] to address the computational challenges encountered in the Taylor–Couette flow, a classical fluid mechanics problem. The Taylor–Couette flow exhibits, depending on the Reynolds number, a range of flow patterns, with the interesting ones having small-scale structures, and sometimes even wavy nature. Accurate representation of these flow patterns in computational flow analysis requires methods that can, with a reasonable computational cost, represent the circular geometry accurately and provide a high-fidelity flow solution. We evaluate how the ST-VMS with ST isogeometric discretization perform under different scenarios of computing the Taylor–Couette flow.

The ST framework provides higher-order accuracy in general, and the VMS feature of the ST-VMS addresses the computational challenges associated with the multiscale nature of the flow. The ST isogeometric discretization enables exact representation of the circular geometry and increased accuracy in the flow solution. In computations where the reference-frame choice leads to rotating the mesh, the ST/NURBS Mesh Update Method (STNMUM) [4, 5], with NURBS basis functions in time, enables exact representation of the mesh rotation, in terms of both the paths of the mesh points and the velocity of the points along their paths. In computations with rotational-periodicity representation of the flow field, the periodicity is enforced with the ST Slip Interface (ST-SI) method [6].

The first reported computation of the Taylor–Couette flow in the ST framework goes back to 1993 [7], where the computations were carried out with the Deforming-Spatial-Domain/Stabilized ST (DSD/SST) method [8] and finite element discretization. The stabilization components of the DSD/SST constituted a VMS precursor. The Reynolds numbers used in the computations led to the Taylor vortex flow and the wavy vortex flow, where the waves are in motion. The first reported computation with the isogeometric discretization, in the semi-discrete framework, was in 2010 [9]. We see advantages in computing the Taylor–Couette flow, and other classes of flow problems that have similar features, with the isogeometric discretization in the ST framework, for reasons summarized in the previous paragraph and explained more in the later parts of the article.

1.1 Stabilized and VMS ST computational methods

The stabilized and VMS ST computational methods started with the DSD/SST, which was introduced for computation of flows with moving boundaries and interfaces (MBI), including fluid–structure interaction (FSI). In flow computations with MBI, 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, high-resolution boundary-layer representation near fluid–solid interfaces.

Stabilized and VMS methods have for decades been playing a core-method role in flow analysis with semi-discrete and ST computational methods. The incompressible-flow Streamline-Upwind/Petrov-Galerkin (SUPG) [10, 11] and compressible-flow SUPG [12,13,14] methods are two of the earliest and most widely used stabilized methods. The Pressure-Stabilizing/Petrov-Galerkin (PSPG) method [8], with its Stokes-flow version introduced in [15], is also among the earliest and most widely used. These methods bring numerical stability in computation of flow problems at high Reynolds or Mach numbers and when using equal-order basis functions for velocity and pressure in incompressible flows. Because the methods are residual-based, the stabilization is achieved without loss of accuracy. The residual-based VMS (RBVMS) [9, 16,17,18], which is also widely used now, subsumes its precursor SUPG/PSPG.

Because the stabilization components of the original DSD/SST are the SUPG and PSPG stabilizations, it is now also called “ST-SUPS.” The ST-VMS is the VMS version of the DSD/SST. The VMS components of the ST-VMS are from the RBVMS. The ST-VMS, which subsumes its precursor ST-SUPS, 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 [1, 19]), are desirable also in computations without MBI.

As a moving-mesh method, the DSD/SST is an alternative to the Arbitrary Lagrangian–Eulerian (ALE) method, which is older (see, for example, [20]) and more commonly used. The ALE-VMS method [21, 22] is the VMS version of the ALE. It succeeded the ST-SUPS and ALE-SUPS [23] and preceded the ST-VMS. 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 [23], wind turbine aerodynamics and FSI [24,25,26,27,28,29,30,31,32,33], more specifically, vertical-axis wind turbines (VAWTs) [31, 34, 35], floating wind turbines [36], wind turbines in atmospheric boundary layers [30, 31, 37,38,39], and fatigue damage in wind turbine blades [40], patient-specific cardiovascular fluid mechanics and FSI [21, 41,42,43,44,45,46], biomedical-device FSI [47,48,49,50,51,52,53,54], ship hydrodynamics with free-surface flow and fluid–object interaction [55, 56], hydrodynamics and FSI of a hydraulic arresting gear [57, 58], hydrodynamics of tidal-stream turbines with free-surface flow [59], bioinspired FSI for marine propulsion [60, 61], bridge aerodynamics and fluid–object interaction [62,63,64], and mixed ALE-VMS/Immersogeometric computations [50,51,52, 65, 66] in the framework of the Fluid–Solid Interface-Tracking/Interface-Capturing Technique. Recent advances in stabilized and multiscale methods may be found for stratified incompressible flows in [67], for divergence-conforming discretizations of incompressible flows in [68], and for compressible flows with emphasis on gas-turbine modeling in [69].

In flow computations with FSI or MBI, the ST-SUPS and ST-VMS require a mesh moving method. 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 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. A number of well-performing mesh moving methods were developed in conjunction with the ST-SUPS and ST-VMS. The first one, introduced in [7, 70], was the Jacobian-based stiffening, which is now called, for reasons explained in [71], “mesh-Jacobian-based stiffening.” The most recent ones are the element-based mesh relaxation [72], where the mesh motion is determined by using the large-deformation mechanics equations and an element-based zero-stress-state (ZSS), a mesh moving method [73] based on fiber-reinforced hyperelasticity and optimized ZSS, and a linear-elasticity-based mesh moving method with no cycle-to-cycle accumulated distortion [71, 74].

The ST-SUPS and ST-VMS have also been applied to many classes of FSI, MBI and fluid mechanics problems (see [75] for a comprehensive summary of the work prior to July 2018). The classes of problems include spacecraft parachute analysis for the landing-stage parachutes [22, 72], cover-separation parachutes [76] and the drogue parachutes [77], wind turbine aerodynamics for horizontal-axis wind turbine (HAWT) rotors [22], full HAWTs [78] and VAWTs [6, 31,32,33, 79], flapping-wing aerodynamics for an actual locust [4, 5, 22], bioinspired MAVs [80] and wing-clapping [81, 82], blood flow analysis of cerebral aneurysms [83], stent-blocked aneurysms [84, 85], aortas [53, 54, 86,87,88], heart valves [53, 54, 81, 87, 89,90,91,92], ventricle-valve-aorta sequences [71], and spacecraft aerodynamics [93], thermo-fluid analysis of ground vehicles and their tires [2, 38, 39, 90], thermo-fluid analysis of disk brakes [94], flow-driven string dynamics in turbomachinery [32, 33, 95, 96], flow analysis of turbocharger turbines [3, 97, 98], flow around tires with road contact and deformation [90, 99,100,101], fluid films [102], ram-air parachutes [38, 39, 103], and compressible-flow spacecraft parachute aerodynamics [104, 105].

1.2 ST-SI

The ST-SI was introduced in [6] 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 [106, 107]. 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 [6]. 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 [94] 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 [6, 31,32,33, 79], thermo-fluid analysis of disk brakes [94], flow-driven string dynamics in turbomachinery [32, 33, 95, 96], flow analysis of turbocharger turbines [3, 97, 98], flow around tires with road contact and deformation [90, 99,100,101], fluid films [102], aerodynamic analysis of ram-air parachutes [38, 39, 103], and flow analysis of heart valves [53, 54, 81, 87, 89,90,91,92] and ventricle-valve-aorta sequences [71]. In the ST-SI version presented in [6] 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 [38, 39, 103].

1.3 ST Isogeometric Analysis

The success with Isogeometric Analysis (IGA) basis functions in space [21, 41, 106, 108] motivated the integration of the ST methods with isogeometric discretization, which is broadly called “ST-IGA.” The ST-IGA was introduced in [1]. Computations with the ST-VMS and ST-IGA were first reported in [1] 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 [1] 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 [1, 19] and demonstrated in [4, 5]. 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 STNMUM [4, 5], with the name coined in [78]. The STNMUM has a wide scope that includes spinning solid surfaces. With the spinning motion represented by quadratic NURBS in time, and with sufficient number of temporal patches for a full rotation, the circular paths are represented exactly. A “secondary mapping” [1, 4, 19, 22] 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 [2, 90, 94,95,96, 109]. 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 [4, 5, 22], bioinspired MAVs [80] and wing-clapping [81, 82], separation aerodynamics of spacecraft [76], aerodynamics of horizontal-axis [22, 78] and vertical-axis [6, 31,32,33, 79] wind turbines, thermo-fluid analysis of ground vehicles and their tires [2, 38, 39, 90], thermo-fluid analysis of disk brakes [94], flow-driven string dynamics in turbomachinery [32, 33, 95, 96], and flow analysis of turbocharger turbines [3, 97, 98].

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 [3, 97, 98], flow-driven string dynamics in turbomachinery [32, 33, 95, 96], ram-air parachutes [38, 39, 103], spacecraft parachutes [105], aortas [53, 54, 87, 88], heart valves [53, 54, 87, 91, 92], ventricle-valve-aorta sequences [71], tires with road contact and deformation [99,100,101], fluid films [102], and VAWTs [6, 31,32,33, 79]. The image-based arterial geometries used in patient-specific arterial FSI computations do not come from the ZSS of the artery. Using IGA basis functions in space is now a key part of some of the newest ZSS estimation methods [53, 110,111,112] and related shell analysis [113]. The IGA has also been successfully applied to the structural analysis of wind turbine blades [114,115,116,117,118].

1.4 Stabilization parameters and element lengths

In all the semi-discrete and ST stabilized and VMS methods discussed in Sect. 1.1, an embedded stabilization parameter, known as “$\tau$,” plays a significant role (see [22]). This parameter involves a measure of the local length scale (also known as “element length”) and other parameters such as the element Reynolds and Courant numbers. The interface terms in the ST-SI also involve element length, in the direction normal to the interface. Various element lengths and $\tau$s were proposed, starting with those in [10, 11] and [12,13,14], followed by the ones introduced in [119, 120]. In many cases, the element length was seen as an advection length scale, in the flow-velocity direction. The $\tau$ definition introduced in [120], which is for the advective limit and is now called “$\tau _{\mathrm {SUGN1}}$” (and the corresponding element length is now called “$h_{\mathrm {UGN}}$”), automatically yields lower values for higher-order finite element basis functions.

Calculating the $\tau$s based on the element-level matrices and vectors was introduced in [121] in the context of the advection–diffusion equation and the Navier–Stokes equations of incompressible flows. These definitions are expressed in terms of the ratios of the norms of the matrices or vectors. They automatically take into account the local length scales, advection field and the element Reynolds number. The definitions based on the element-level vectors were shown [121] to address the difficulties reported at small time-step sizes. A second element length scale, in the solution-gradient direction and called “$h_{\mathrm {RGN}}$,” was introduced in [122, 123]. Recognizing this as a diffusion length scale, a new stabilization parameter for the diffusive limit, “$\tau _{\mathrm {SUGN3}}$,” was introduced in [123, 124], to be used together with $\tau _{\mathrm {SUGN1}}$ and “$\tau _{\mathrm {SUGN2}}$,” the parameters for the advective and transient limits. For the stabilized ST methods, “$\tau _{\mathrm {SUGN12}}$,” representing both the advective and transient limits, was also introduced in [122, 123].

Some new options for the stabilization parameters used with the SUPS and VMS were proposed in [2, 4, 78, 125]. These include a fourth $\tau$ component, “$\tau _{\mathrm {SUGN4}}$” [2], which was introduced for the VMS, considering one of the two extra stabilization terms the VMS has compared to the SUPS. They also include stabilization parameters [2] for the thermal-transport part of the VMS for the coupled incompressible-flow and thermal-transport equations.

Some of the stabilization parameters described in this subsection were also used in computations with other SUPG-like methods, such as the computations reported in [126, 127].

The stabilization parameters and element lengths discussed in this subsection so far were all originally intended for finite element discretization but quite often used also for isogeometric discretization. The element lengths and stabilization parameters introduced in [128] target isogeometric discretization but are also applicable to finite element discretization. They were introduced in the context of the advection–diffusion equation and the Navier–Stokes equations of incompressible flows. The direction-dependent element length expression was outcome of a conceptually simple derivation. The key components of the derivation are mapping the direction vector from the physical ST element to the parent ST element, accounting for the discretization spacing along each of the parametric coordinates, and mapping what has been obtained in the parent element back to the physical element. The test computations presented in [128] for pure-advection cases, including those with discontinuous solution, showed that the element lengths and stabilization parameters proposed result in good solution profiles. The test computations also showed that the “UGN” parameters give reasonably good solutions even with NURBS basis functions. The stabilization parameters given in [100], which were mostly from [128], were the latest ones designed in conjunction with the ST-VMS.

In general, we decide what parametric space to use based on reasons like numerical integration efficiency or implementation convenience. Obviously, choices based on such reasons should not influence the method in substance. We require the element lengths, including the direction-dependent element lengths, to have node-numbering invariance for all element types, including simplex elements. The direction-dependent element length expression introduced in [129] meets that requirement. This is accomplished by using in the element length calculations for simplex elements a preferred parametric space instead of the standard integration parametric space. The element length expressions based on the two parametric spaces were evaluated in [129] in the context of simplex elements. It was shown that when the element length expression is based on the integration parametric space, the variation with the node numbering could be by a factor as high as 1.9 for 3D elements and 2.2 for ST elements. It was also shown that the element length expression based on the integration parametric space could overestimate the element length by a factor as high as 2.8 for 3D elements and 3.2 for ST elements.

Targeting B-spline meshes for complex geometries, new direction-dependent element length expressions were introduced in [130]. These latest element length expressions are outcome of a clear and convincing derivation and more suitable for element-level evaluation. The new expressions are based on a preferred parametric space, instead of the standard integration parametric space, and a transformation tensor that represents the relationship between the integration and preferred parametric spaces. We do not want the element splitting to influence the actual discretization, which is represented by the control or nodal points. Therefore, the local length scale should be invariant with respect to element splitting. That invariance is a crucial requirement in element definition, because unlike the element definition choices based on implementation convenience or computational efficiency, it influences the solution. It was proven in [131] that the local-length-scale expressions introduced in [130] meet that requirement.

The direction-dependent local-length-scale expressions introduced in [128, 130] have been used in computational flow analysis of turbocharger turbines [98], compressible-flow spacecraft parachutes [105], tires with road contact and deformation [100, 101], fluid films [102], ventricle-valve-aorta sequences [71], and tsunami-shelter VAWTs [79].

1.5 Taylor–Couette flow

We conduct the computational analysis with different combinations of the Reynolds numbers based on the inner and outer cylinder rotation speeds, with different choices of the reference frame, one of which leads to rotating the mesh, with the full-domain and rotational-periodicity representations of the flow field, with both the convective and conservative forms of the ST-VMS, with both the strong and weak enforcement of the prescribed velocities on the cylinder surfaces, and with different mesh refinements.

With the combinations of the Reynolds numbers used in the computations, we cover the cases leading to the Taylor vortex flow and the wavy vortex flow, where the waves are in motion. The computations show that all these ST methods, integrated together, offer a high-fidelity computational analysis platform for the Taylor–Couette flow and for other classes of flow problems with similar features.

1.6 Outline of the remaining sections

The computations are presented in Sect. 2, and the concluding remarks are given in Sect. 3. The ST-VMS and ST-SI methods and the stabilization parameters are described in Appendices A and B. The methods for the exact representation of the mesh rotation and rotation-generated velocities are described in Appendix C.

2 Taylor–Couette flow

2.1 Problem setup

Fig. 1

Flow domain. We use the local orthonormal basis set ${\mathbf {e}}_r$, ${\mathbf {e}}_\theta$ and ${\mathbf {e}}_z$ associated with the cylindrical coordinate system

The flow domain is shown in Fig. 1. The inner and outer cylinder radii are $r_\mathrm {i}$ and $r_\mathrm {o}$, and $\eta = \frac{r_\mathrm {i}}{r_\mathrm {o}}$ is one of the parameters representing the flow conditions. We use the local orthonormal basis set ${\mathbf {e}}_r$, ${\mathbf {e}}_\theta$ and ${\mathbf {e}}_z$ associated with the cylindrical coordinate system, and $u_r$, $u_\theta$ and $u_z$ denote the components of the velocity ${\mathbf {u}}$. The angular velocity vector is given as $\pmb {\omega } = \omega {\mathbf {e}}_{z}$. From that, the inner and outer surface speeds are $U_\mathrm {i}=\omega _\mathrm {i} r_\mathrm {i}$ and $U_\mathrm {o}=\omega _\mathrm {o} r_\mathrm {o}$. We define two Reynolds numbers, $\mathrm {Re}_\mathrm {i}=\frac{U_\mathrm {i}(r_\mathrm {o}-r_\mathrm {i})}{\nu }$ and $\mathrm {Re}_\mathrm {o}=\frac{U_\mathrm {o}(r_\mathrm {o}-r_\mathrm {i})}{\nu }$, where $\nu =\frac{\mu }{\rho }$ is the kinematic viscosity, $\rho$ is the density, and $\mu$ is the dynamic viscosity.

Remark 1

We note that, depending on the rotation direction, a Reynolds number could be negative.

The Taylor number is

$$\begin{aligned} \mathrm {Ta}&= \frac{\omega _\mathrm {i}^2 r_\mathrm {i} \left( r_\mathrm {o} - r_\mathrm {i}\right) ^3}{\nu ^2} \end{aligned}$$

(1)

$$\begin{aligned}&= \mathrm {Re}_\mathrm {i}^2 \left( \frac{1}{\eta } - 1\right) . \end{aligned}$$

(2)

The time scale is defined as $T = \frac{2 \pi }{\omega _\mathrm {i} - \omega _\mathrm {o}}$.

Fig. 2

Flow patterns depending on the two Reynolds numbers. The blue, red and gray mark the regions of wavy vortex flow, Taylor vortex flow and Couette flow. The red points mark the cases. The negative Reynolds number represents clockwise rotation of the outer cylinder. The figure has been adapted from [132]

In our computations, $\eta = 0.883$, which has been studied experimentally in [132]. Figure 2 shows the flow patterns depending on the two Reynolds numbers. The figure has been adapted from [132]. Table 1 shows the combinations of the Reynolds numbers in our computations, all with $\mathrm {Re}_\mathrm {i}- \eta \mathrm {Re}_\mathrm {o} = 250$, making the relative angular speed the same in all combinations.

Table 1 Combinations of the Reynolds numbers in the computations, all with $\mathrm {Re}_\mathrm {i}- \eta \mathrm {Re}_\mathrm {o} = 250$, and the corresponding Taylor numbers and flow patterns

In the computations, we attach the mesh to the inner cylinder to test the methods, the effect of the reference frame selected, and the effect of using a moving mesh. The outer surface motion is represented by the prescribed-velocity condition. We have axial periodicity over the length $\pi {r_\mathrm {i}}$. We consider two choices for the reference frame, the inertial reference frame (IRF) and the noninertial reference frame (NRF) rotating with the inner cylinder. We compute with both the full-domain and rotational-periodicity representations of the flow field.

2.2 Meshes

We represent the circular problem geometry exactly by using quadratic NURBS patches. However, we cannot represent a circle with $C^1$-continuous functions. In the computations, we use patches that have quarter, one-third and half-domain sizes. In the rotational-periodicity computations, we use those patches. In the full-domain computations, we use two half-domain patches and connect them in a $C^0$-continuous fashion. Figure 3 shows how the circular geometry represented by quadratic NURBS patches with the minimum number of elements.

Fig. 3

Representation of the circular geometry by using quadratic NURBS patches with the minimum number of elements. Patches that have quarter, one-third and half-domain sizes

Remark 2

We note that in the patch that has half-domain size, we need two elements, otherwise the NURBS weights would have negative values.

With these and by using the knot-insertion technique, we obtain the meshes shown in Fig. 4, which we call “coarse.” For the patch that has half-domain size, we refine twice to obtain the meshes that we call “medium” and “fine” (see Fig. 5). The full-domain coarse mesh is shown in Fig. 6. The number of control points and elements are given in Table 2.

Fig. 4

Coarse control meshes. Patches that have quarter, one-third and half-domain sizes. The red circles represent the control points

Fig. 5

Medium and fine control meshes for the patch that has half-domain size. The red circles represent the control points

Fig. 6

Full-domain coarse control mesh. The red circles represent the control points

Table 2 Number of control points (nc) and elements (ne) in the radial, circumferential and axial directions

2.3 Boundary conditions and mesh motions

The boundary conditions are no-slip on both the inner and outer surfaces and periodic in the axial direction. The no-slip conditions are enforced strongly or weakly, depending on the test computation, and the periodicity is enforced with the ST-SI. In computations with the IRF, $U_\mathrm {i} {\mathbf {e}}_{\theta }$ and $U_\mathrm {o} {\mathbf {e}}_{\theta }$ are prescribed on the inner and outer surfaces, and the mesh is rotating with $\pmb {\omega }_\mathrm {i}$. In computations with the NRF, zero and $\left( U_\mathrm {o} -U_\mathrm {i} \frac{1}{\eta } \right) {\mathbf {e}}_\theta$ are prescribed on the inner and outer surfaces, and the mesh is stationary. The Coriolis and centrifugal forces are added to the governing equations based on the steady rotation speed of $\omega _\mathrm {i}$.

Remark 3

The velocity relative to the mesh is comparable between all computations.

2.4 Computational conditions

We use the convective or conservative form of the ST-VMS (see Appendix A), depending on the test computation. The stabilization parameters and elements are given in Appendix B. The mesh motion is represented by quadratic NURBS in time (see Appendix C), and the flow solution is based on linear functions in time. The time-step size, $\varDelta t$, is $\frac{T}{60}$, $\frac{T}{120}$ and $\frac{T}{240}$ for the coarse, medium and fine meshes. Table 3 shows the mesh rotation per time step when the reference frame choice is the IRF. The rotation per time step has a computational significance related to the accuracy in representing the prescribed velocity in the flow solution, and that will become clear in Remark 7. The number of nonlinear iterations per time step is 3 (except in Sect. 2.5.5, where it is 5), and the number of GMRES iterations per nonlinear iteration is 500.

Table 3 Mesh rotation per time step when the reference frame choice is the IRF, for all possible combinations of the case numbers and mesh refinement levels

2.5 Results

We present the results from different scenarios of computing the Taylor–Couette flow. We organize the results under the categories implied by the headings of Sects. 2.5.1–2.5.5.

Except for the computations reported in Sect. 2.5.5, we use the conservative form of the ST-VMS and the strong enforcement of the prescribed velocities.

In flow visualization, we use the configurations shown in Fig. 7. In the volume visualization, we show the isosurfaces of $u_z$. In the section visualization, we show $\left\| {\mathbf {u}}\right\|$ on the sections indicated in Fig. 7. In computations with rotational-periodicity representation of the flow field, we replicate the solution to form the full flow field. In computations with the NRF, we use the velocity representation in the IRF.

Fig. 7

Volume and section configurations used in the flow visualization. In the volume visualization, we remove $\frac{1}{4}$ of the flow field and show the results for the remaining part (gray). In the section visualization, we use the red and green sections placed at the midpoint along the radial and axial directions

We report two kinds of discrete Fourier decomposition of $u_\theta$. One is based on the temporal frequency, f, and the other based on the spatial wavelength $\lambda$. The Fourier coefficients are represented by $(u_\theta )_f$ and $(u_\theta )_\lambda$. We use the notation ${\mathcal {T}} = (T_1, T_2)$ to indicate the time range used in the Fourier decomposition and in time-averaging. An overbar, such as in ${\overline{u}}$, will indicate averaging in space or time, as specified where it is used, and double overbar, such as in $\overline{\overline{u}}$, will indicate averaging in space and time.

2.5.1 Flow development: Cases 1 and 4

We use the full-domain representation and the coarse mesh. Figures 8 and 9 show, for Case 1, for both the IRF and NRF, the isosurfaces of $u_{z}/\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ at various instants and the f-based Fourier decomposition of $u_\theta /\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ at $r = \frac{r_\mathrm {o}+r_\mathrm {i}}{2}$, $\theta = 0$. In the axial direction, we average the amplitude of the Fourier coefficients, and that diminishes the influence of the wave location in that direction. From both figures, we see that a time-periodic solution is reached after 50T. We see four waves in the circumferential direction. The waves move in the circumferential direction with the same sign as $\omega _\mathrm {i} - \omega _\mathrm {o}$, and the lower-mode solutions clearly reflect the motion of the waves.

Fig. 8

Case 1. Isosurfaces of $u_z/\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$. There are 16 isosurfaces, with values equally spaced from $-0.15$ to 0.15. IRF (left) and NRF (right) at $t=10T$, 30T, 50T and 70T (top to bottom)

Fig. 9

Case 1. f-based Fourier decomposition of $u_\theta /\left( U_\mathrm {i} - \eta U_\mathrm {o}\right)$ at $r =\frac{r_\mathrm {o}+r_\mathrm {i}}{2}$, $\theta = 0$, in the time ranges indicated in the legend, all spanning 10T. The amplitude is averaged over the axial direction using 56 equally-spaced points, numerous enough to diminish the influence of the wave location in that direction

Figures 10 and 11 show, for Case 4, velocity at $t=9T$ and angular velocity $\omega = u_\theta /r$ at various instants. We note that there is no variation in the axial direction in this case. With the initial condition being zero everywhere except the outer surface, the solution reaches the steady state after $t=3T$.

Fig. 10

Case 4. $\left\| {\mathbf {u}}\right\| /\left( U_\mathrm {i} - \eta U_\mathrm {o}\right)$ at $t=9T$

Fig. 11

Case 4. $\left( \omega _\mathrm {i} - {\overline{\omega }}\right) /\left( \omega _\mathrm {i} - \omega _\mathrm {o}\right)$ at instants indicated in the legend. The space averaging is in both axial and circumferential directions. The exact, steady-state solution is also shown

2.5.2 Rotational periodicity: Case 1

We use the coarse mesh. Figures 12 and 13 show, for both the IRF and NRF, for both the rotational-periodicity and full-domain representations of the flow field, the isosurfaces of $u_z/\left( U_\mathrm {i} - \eta U_\mathrm {o}\right)$ after the solutions become time-periodic and the $\lambda$-based Fourier decomposition of $u_\theta /\left( U_\mathrm {i} - \eta U_\mathrm {o}\right)$ at $t=80T$, $r = \frac{r_\mathrm {o}+r_\mathrm {i}}{2}$. In the axial direction, we average the amplitude of the Fourier coefficients. It can be clearly seen that the largest amplitude is at wavelength $\lambda = \pi \frac{r_\mathrm {o} + r_\mathrm {i}}{4}$. As expected, the results from the rotational-periodicity representation with the patch that has one-third-domain size do not match the results from the full-domain representation. In addition, the solutions for the IRF and NRF are somewhat different at higher modes, as can be seen in Fig. 13. However, with the patch that has half-domain size, the solution is in good agreement with what we get from the full-domain representation. Based on that, in the rest of the studies, we will use the rotational-periodicity representation with the patch that has half-domain size.

Fig. 12

Case 1. Isosurfaces of $u_z/\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ after the solutions become time-periodic, for the IRF (left) and NRF (right), for the rotational-periodicity (with patches that have quarter, one-third and half-domain sizes) and full-domain representations of the flow field (top to bottom). There are 16 isosurfaces, with values equally spaced from $-0.15$ to 0.15

Fig. 13

Case 1. $\lambda$-based Fourier decomposition of $u_\theta /\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ at $t=80T$, $r=\frac{r_\mathrm {o}+r_\mathrm {i}}{2}$, for the rotational-periodicity (with patches that have quarter, one-third and half-domain sizes) and full-domain representations of the flow field. The amplitude is averaged over the axial direction using 56 equally-spaced points, numerous enough to diminish the influence of the wave location in that direction

2.5.3 Mesh refinement: Case 1

We use the rotational-periodicity representation with the patch that has half-domain size. Figures 14 and 15 show, after the solutions become time-periodic, for both the IRF and NRF, for the three meshes, the isosurfaces of $u_z/\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ and section values of $\left\| {\mathbf {u}}\right\| /\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$. Figure 16 shows, for both the IRF and NRF, for the three meshes, $(\omega _\mathrm {i} - \overline{\overline{\omega }})/\left( \omega _\mathrm {i} - \omega _\mathrm {o}\right)$. The convergence with the mesh refinement is evident.

Fig. 14

Case 1. Isosurfaces of $u_z/\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ after the solutions become time-periodic. There are 16 isosurfaces, with values equally spaced from $-0.15$ to 0.15. IRF (left) and NRF (right). Coarse, medium and fine meshes (top to bottom)

Fig. 15

Case 1. $\left\| {\mathbf {u}}\right\| /\left( U_\mathrm {i} - \eta U_\mathrm {o}\right)$ after the solutions become time-periodic, on the sections indicated in Fig. 7. IRF (left) and NRF (right). Coarse, medium and fine meshes (top to bottom)

Fig. 16

Case 1. $(\omega _\mathrm {i} - \overline{\overline{\omega }} )/\left( \omega _\mathrm {i} - \omega _\mathrm {o}\right)$ over the range ${\mathcal {T}} = (97T, 100 T)$. The space averaging is in both axial and circumferential directions

Fig. 17

Case 2. Isosurfaces of $u_z/\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$. There are 16 isosurfaces, with values equally spaced from $-0.15$ to 0.15. IRF (left) and NRF (right) at $t=5T$, 10T, 15T and 20T (top to bottom)

Fig. 18

Case 2. $\left\| {\mathbf {u}}\right\| /\left( U_\mathrm {i} - \eta U_\mathrm {o}\right)$ on the sections indicated in Fig. 7. IRF (left) and NRF (right) at $t=5T$, 10T, 15T and 20T (top to bottom)

2.5.4 Flow patterns: Cases 2 and 3

We know from Fig. 2 that Case 3 has no waves in the circumferential direction. We also know from a prior computation we conducted for Case 2 with the full-domain representation, which we do not report here, that the number of waves in the circumferential direction is even. Therefore we compute both cases with the rotational-periodicity representation with the patch that has half-domain size. We use the medium mesh.

Figures 17 and 18 show, for Case 2, for both the IRF and NRF, the isosurfaces of $u_{z}/\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ and section values of $\left\| {\mathbf {u}}\right\| /\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ at various instants. For both the IRF and NRF, we see eight waves in the circumferential direction. They move in the circumferential direction with the same sign as $\omega _\mathrm {i}- \omega _\mathrm {o}$. Figure 19 shows, for both the IRF and NRF, $(\omega _\mathrm {i} - \overline{\overline{\omega }})/\left( \omega _\mathrm {i} - \omega _\mathrm {o}\right)$. It clearly indicates that the IRF and NRF solutions are in very good agreement.

Fig. 19

Case 2. $(\omega _\mathrm {i} - \overline{\overline{\omega }} )/\left( \omega _\mathrm {i} - \omega _\mathrm {o}\right)$ over the range ${\mathcal {T}} = (15T, 20 T)$. The space averaging is in both axial and circumferential directions

Figure 20 shows, for Case 3, for both the IRF and NRF, the isosurfaces of $u_{z}/\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ at $t = 15T$. For both the IRF and NRF, we see the Taylor vortex flow pattern. Figure 21 shows, for both the IRF and NRF, $(\omega _\mathrm {i} - \overline{\overline{\omega }})/\left( \omega _\mathrm {i} - \omega _\mathrm {o}\right)$, which clearly indicates that the two solutions are in very good agreement.

Fig. 20

Case 3. Isosurfaces of $u_z/\left( U_\mathrm {i}-\eta U_\mathrm {o}\right)$ at $t=15T$. There are 16 isosurfaces, with values equally spaced from $-0.15$ to 0.15. IRF and NRF

Fig. 21

Case 3. $(\omega _\mathrm {i} - \overline{\overline{\omega }} )/\left( \omega _\mathrm {i} - \omega _\mathrm {o}\right)$ over the range ${\mathcal {T}} = (10T, 15 T)$. The space averaging is in both axial and circumferential directions

2.5.5 Methods: Case 1

There are total eight test cases, which are all combinations of the conservative and convective forms of the ST-VMS, strong and weak enforcement of the prescribed velocities, and the IRF and NRF choices of the reference frame. We use the rotational-periodicity representation with the patch that has half-domain size and the medium mesh. For detailed comparison purposes, the number of nonlinear iterations per time step is increased to 5. The flow patterns are almost the same in all test cases, therefore we compare $(\omega _\mathrm {i} - \overline{\overline{\omega }})/\left( \omega _\mathrm {i} - \omega _\mathrm {o}\right)$ in Fig. 22. The solutions are very close.

Fig. 22

Case 1. $(\omega _\mathrm {i} - \overline{\overline{\omega }} )/\left( \omega _\mathrm {i} - \omega _\mathrm {o}\right)$ over the range ${\mathcal {T}} = (105T, 110 T)$. The space averaging is in both axial and circumferential directions. We also show, as a reference, the solution computed over the fine mesh, using the conservative form of the ST-VMS, strong enforcement of the prescribed velocities, and the IRF. For that, as can be deduced from Table 2, we use the rotational-periodicity representation with the patch that has half-domain size

We also check the global conservation of the angular momentum and the torque acting on the inner boundary. We define the angular momentum at time level n as

$$\begin{aligned} {\mathbf {L}}_{n}^-&= \int _{\varOmega _n} {\mathbf {x}} \times {\mathbf {u}}_n^- \mathrm {d}\varOmega , \end{aligned}$$

(3)

and the angular momentum time-averaged over the range ${\mathcal {T}} =(T_1, T_2)$ as

$$\begin{aligned} {\overline{{\mathbf {L}}}}&= \frac{1}{T_2 - T_1} \int _{Q} {\mathbf {x}} \times {\mathbf {u}}\mathrm {d}Q, \end{aligned}$$

(4)

where Q is the ST domain between the instants $T_1$ and $T_2$, and the definition of $(~)^-_{n}$ is as given in Appendix A. We use the symbols $L_n^-$ and ${\overline{L}}$ to denote the axial components. We define the torque acting on the inner and outer boundaries as

$$\begin{aligned} (\pmb {\varUpsilon }_{n})_\mathrm {IB}&= \frac{1}{\varDelta t} \int _{(P_n)_\mathrm {IB}} {\mathbf {x}} \times {\mathbf {h}} \mathrm {d}P, \end{aligned}$$

(5)

$$\begin{aligned} (\pmb {\varUpsilon }_{n})_\mathrm {OB}&= \frac{1}{\varDelta t} \int _{(P_n)_\mathrm {OB}} {\mathbf {x}} \times {\mathbf {h}} \mathrm {d}P, \end{aligned}$$

(6)

with the definition of $P_n$ as given in Appendix A. We note that ${\mathbf {h}}$ (defined in Appendix A) is computed on the prescribed-velocity boundary as part of the flow solution. Therefore the way it is computed is different between the strong and weak enforcement of the velocity boundary condition. We use the symbols $(\varUpsilon _n)_\mathrm {IB}$ and $(\varUpsilon _n)_\mathrm {OB}$ to denote the axial components. The global conservation of the angular momentum can be written as

$$\begin{aligned} \frac{ {\mathbf {L}}_{n+1}^- - {\mathbf {L}}_{n}^- }{\varDelta t}&= (\pmb {\varUpsilon }_{n})_\mathrm {IB} + (\pmb {\varUpsilon }_{n})_\mathrm {OB} \end{aligned}$$

(7)

for every time step. Because of the stabilization terms, only the conservative form of the ST-VMS satisfies Eq. (7) (see Remark 5).

Table 4 shows ${\overline{L}}/{\overline{L}}_\mathrm {FINE}$ and ${\overline{\varUpsilon }}_\mathrm {IB} T/{\overline{L}}_\mathrm {FINE}$, for all eight test cases, over the range ${\mathcal {T}} = (110T, 115T)$. The time-averaged torque is defined using the torque counterpart of Eq. (4). We obtain ${\overline{L}}_\mathrm {FINE}$ from the solution computed over the fine mesh, using the conservative form of the ST-VMS, strong enforcement of the prescribed velocities, and the IRF, and by time-averaging over the range ${\mathcal {T}} = (97 T, 102 T)$. All ${\overline{L}}$ values are very close to ${\overline{L}}_\mathrm {FINE}$, and the ${\overline{\varUpsilon }}_\mathrm {IB}$ values vary slightly between the cases. Figures 23, 24 and 25 show the global conservation of the angular momentum. For clear comparisons with amplified values, we use the expressions $\frac{(L^-_{n+1} - L^-_{n})T}{\varDelta t {\overline{L}}_\mathrm {FINE}}$, $\frac{(\varUpsilon _\mathrm {IB} + \varUpsilon _\mathrm {OB}) T}{{\overline{L}}_\mathrm {FINE}}$. We see that all the solutions practically have global conservation at every time step. In all test cases, there are fluctuations with period T, and the fluctuation magnitude varies between the cases. The prescribed velocities are represented exactly in the computations with the NRF and in the computations with weak enforcement (see Remark 7). Therefore, the fluctuations are not coming from how the prescribed velocities are represented but possibly from the inertia itself or nonuniform element lengths in the circumferential direction (see Fig. 5). Even in the computation with convective form of the ST-VMS, weak enforcement of the prescribed velocities, and the NRF, where we see the largest fluctuations in $\frac{\left( L_{n+1}^- - L_n^-\right) T}{\varDelta T {\overline{L}}_\mathrm {FINE}}$, the maximum fluctuation is $0.3~\%$. Therefore, we consider these fluctuations to be acceptable in computation of the wavy vortex flows.

3 Concluding remarks

We have shown how the ST-VMS with ST isogeometric discretization can be used to address the computational challenges encountered in the Taylor–Couette flow, a classical fluid mechanics problem. The Taylor–Couette flow exhibits, depending on the Reynolds number, a range of flow patterns, with the interesting ones having small-scale structures, and sometimes even wavy nature. Accurate representation of these flow patterns in computational flow analysis requires methods that can, with a reasonable computational cost, represent the circular geometry accurately and provide a high-fidelity flow solution. The ST framework provides higher-order accuracy in general, and the VMS feature of the ST-VMS addresses the computational challenges associated with the multiscale nature of the flow. The ST isogeometric discretization enables exact representation of the circular geometry and increased accuracy in the flow solution. In computations where the reference-frame choice leads to rotating the mesh, the STNMUM, with NURBS basis functions in time, enables exact representation of the mesh rotation, in terms of both the paths of the mesh points and the velocity of the points along their paths. In computations with rotational-periodicity representation of the flow field, the periodicity is enforced with the ST-SI.

We conducted the computational analysis with different combinations of the Reynolds numbers based on the inner and outer cylinder rotation speeds, with different choices of the reference frame, one of which leads to rotating the mesh, with the full-domain and rotational-periodicity representations of the flow field, with both the convective and conservative forms of the ST-VMS, with both the strong and weak enforcement of the prescribed velocities on the cylinder surfaces, and with different mesh refinements. With the combinations of the Reynolds numbers used in the computations, we covered the cases leading to the Taylor vortex flow and the wavy vortex flow, where the waves are in motion. We showed that the ST methods described in the article, integrated together, offer a high-fidelity computational analysis platform for the Taylor–Couette flow and for other classes of flow problems with similar features.

Table 4 Case 1. ${\overline{L}}/{\overline{L}}_\mathrm {FINE}$ and ${\overline{\varUpsilon }}_\mathrm {IB} T/{\overline{L}}_\mathrm {FINE}$ for the eight test cases

Fig. 23

Case 1. Global conservation of the angular momentum for the IRF

Fig. 24

Case 1. Global conservation of the angular momentum for the NRF

Fig. 25

Case 1. Global conservation of the angular momentum for the solution computed over the fine mesh, using the conservative form of the ST-VMS, strong enforcement of the prescribed velocities, and the IRF and NRF

Appendices

ST-VMS and ST-SI

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

1.1 ST-VMS

The convective form of the ST-VMS is given as

$$\begin{aligned}&\int _{Q_n} {\mathbf {w}}^h \cdot \rho \left( \frac{\partial {\mathbf {u}}^h}{\partial t} + {\mathbf {u}}^h \cdot \pmb {\nabla }{\mathbf {u}}^h - {\mathbf {f}}^h \right) \mathrm {d}Q \nonumber \\&\quad + \int _{Q_n} \pmb {\varepsilon }({\mathbf {w}}^h) : \pmb {\sigma }^h \mathrm {d}Q-\int _{\left( P_n\right) _\mathrm {h}} {\mathbf {w}}^h \cdot {{\mathbf {h}}}^h \mathrm {d}P\nonumber \\&\quad +\int _{Q_n} q^h \pmb {\nabla }\cdot {\mathbf {u}}^h \mathrm {d}Q+ \int _{\varOmega _n} ({\mathbf {w}}^h)_{n}^{+} \cdot \rho \left( ({\mathbf {u}}^h)_{n}^+ - ({\mathbf {u}}^h)_{n}^- \right) \mathrm {d}\varOmega \nonumber \\&\quad + \sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \frac{\tau _{\mathrm {SUPS}}}{\rho } \left( \rho \left( \frac{\partial {\mathbf {w}}^h}{\partial t} + {\mathbf {u}}^h \cdot \pmb {\nabla }{\mathbf {w}}^h \right) +\pmb {\nabla }q^h \right) \cdot {\mathbf {r}}_{\mathrm {M}}^h \mathrm {d}Q\nonumber \\&\quad + \sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \nu _{\mathrm {LSIC}}\pmb {\nabla }\cdot {\mathbf {w}}^h \rho {r}_{\mathrm {C}}^h \mathrm {d}Q\nonumber \\&\quad -\sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \tau _{\mathrm {SUPS}}{\mathbf {w}}^h \cdot \left( {\mathbf {r}}_{\mathrm {M}}^h \cdot \pmb {\nabla }{\mathbf {u}}^h \right) \mathrm {d}Q\nonumber \\&\quad -\sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \frac{\tau _{\mathrm {SUPS}}^2}{\rho } {\mathbf {r}}_{\mathrm {M}}^h \cdot \left( \pmb {\nabla }{\mathbf {w}}^h\right) \cdot {\mathbf {r}}_{\mathrm {M}}^h \mathrm {d}Q\nonumber \\&= 0, \end{aligned}$$

(8)

where

$$\begin{aligned} {\mathbf {r}}_{\mathrm {M}}^h&= \rho \left( \frac{\partial {\mathbf {u}}^h}{\partial t} + {\mathbf {u}}^h \cdot \pmb {\nabla }{\mathbf {u}}^h - {\mathbf {f}}^h \right) - \pmb {\nabla }\cdot \pmb {\sigma }^h, \end{aligned}$$

(9)

$$\begin{aligned} {r}_{\mathrm {C}}^h&= \pmb {\nabla }\cdot {\mathbf {u}}^h \end{aligned}$$

(10)

are the residuals of the momentum equation and incompressibility constraint. Here, $p$, ${\mathbf {f}}$, and ${\mathbf {h}}$ are the pressure, body force, and the stress specified at the boundary. In the case of an NRF with angular velocity $\pmb {\omega }$, ${\mathbf {f}}^h$ with gravity ${\mathbf {a}}_\mathrm {GRAV}$ becomes

$$\begin{aligned} {\mathbf {f}}^h = {\mathbf {a}}_\mathrm {GRAV}- \dot{\pmb {\omega }}\times {\mathbf {r}}^h- 2 \pmb {\omega } \times {\mathbf {u}}^h - \pmb {\omega } \times \left( \pmb {\omega } \times {\mathbf {r}}^h\right) , \end{aligned}$$

(11)

where ${\mathbf {r}}^h = {\mathbf {x}}^h - {\mathbf {x}}_0$ is the position vector from the rotation center ${\mathbf {x}}_0$. The stress tensor is defined as $\pmb {\sigma }= - p{\mathbf {I}}+ 2 \mu \pmb {\varepsilon }({\mathbf {u}})$, where ${\mathbf {I}}$ is the identity tensor and $\pmb {\varepsilon }\left( {\mathbf {u}}\right) = \left( \left( \pmb {\nabla }{\mathbf {u}}\right) +\left( \pmb {\nabla }{\mathbf {u}}\right) ^{\mathop \intercal }\right) /2$ is the strain-rate tensor. The test functions associated with the ${\mathbf {u}}$ and $p$ are ${\mathbf {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$, $\left( P_n\right) _\mathrm {h}$ is the part of the slice lateral boundary associated with the boundary condition ${\mathbf {h}}$, and $\varOmega _n$ is the spatial domain at time level n. The superscript “e” is the ST element counter, and $n_{\mathrm {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 4

The ST-SUPS can be obtained from the ST-VMS by dropping the eighth and ninth integrations.

The stabilization parameters, $\tau _{\mathrm {SUPS}}$ and $\nu _{\mathrm {LSIC}}$, will be given in Appendix B.1.

The conservative form of the ST-VMS [1, 93] is given as

$$\begin{aligned}&\int _{Q_n} {\mathbf {w}}^h \cdot \rho \left( \frac{\partial {\mathbf {u}}^h}{\partial t} + \pmb {\nabla }\cdot \left( {\mathbf {u}}^h {\mathbf {u}}^h \right) - {\mathbf {f}}^h \right) \mathrm {d}Q \nonumber \\&\quad + \int _{Q_n} \pmb {\varepsilon }({\mathbf {w}}^h) : \pmb {\sigma }^h \mathrm {d}Q-\int _{\left( P_n\right) _\mathrm {h}} {\mathbf {w}}^h \cdot {{\mathbf {h}}}^h \mathrm {d}P\nonumber \\&\quad +\int _{Q_n} q^h \pmb {\nabla }\cdot {\mathbf {u}}^h \mathrm {d}Q+ \int _{\varOmega _n} ({\mathbf {w}}^h)_{n}^{+} \cdot \rho \left( ({\mathbf {u}}^h)_{n}^+ - ({\mathbf {u}}^h)_{n}^- \right) \mathrm {d}\varOmega \nonumber \\&\quad + \sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \frac{\tau _{\mathrm {SUPS}}}{\rho } \left( \rho \left( \frac{\partial {\mathbf {w}}^h}{\partial t} + {\mathbf {u}}^h \cdot \pmb {\nabla }{\mathbf {w}}^h \right) +\pmb {\nabla }q^h \right) \cdot {\mathbf {r}}_{\mathrm {M}}^h \mathrm {d}Q\nonumber \\&\quad + \sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \nu _{\mathrm {LSIC}}\pmb {\nabla }\cdot {\mathbf {w}}^h \rho {r}_{\mathrm {C}}^h \mathrm {d}Q\nonumber \\&\quad +\sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \tau _{\mathrm {SUPS}}{\mathbf {r}}_{\mathrm {M}}^h \cdot \pmb {\nabla }{\mathbf {w}}^h \cdot {\mathbf {u}}^h \mathrm {d}Q\nonumber \\&\quad -\sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \frac{\tau _{\mathrm {SUPS}}^2}{\rho } {\mathbf {r}}_{\mathrm {M}}^h \cdot \left( \pmb {\nabla }{\mathbf {w}}^h\right) \cdot {\mathbf {r}}_{\mathrm {M}}^h \mathrm {d}Q\nonumber \\&= 0. \end{aligned}$$

(12)

Remark 5

Combining the SUPG term in the sixth integration and the eighth integration gives the term

$$\begin{aligned}&\sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \tau _{\mathrm {SUPS}}\left( {\mathbf {u}}^h \cdot \pmb {\nabla }{\mathbf {w}}^h \cdot {\mathbf {r}}_{\mathrm {M}}^h + {\mathbf {r}}_{\mathrm {M}}^h \cdot \pmb {\nabla }{\mathbf {w}}^h \cdot {\mathbf {u}}^h \right) \mathrm {d}Q\nonumber \\&\quad = \sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \tau _{\mathrm {SUPS}}\pmb {\varepsilon }({\mathbf {w}}^h) : \left( {\mathbf {u}}^h {\mathbf {r}}_{\mathrm {M}}^h + {\mathbf {r}}_{\mathrm {M}}^h {\mathbf {u}}^h \right) \mathrm {d}Q, \end{aligned}$$

(13)

and the tenth integration can be rearranged as

$$\begin{aligned}&-\sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \frac{\tau _{\mathrm {SUPS}}^2}{\rho } {\mathbf {r}}_{\mathrm {M}}^h \cdot \left( \pmb {\nabla }{\mathbf {w}}^h\right) \cdot {\mathbf {r}}_{\mathrm {M}}^h \mathrm {d}Q\nonumber \\&\quad = -\sum _{e=1}^{(n_{\mathrm {el}})_n} \int _{Q_n^e} \frac{\tau _{\mathrm {SUPS}}^2}{\rho } \pmb {\varepsilon }({\mathbf {w}}^h) : {\mathbf {r}}_{\mathrm {M}}^h {\mathbf {r}}_{\mathrm {M}}^h \mathrm {d}Q. \end{aligned}$$

(14)

These symmetric terms do not produce angular momentum (see [9]). This is not the case for the convective form of the ST-VMS.

1.2 ST-SI

1.2.1 Two-side formulation (fluid–fluid 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 Eqs. (8) or (12) include added boundary terms corresponding to the SI. When the periodicity is enforced with an SI, the unit normal vectors ${{\mathbf {n}}}_\mathrm {A}$ and ${{\mathbf {n}}}_\mathrm {B}$ appearing in the SI terms are the normal vectors associated with the periodicity surfaces, and we define ${{\mathbf {t}}}_\mathrm {A}$ and ${{\mathbf {t}}}_\mathrm {B}$ as unit vectors tangent to those surfaces. Then we define ${\mathbf {R}}$ as the tensor transforming the basis set $({{\mathbf {n}}}_\mathrm {A},\, {{\mathbf {t}}}_\mathrm {A},\, {{\mathbf {n}}}_\mathrm {A}\times {{\mathbf {t}}}_\mathrm {A})$ to $(-{{\mathbf {n}}}_\mathrm {B},\, {{\mathbf {t}}}_\mathrm {B},\, -{{\mathbf {n}}}_\mathrm {B}\times {{\mathbf {t}}}_\mathrm {B})$, which is expressed as

$$\begin{aligned} {\mathbf {R}}&= - {{\mathbf {n}}}_\mathrm {B} {{\mathbf {n}}}_\mathrm {A} + {{\mathbf {t}}}_\mathrm {B} {{\mathbf {t}}}_\mathrm {A} - \left( {{\mathbf {n}}}_\mathrm {B} \times {{\mathbf {t}}}_\mathrm {B} \right) \left( {{\mathbf {n}}}_\mathrm {A} \times {{\mathbf {t}}}_\mathrm {A} \right) . \end{aligned}$$

(15)

With that, we now describe the ST-SI in its form generalized to enforce periodicity. The boundary terms for the two sides are first added separately, using the test functions ${{\mathbf {w}}}_\mathrm {A}^h$ and ${q}_\mathrm {A}^h$ and ${{\mathbf {w}}}_\mathrm {B}^h$ and ${q}_\mathrm {B}^h$. Then, putting together the terms added to each side, the complete set of terms added becomes

$$\begin{aligned}&- \int _{\left( P_n\right) _\mathrm {SI}} \left( {q}_\mathrm {B}^h {{\mathbf {n}}}_\mathrm {B} - {\mathbf {R}}\cdot {q}_\mathrm {A}^h {{\mathbf {n}}}_\mathrm {A} \right) \cdot \frac{1}{2} \left( {{\mathbf {u}}}_\mathrm {B}^h - {\mathbf {R}}\cdot {{\mathbf {u}}}_\mathrm {A}^h \right) \mathrm {d}P \nonumber \\&\quad - \int _{\left( P_n\right) _\mathrm {SI}} \rho {{\mathbf {w}}}_\mathrm {B}^h \cdot \frac{1}{2} \left( \left( {{\mathcal {F}}}_\mathrm {B}^h - \left| {{\mathcal {F}}}_\mathrm {B}^h \right| \right) {{\mathbf {u}}}_\mathrm {B}^h \right. \nonumber \\&\left. \quad - \left( {{\mathcal {F}}}_\mathrm {B}^h - \left| {{\mathcal {F}}}_\mathrm {B}^h \right| \right) {\mathbf {R}}\cdot {{\mathbf {u}}}_\mathrm {A}^h \right) \mathrm {d}P \nonumber \\&\quad - \int _{\left( P_n\right) _\mathrm {SI}} \rho \left( {\mathbf {R}}\cdot {{\mathbf {w}}}_\mathrm {A}^h \right) \cdot \frac{1}{2} \left( \left( {{\mathcal {F}}}_\mathrm {A}^h - \left| {{\mathcal {F}}}_\mathrm {A}^h \right| \right) {\mathbf {R}}\cdot {{\mathbf {u}}}_\mathrm {A}^h \right. \nonumber \\&\left. \quad - \left( {{\mathcal {F}}}_\mathrm {A}^h - \left| {{\mathcal {F}}}_\mathrm {A}^h \right| \right) {{\mathbf {u}}}_\mathrm {B}^h \right) \mathrm {d}P \nonumber \\&\quad + \int _{\left( P_n\right) _\mathrm {SI}} \left( {{\mathbf {n}}}_\mathrm {B} \cdot {{\mathbf {w}}}_\mathrm {B}^h + {{\mathbf {n}}}_\mathrm {A} \cdot {{\mathbf {w}}}_\mathrm {A}^h \right) \frac{1}{2} \left( {p}_\mathrm {B}^h + {p}_\mathrm {A}^h \right) \mathrm {d}P \nonumber \\&\quad - \int _{\left( P_n\right) _\mathrm {SI}} \left( {{\mathbf {n}}}_\mathrm {B} \cdot {{\mathbf {w}}}_\mathrm {B}^h - {{\mathbf {n}}}_\mathrm {A} \cdot {{\mathbf {w}}}_\mathrm {A}^h \right) \frac{1}{2} \varDelta p\mathrm {d}P \nonumber \\&\quad - \int _{\left( P_n\right) _\mathrm {SI}} \left( {{\mathbf {w}}}_\mathrm {B}^h - {\mathbf {R}}\cdot {{\mathbf {w}}}_\mathrm {A}^h \right) \cdot \left( {{\hat{{\mathbf {n}}}}}_\mathrm {B} \cdot \mu \left( \pmb {\varepsilon }({{\mathbf {u}}}_\mathrm {B}^h) \right. \right. \nonumber \\&\quad \left. \left. + {\mathbf {R}}\cdot \pmb {\varepsilon }({{\mathbf {u}}}_\mathrm {A}^h) \cdot {\mathbf {R}}^{\mathop \intercal } \right) \right) \mathrm {d}P \nonumber \\&\quad - \gamma \int _{\left( P_n\right) _\mathrm {SI}} {{\hat{{\mathbf {n}}}}}_\mathrm {B} \cdot \mu \left( \pmb {\varepsilon }\left( {{\mathbf {w}}}_\mathrm {B}^h \right) + {\mathbf {R}}\cdot \pmb {\varepsilon }\left( {{\mathbf {w}}}_\mathrm {A}^h \right) \cdot {\mathbf {R}}^{\mathop \intercal } \right) \cdot \left( {{\mathbf {u}}}_\mathrm {B}^h \right. \nonumber \\&\quad \left. - {\mathbf {R}}\cdot {{\mathbf {u}}}_\mathrm {A}^h \right) \mathrm {d}P \nonumber \\&\quad + \int _{\left( P_n\right) _\mathrm {SI}} \frac{\mu C}{h} \left( {{\mathbf {w}}}_\mathrm {B}^h - {\mathbf {R}}\cdot {{\mathbf {w}}}_\mathrm {A}^h \right) \cdot \left( {{\mathbf {u}}}_\mathrm {B}^h - {\mathbf {R}}\cdot {{\mathbf {u}}}_\mathrm {A}^h \right) \mathrm {d}P, \end{aligned}$$

(16)

where

$$\begin{aligned} {{\mathcal {F}}}_\mathrm {B}^h&= {{\mathbf {n}}}_\mathrm {B} \cdot \left( {{\mathbf {u}}}_\mathrm {B}^h - {{\mathbf {v}}}_\mathrm {B}^h \right) , \end{aligned}$$

(17)

$$\begin{aligned} {{\mathcal {F}}}_\mathrm {A}^h&= {{\mathbf {n}}}_\mathrm {A} \cdot \left( {{\mathbf {u}}}_\mathrm {A}^h - {{\mathbf {v}}}_\mathrm {A}^h \right) , \end{aligned}$$

(18)

$$\begin{aligned} {{\hat{{\mathbf {n}}}}}_\mathrm {B}&= \frac{ {{\mathbf {n}}}_\mathrm {B} - {\mathbf {R}}\cdot {{\mathbf {n}}}_\mathrm {A} }{ \left\| {{\mathbf {n}}}_\mathrm {B} - {\mathbf {R}}\cdot {{\mathbf {n}}}_\mathrm {A} \right\| } . \end{aligned}$$

(19)

Here, $\left( P_n\right) _\mathrm {SI}$ is the SI in the ST domain, ${\mathbf {n}}$ is the unit normal vector, ${\mathbf {v}}$ is the mesh velocity, $\gamma =1$, and C is a nondimensional constant. The element length h will be defined in Appendix B.2. In Remark 3 of [133], it was proposed to use ST-SI also as a way of imposing the periodicity in ST-VMS computations with spatial periodicity. As mentioned in the same remark, in such cases, an SI is placed where the periodicity is imposed and the computations are carried out while including the corresponding ST-SI term. The fifth term in Eq. (16), with $\varDelta p$ representing the pressure jump across the SI, is the ST-SI term corresponding to the periodicity. We note that, consistent with the notation in [6, 103,104,105], positive $\varDelta p$ means higher pressure on Side A.

Remark 6

When the SI is a standard fluid–fluid SI, ${{\mathbf {n}}}_\mathrm {B} = -{{\mathbf {n}}}_\mathrm {A}$, ${{\mathbf {t}}}_\mathrm {B} = {{\mathbf {t}}}_\mathrm {A}$, and ${\mathbf {R}}$ becomes the identity tensor.

1.2.2 One-side formulation (fluid–solid SI)

On solid surfaces where we prefer weak enforcement of the Dirichlet conditions [134, 135] for the fluid, we use the ST-SI version where the SI is between the fluid and solid domains. This version is obtained (see [6]) by starting with the terms added to Side B and replacing the Side A velocity with the velocity ${\mathbf {g}}^h$ coming from the solid domain. Then the SI terms added to Eq. (8) to represent the weakly-enforced Dirichlet conditions become

$$\begin{aligned}&- \int _{\left( P_n\right) _\mathrm {SI}} {q}_\mathrm {B}^h {{\mathbf {n}}}_\mathrm {B} \cdot {{\mathbf {u}}}_\mathrm {B}^h \mathrm {d}P - \int _{\left( P_n\right) _\mathrm {SI}} \rho {{\mathbf {w}}}_\mathrm {B}^h \cdot {{\mathcal {F}}}_\mathrm {B}^h {{\mathbf {u}}}_\mathrm {B}^h \mathrm {d}P \nonumber \\&\quad + \int _{\left( P_n\right) _\mathrm {SI}} {q}_\mathrm {B}^h {{\mathbf {n}}}_\mathrm {B} \cdot {\mathbf {g}}^h \mathrm {d}P \nonumber \\&\quad + \int _{\left( P_n\right) _\mathrm {SI}} \rho {{\mathbf {w}}}_\mathrm {B}^h \cdot \frac{1}{2} \left( \left( {{\mathcal {F}}}_\mathrm {B}^h + \left| {{\mathcal {F}}}_\mathrm {B}^h \right| \right) {{\mathbf {u}}}_\mathrm {B}^h + \left( {{\mathcal {F}}}_\mathrm {B}^h - \left| {{\mathcal {F}}}_\mathrm {B}^h \right| \right) {\mathbf {g}}^h \right) \mathrm {d}P \nonumber \\&\quad - \int _{\left( P_n\right) _\mathrm {SI}} {{\mathbf {w}}}_\mathrm {B}^h \cdot \left( {{\mathbf {n}}}_\mathrm {B} \cdot {\pmb {\sigma }}_\mathrm {B}^h \right) \mathrm {d}P \nonumber \\&\quad - \gamma \int _{\left( P_n\right) _\mathrm {SI}} {{\mathbf {n}}}_\mathrm {B} \cdot 2\mu \pmb {\varepsilon }\left( {{\mathbf {w}}}_\mathrm {B}^h \right) \cdot \left( {{\mathbf {u}}}_\mathrm {B}^h - {\mathbf {g}}^h \right) \mathrm {d}P \nonumber \\&\quad + \int _{\left( P_n\right) _\mathrm {SI}} \frac{\mu C}{{h}_\mathrm {B}} {{\mathbf {w}}}_\mathrm {B}^h \cdot \left( {{\mathbf {u}}}_\mathrm {B}^h - {\mathbf {g}}^h \right) \mathrm {d}P. \end{aligned}$$

(20)

The element length $h_\mathrm {B}$ will be given in Appendix B.2.

Stabilization parameters and element lengths

1.1 ST-VMS

There are various ways of defining the stabilization parameters $\tau _{\mathrm {SUPS}}$ and $\nu _{\mathrm {LSIC}}$. Here, $\tau _{\mathrm {SUPS}}$ is mostly from [128]:

$$\begin{aligned} \tau _{\mathrm {SUPS}}&= \left( \tau _{\mathrm {SUGN12}}^{-2} + \tau _{\mathrm {SUGN3}}^{-2} + \tau _{\mathrm {SUGN4}}^{-2}\right) ^{-\frac{1}{2}}. \end{aligned}$$

(21)

The first and second components are given as

$$\begin{aligned} \tau _{\mathrm {SUGN12}}^{-2}&= \begin{bmatrix} 1 \\ {\mathbf {u}}\\ \end{bmatrix} \begin{bmatrix} 1 \\ {\mathbf {u}}\\ \end{bmatrix} : {\mathbf {G}}^\mathrm {ST} \end{aligned}$$

(22)

and

$$\begin{aligned} \tau _{\mathrm {SUGN3}}^{-1}&= \nu {\mathbf {r}}{\mathbf {r}}:{\mathbf {G}}, \end{aligned}$$

(23)

where ${\mathbf {r}}$ is the solution-gradient direction:

$$\begin{aligned} {\mathbf {r}}&= \frac{\pmb {\nabla }\left\| {\mathbf {u}}\right\| }{\left\| \pmb {\nabla }\left\| {\mathbf {u}}\right\| \right\| }. \end{aligned}$$

(24)

Here ${\mathbf {G}}^\mathrm {ST}$ and ${\mathbf {G}}$ are the ST and space-only element metric tensors:

$$\begin{aligned} {\mathbf {G}}^\mathrm {ST}&= \left( \hat{{\mathbf {Q}}}^\mathrm {ST}\right) ^{-T} \cdot \left( \hat{{\mathbf {Q}}}^\mathrm {ST}\right) ^{-1}, \end{aligned}$$

(25)

$$\begin{aligned} {\mathbf {G}}&= \hat{{\mathbf {Q}}}^{-T} \cdot \hat{{\mathbf {Q}}}^{-1}, \end{aligned}$$

(26)

where

$$\begin{aligned} \hat{{\mathbf {Q}}}^\mathrm {ST}&= {\mathbf {Q}}^\mathrm {ST} \cdot \left( {\mathbf {D}}^\mathrm {ST}\right) ^{-1}, \end{aligned}$$

(27)

$$\begin{aligned} \hat{{\mathbf {Q}}}&= {\mathbf {Q}} \cdot {\mathbf {D}}^{-1}. \end{aligned}$$

(28)

The ST and space-only Jacobian tensors are

$$\begin{aligned} {\mathbf {Q}}^\mathrm {ST}&= \begin{bmatrix} \frac{\partial t}{\partial \theta } &{}\quad \frac{\partial t}{\partial \pmb {\xi }} \\ \frac{\partial {\mathbf {x}}}{\partial \theta } &{}\quad {\mathbf {Q}} \\ \end{bmatrix} \end{aligned}$$

(29)

and

$$\begin{aligned} {\mathbf {Q}}&= \frac{\partial {\mathbf {x}}}{\partial \pmb {\xi }}, \end{aligned}$$

(30)

where $\theta$ and $\pmb {\xi }$ are the temporal and spatial parametric coordinates. The transformation tensor ${\mathbf {D}}^\mathrm {ST}$ is defined as

$$\begin{aligned} {\mathbf {D}}^\mathrm {ST}&= \begin{bmatrix} D_\theta &{}\quad {\mathbf {0}}^T \\ {\mathbf {0}}&{}\quad {\mathbf {D}} \\ \end{bmatrix}. \end{aligned}$$

(31)

The definitions used for $D_\theta$ and ${\mathbf {D}}$ play an important role, especially for higher-order isogeometric discretization [128, 130] and simplex elements [129]. In this article, we set $D_\theta =1$ and set ${\mathbf {D}}$ to its “RQD-MAX” version [130].

The third component, originating from [2], is defined as

$$\begin{aligned} \tau _{\mathrm {SUGN4}}&= \left\| \pmb {\nabla }{\mathbf {u}}^h\right\| _F^{-1} , \end{aligned}$$

(32)

where $\Vert \cdot \Vert _F$ represents the Frobenius norm.

The stabilization parameter $\nu _{\mathrm {LSIC}}$ is from [78]:

$$\begin{aligned} \nu _{\mathrm {LSIC}}&= \frac{h^2_\mathrm {LSIC}}{\tau _{\mathrm {SUPS}}}, \end{aligned}$$

(33)

where $h_\mathrm {LSIC}$ is set equal to the minimum element length $h_\mathrm {MIN}$:

$$\begin{aligned} h_\mathrm {MIN}&= 2 \left( \max _{{\mathbf {r}}} \left( {\mathbf {r}} {\mathbf {r}} : {\mathbf {G}} \right) \right) ^{-\frac{1}{2}} . \end{aligned}$$

(34)

1.2 ST-SI

The element length used in the ST-SI is given as

$$\begin{aligned} h&= \left( \frac{{h}_\mathrm {B}^{-2}+ {h}_\mathrm {A}^{-2}}{2}\right) ^{-\frac{1}{2}}, \end{aligned}$$

(35)

$$\begin{aligned} {h}_\mathrm {B}&= 2 \left( {{\mathbf {n}}}_\mathrm {B} {{\mathbf {n}}}_\mathrm {B} : {\mathbf {G}} \right) ^{-\frac{1}{2}}&\text {(for Side B)}, \end{aligned}$$

(36)

$$\begin{aligned} {h}_\mathrm {A}&= 2 \left( {{\mathbf {n}}}_\mathrm {A} {{\mathbf {n}}}_\mathrm {A} : {\mathbf {G}} \right) ^{-\frac{1}{2}}&\text {(for Side A)}, \end{aligned}$$

(37)

$$\begin{aligned} {{\hat{{\mathbf {n}}}}}_\mathrm {B}&= \frac{ {{\mathbf {n}}}_\mathrm {B} - {{\mathbf {n}}}_\mathrm {A} }{ \left\| {{\mathbf {n}}}_\mathrm {B} - {{\mathbf {n}}}_\mathrm {A} \right\| }. \end{aligned}$$

(38)

Circular-arc geometry, velocity and motion

1.1 Geometry

An arc with its angle less than $\pi$ can be represented exactly by a quadratic NURBS element. Using the half angle $0< q< \frac{\pi }{2}$, we place the control points as shown in Fig. 26.

Fig. 26

A circular arc represented by a quadratic NURBS element in space

The weights are $w_1 = w_3 = 1$, and $w_2 = \cos q$, where

$$\begin{aligned} \cos 2 q&= \frac{{{\mathbf {x}}}_1 \cdot {{\mathbf {x}}}_3}{r^2},&r&= \Vert {{\mathbf {x}}}_1\Vert = \Vert {{\mathbf {x}}}_3\Vert , \end{aligned}$$

(39)

and

$$\begin{aligned} \Vert {{\mathbf {x}}}_2\Vert = \frac{r}{\cos q}. \end{aligned}$$

(40)

With that, the circular arc can be represented with NURBS basis functions $N_a$:

$$\begin{aligned} {\mathbf {x}}(\xi )&= \sum _{a=1}^3 N_a (\xi ) {\mathbf {x}}_a. \end{aligned}$$

(41)

In the case of $q= \frac{\pi }{2}$, we can represent the circular arc exactly by using two elements. With a uniform knot space, and $w_1 = w_4 = 1$, and $w_2 = w_3 = \frac{1}{2}$, the control points are placed as shown in Fig. 27, and the control points satisfy the conditions

$$\begin{aligned} \Vert {{\mathbf {x}}}_1\Vert&= \Vert {{\mathbf {x}}}_4\Vert = r, \end{aligned}$$

(42)

$$\begin{aligned} \Vert {{\mathbf {x}}}_2\Vert&= \Vert {{\mathbf {x}}}_3\Vert = \frac{\sqrt{2}}{2} r. \end{aligned}$$

(43)

Fig. 27

A circular arc with $q= \frac{\pi }{2}$ represented by two quadratic NURBS elements in space

1.2 Prescribed velocity

We can represent the rotation-generated prescribed velocity ${\mathbf {g}}^h$ exactly along the circular arc. Because the circular arc is represented exactly, the velocity can be also represented exactly by

$$\begin{aligned} {\mathbf {g}}^h(\xi )&= \pmb {\omega } \times {\mathbf {x}}^h(\xi ) \end{aligned}$$

(44)

$$\begin{aligned}&= \pmb {\omega } \times \sum _{a=1}^3 N_a (\xi ) {\mathbf {x}}_a \end{aligned}$$

(45)

$$\begin{aligned}&= \sum _{a=1}^3 N_a (\xi ) \left( \pmb {\omega } \times {\mathbf {x}}_a\right) . \end{aligned}$$

(46)

This implies $\displaystyle {\mathbf {g}}^h(\xi ) = \sum _{a=1}^3 N_a(\xi ) {\mathbf {g}}_a$, where

$$\begin{aligned} {\mathbf {g}}_a&= \pmb {\omega } \times {\mathbf {x}}_a. \end{aligned}$$

(47)

This is applicable also for the case $q= \frac{\pi }{2}$.

1.3 An example: circular-arc path

This section is mostly from [4, 19, 22].

1.3.1 Path representation

NURBS temporal basis functions can be used in representing the circular-arc path of a moving point. Let us select as the origin the center of the circle to which the arc belongs. The point moves from ${{\mathbf {x}}}^1$ to ${{\mathbf {x}}}^3$, with $\Vert {{\mathbf {x}}}^1\Vert = \Vert {{\mathbf {x}}}^3\Vert$, as shown in Fig. 28. As it was done in space in Appendix C.1, the circular-arc path can be represented in time exactly by three control points with quadratic NURBS basis functions (it is valid only for $q< \frac{\pi }{2}$). The weights are $w^1 = w^3 = 1$, and $w^2 = \cos q$, where

$$\begin{aligned} \cos 2 q&= \frac{{{\mathbf {x}}}^1 \cdot {{\mathbf {x}}}^3}{r^2},&r&= \Vert {{\mathbf {x}}}^1\Vert = \Vert {{\mathbf {x}}}^3\Vert . \end{aligned}$$

(48)

Fig. 28

A circular-arc path represented by a quadratic NURBS element in time

This results in the temporal basis functions

$$\begin{aligned} T^1\left( \varTheta \right)&= \frac{(1-\varTheta )^2}{2\left( \left( 1+\varTheta ^2\right) +w^2\left( 1-\varTheta ^2\right) \right) }, \end{aligned}$$

(49)

$$\begin{aligned} T^2\left( \varTheta \right)&= \frac{w^2 (1-\varTheta ^2)}{\left( 1+\varTheta ^2\right) +w^2\left( 1-\varTheta ^2\right) }, \end{aligned}$$

(50)

$$\begin{aligned} T^3\left( \varTheta \right)&= \frac{(1+\varTheta )^2}{2\left( \left( 1+\varTheta ^2\right) +w^2\left( 1-\varTheta ^2\right) \right) }. \end{aligned}$$

(51)

and the control points are ${{\mathbf {x}}}^1$,

$$\begin{aligned} {{\mathbf {x}}}^2&= \frac{r}{w^2} \frac{{{\mathbf {x}}}^1 + {{\mathbf {x}}}^3}{\Vert {{\mathbf {x}}}^1 + {{\mathbf {x}}}^3 \Vert } \end{aligned}$$

(52)

$$\begin{aligned}&= \frac{1}{2 (w^2)^2} \left( {{\mathbf {x}}}^1 + {{\mathbf {x}}}^3\right) , \end{aligned}$$

(53)

and ${{\mathbf {x}}}^3$. Thus, the arc can be represented as

$$\begin{aligned} {{\mathbf {x}}}\left( \varTheta _{\mathbf {x}}\right) = {{\mathbf {x}}}^1 T^1\left( \varTheta _{\mathbf {x}}\right) +{{\mathbf {x}}}^2 T^2\left( \varTheta _{\mathbf {x}}\right) +{{\mathbf {x}}}^3 T^3\left( \varTheta _{\mathbf {x}}\right) . \end{aligned}$$

(54)

1.3.2 Constant angular velocity

First, by using Eq. (53), we rearrange Eq. (54) as

$$\begin{aligned} {{\mathbf {x}}}\left( \varTheta _{\mathbf {x}}\right)&= \underbrace{ \left( T^1(\varTheta _{\mathbf {x}}) + \frac{1}{2w_2^2}T^2(\varTheta _{\mathbf {x}})\right) }_{Q^1\left( \varTheta _{\mathbf {x}}\right) } {{\mathbf {x}}}^1 \nonumber \\&\quad + \underbrace{ \left( T^3(\varTheta _{\mathbf {x}}) + \frac{1}{2w_2^2}T^2(\varTheta _{\mathbf {x}})\right) }_{Q^3\left( \varTheta _{\mathbf {x}}\right) } {{\mathbf {x}}}^3 , \end{aligned}$$

(55)

where $Q^1$ and $Q^3$ are introduced for notation convenience. Taking the cross product with the unit vector along ${{\mathbf {x}}}^2$, we get

$$\begin{aligned} \frac{{{\mathbf {x}}}^1 + {{\mathbf {x}}}^3}{\Vert {{\mathbf {x}}}^1 + {{\mathbf {x}}}^3\Vert } \times {{\mathbf {x}}}({\varTheta _{\mathbf {x}}}) = \frac{{{\mathbf {x}}}^1\times {{\mathbf {x}}}^3}{r^2\sin (2 q)}{r \sin (\omega t)}, \end{aligned}$$

(56)

where

$$\begin{aligned} -\frac{\varDelta t}{2} \le t \le \frac{\varDelta t}{2}, \end{aligned}$$

(57)

and $\omega \varDelta t = 2 q$ for notation convenience. From Eq. (56), we get

$$\begin{aligned} \frac{{{\mathbf {x}}}^1 + {{\mathbf {x}}}^3}{2 r \cos q} \times {{\mathbf {x}}}({\varTheta _{\mathbf {x}}}) = \frac{{{\mathbf {x}}}^1\times {{\mathbf {x}}}^3}{2 r\sin q\cos q}{\sin (\omega t)}. \end{aligned}$$

(58)

From Eqs. (55) and (58), we get

$$\begin{aligned} Q^3 - Q^1 = \frac{\sin (\omega t)}{\sin q} . \end{aligned}$$

(59)

From Eq. (59), it can be shown that

$$\begin{aligned} \varTheta _{\mathbf {x}}= \frac{\sin q}{1-\cos q}\frac{\sin (\omega t)}{1+\cos (\omega t)}. \end{aligned}$$

(60)

We assume that time is represented with the same basis functions, and therefore

$$\begin{aligned} t(\varTheta _t)&= \frac{\varDelta t}{2}\left( T^3(\varTheta _t) - T^1(\varTheta _t)\right) \end{aligned}$$

(61)

$$\begin{aligned}&=\frac{\varDelta t \varTheta _t}{1 + \varTheta _t^2 + (1- \varTheta _t^2) \cos q}. \end{aligned}$$

(62)

We select $\frac{\mathrm {d}t}{\mathrm {d}\theta } = \frac{\varDelta t}{2}$, and with Eq. (57), we get $t=\frac{\varDelta t}{2} \theta$. We solve Eq. (61) with that on the left-hand side, and obtain $\varTheta _t$ as a function of $\theta$:

$$\begin{aligned} \varTheta _t(\theta )&= \frac{\theta (1+\cos q)}{1 + \sqrt{1 - (\theta \sin q)^2}}. \end{aligned}$$

(63)

Differentiating that with respect to $\theta$, we reach

$$\begin{aligned} \frac{\mathrm {d}\varTheta _t}{\mathrm {d}\theta }&= \frac{1 + \cos q}{\sqrt{1 - (\theta \sin q)^2} \left( 1 + \sqrt{1 - (\theta \sin q)^2}\right) }. \end{aligned}$$

(64)

The secondary mapping for the circular-arc path of the moving point becomes

$$\begin{aligned} \varTheta _{\mathbf {x}}= \frac{\sin q}{1-\cos q}\frac{\sin (q\theta )}{1+\cos (q\theta )}, \end{aligned}$$

(65)

and

$$\begin{aligned} \frac{\mathrm {d}\varTheta _{{\mathbf {x}}}}{\mathrm {d}\theta }&= \frac{\sin q}{1-\cos q} \frac{1 }{1+\cos (q\theta )} q . \end{aligned}$$

(66)

From that, the velocity of the point is calculated as

$$\begin{aligned} {\mathbf {v}}(\theta )&= \frac{\mathrm {d}{\mathbf {x}}}{\mathrm {d}\varTheta _{\mathbf {x}}} \frac{\mathrm {d}\varTheta _{\mathbf {x}}}{\mathrm {d}\theta } \frac{2}{\varDelta t}, \end{aligned}$$

(67)

which is an exact representation.

1.4 Prescribed velocity in the ST domain

We can also prescribe a rotation-generated velocity along the circular path exactly. Let us assume that we have a set of control points ${\mathbf {x}}_a^\alpha$ with corresponding basis functions $T^\alpha$ and $N_a$. A position in the ST domain is represented as

$$\begin{aligned} {\mathbf {x}}^h(\theta , \xi )&= \sum _{\alpha =1}^3 \sum _{a=1}^3 T^\alpha (\varTheta _{\mathbf {x}}(\theta )) N_a (\xi ) {\mathbf {x}}_a^\alpha . \end{aligned}$$

(68)

With that, we can express the prescribed velocity as

$$\begin{aligned} {\mathbf {g}}^h(\theta , \xi )&= \pmb {\omega } \times {\mathbf {x}}^h(\theta , \xi ) \end{aligned}$$

(69)

$$\begin{aligned}&= \pmb {\omega } \times \sum _{\alpha =1}^3 \sum _{a=1}^3 T^\alpha (\varTheta _{\mathbf {x}}(\theta )) N_a (\xi ) {\mathbf {x}}_a^\alpha \end{aligned}$$

(70)

$$\begin{aligned}&= \sum _{\alpha =1}^3 \sum _{a=1}^3 T^\alpha (\varTheta _{\mathbf {x}}(\theta )) N_a (\xi ) \left( \pmb {\omega } \times {\mathbf {x}}_a^\alpha \right) . \end{aligned}$$

(71)

This implies $\displaystyle {\mathbf {g}}^h(\theta , \xi ) = \sum _{\alpha =1}^3 \sum _{a=1}^3 T^\alpha (\varTheta _{\mathbf {x}}(\theta ) ) N_a(\xi ) {\mathbf {g}}_a^\alpha$, where

$$\begin{aligned} {\mathbf {g}}_a^\alpha&= \pmb {\omega } \times {\mathbf {x}}_a^\alpha . \end{aligned}$$

(72)

Remark 7

As shown above, the exact representation of the prescribed velocity when the mesh is rotating is possible with the secondary mappings $\varTheta _t$ and $\varTheta _{\mathbf {x}}$ based on quadratic NURBS in time. However, if the flow solution is based on linear functions in time, then the representation of the prescribed velocity will not be exact in the flow solution. The accuracy of the representation will decrease with increasing rotation per time step.