U-duct turbulent-flow computation with the ST-VMS method and isogeometric discretization

Abstract

The U-duct turbulent flow is a known benchmark problem with the computational challenges of high Reynolds number, high curvature and strong flow dependence on the inflow profile. We use this benchmark problem to test and evaluate the Space–Time Variational Multiscale (ST-VMS) method with ST isogeometric discretization. A fully-developed flow field in a straight duct with periodicity condition is used as the inflow profile. The ST-VMS serves as the core method. 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 unsteady flow. The ST isogeometric discretization enables more accurate representation of the duct geometry and increased accuracy in the flow solution. In the straight-duct computations to obtain the inflow velocity, the periodicity condition is enforced with the ST Slip Interface method. All computations are carried out with quadratic NURBS meshes, which represent the circular arc of the duct exactly in the U-duct computations. We investigate how the results vary with the time-averaging range used in reporting the results, mesh refinement, and the Courant number. The results are compared to experimental data, showing that the ST-VMS with ST isogeometric discretization provides good accuracy in this class of flow problems.

1 Introduction

In the benchmarking context of the U-duct turbulent flow, which has a number of computational challenges, we conduct test and evaluation of the Space–Time Variational Multiscale (ST-VMS) method [1,2,3] with ST isogeometric discretization [1, 4,5,6]. Turbulent-flow test and evaluation studies were conducted earlier for the ST-VMS (see [5, 7, 8]), but the computations in [7, 8] were with finite element discretization, and the computation in [5] was with a significantly milder test problem, a straight pipe with circular cross-section. Furthermore, here we use the latest stabilization parameters [9] designed in conjunction with the ST-VMS, with the latest element length expressions [10].

1.1 Stabilized and VMS ST computational methods

The stabilized and VMS ST computational methods started with the Deforming-Spatial-Domain/Stabilized ST (DSD/SST) method [11,12,13]. The DSD/SST 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) [14, 15] and compressible-flow SUPG [16,17,18] methods are two of the earliest and most widely used stabilized methods. The incompressible-flow Pressure-Stabilizing/Petrov-Galerkin (PSPG) method [11, 19], with its Stokes-flow version introduced in [20], 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) [21,22,23,24], 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 [1,2,3] 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, 2]), 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, [25]) and more commonly used. The ALE-VMS method [26,27,28,29,30,31,32] is the VMS version of the ALE. It succeeded the ST-SUPS and ALE-SUPS [33] 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 [33], wind turbine aerodynamics and FSI [34,35,36,37,38,39,40,41,42,43,44,45,46], more specifically, vertical-axis wind turbines (VAWTs) [43, 44, 47, 48], floating wind turbines [49], wind turbines in atmospheric boundary layers [42,43,44, 50,51,52], and fatigue damage in wind turbine blades [53], patient-specific cardiovascular fluid mechanics and FSI [26, 54,55,56,57,58,59], biomedical-device FSI [60,61,62,63,64,65,66,67], ship hydrodynamics with free-surface flow and fluid–object interaction [68, 69], hydrodynamics and FSI of a hydraulic arresting gear [70, 71], hydrodynamics of tidal-stream turbines with free-surface flow [72], passive-morphing FSI in turbomachinery [73], bioinspired FSI for marine propulsion [74, 75], bridge aerodynamics and fluid–object interaction [76,77,78], and mixed ALE-VMS/Immersogeometric computations [63,64,65, 79, 80] in the framework of the Fluid–Solid Interface-Tracking/Interface-Capturing Technique [81]. Recent advances in stabilized and multiscale methods may be found for stratified incompressible flows in [82], for divergence-conforming discretizations of incompressible flows in [83], and for compressible flows with emphasis on gas-turbine modeling in [84, 85].

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 [86,87,88,89], was the Jacobian-based stiffening (JBS), which is now called, for reasons explained in [90, 91], “mesh-Jacobian-based stiffening.” The most recent ones are the element-based mesh relaxation [92], 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 [91] based on fiber-reinforced hyperelasticity and optimized ZSS, and a linear-elasticity-based mesh moving method with no cycle-to-cycle accumulated distortion [93].

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

1.2 ST Slip Interface method

The ST Slip Interface (ST-SI) method was introduced in [106], 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 [141, 142]. 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 [106]. 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 [127] 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 [43,44,45,46, 106, 107], thermo-fluid analysis of disk brakes [127], flow-driven string dynamics in turbomachinery [45, 46, 128,129,130], flow analysis of turbocharger turbines [5, 131,132,133,134], flow around tires with road contact and deformation [9, 122, 135,136,137], fluid films [137, 138], aerodynamic analysis of ram-air parachutes [6, 51, 52], and flow analysis of heart valves [66, 67, 119, 123,124,125] and ventricle-valve-aorta sequences [90]. In the ST-SI version presented in [106] 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 [6, 51, 52].

1.3 ST Isogeometric Analysis

The success with Isogeometric Analysis (IGA) basis functions in space [26, 54, 141, 143] 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, 2] and demonstrated in [4, 108, 110]. It also enables more efficient temporal representation of the motion and deformation of the volume meshes, and more efficient remeshing. These motivated the development of the ST/NURBS Mesh Update Method (STNMUM) [4, 108, 110], with the name coined in [103]. 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, 2, 4, 29] 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 [3, 122, 127,128,129,130, 144]. 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, 29, 108, 109], bioinspired MAVs [104, 105, 110, 111] and wing-clapping [112, 113], separation aerodynamics of spacecraft [98], aerodynamics of horizontal-axis [40, 103,104,105] and vertical-axis [43,44,45,46, 106, 107] wind turbines, thermo-fluid analysis of ground vehicles and their tires [3, 51, 122], thermo-fluid analysis of disk brakes [127], flow-driven string dynamics in turbomachinery [45, 46, 128,129,130], flow analysis of turbocharger turbines [5, 131,132,133,134], and flow analysis of coronary arteries in motion [126].

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 [5, 131,132,133,134], flow-driven string dynamics in turbomachinery [45, 46, 129, 130], ram-air parachutes [6, 51, 52], spacecraft parachutes [140], aortas [66, 67, 119, 120], heart valves [66, 67, 119, 123,124,125], ventricle-valve-aorta sequences [90], coronary arteries in motion [126], tires with road contact and deformation [9, 136, 137], fluid films [137, 138], and VAWTs [45, 46, 107]. The image-based arterial geometries used in patient-specific arterial FSI computations do not come from the ZSS of the artery. A number of methods [27, 29, 145,146,147,148,149,150,151,152,153,154] have been proposed for estimating the ZSS required in the computations. Using IGA basis functions in space is now a key part of some of the newest ZSS estimation methods [66, 152,153,154,155] and related shell analysis [156]. The IGA has also been successfully applied to the structural analysis of wind turbine blades [157,158,159,160,161].

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 [29]). 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 [14, 15] and [16,17,18], followed by the ones introduced in [162, 163]. In many cases, the element length was seen as an advection length scale, in the flow-velocity direction. The $\tau$ definition introduced in [163], 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 (see [164, 165]).

Calculating the $\tau$s based on the element-level matrices and vectors was introduced in [166] 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 [166, 167] 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 [168, 169]. Recognizing this as a diffusion length scale, a new stabilization parameter for the diffusive limit, “$\tau _{\mathrm {SUGN3}}$,” was introduced in [169, 170], 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 [168, 169].

Some new options for the stabilization parameters used with the SUPS and VMS were proposed in [3, 4, 7, 103, 171]. These include a fourth $\tau$ component, “$\tau _{\mathrm {SUGN4}}$” [3], 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 [3] 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 [73, 172,173,174,175,176,177,178,179,180,181,182,183,184].

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 [185] 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 [185] 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 [9], which were mostly from [185], 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 [186] 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 [186] 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 [10]. 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 [187] that the local-length-scale expressions introduced in [10] meet that requirement.

The direction-dependent local-length-scale expressions introduced in [10, 185] have been used in computational flow analysis of turbocharger turbines [133, 134], compressible-flow spacecraft parachutes [140], tires with road contact and deformation [9, 137], fluid films [137, 138], heart valves [125], ventricle-valve-aorta sequences [90], and tsunami-shelter VAWTs [107]. They have also been used in [85], in the context of gas turbine computational flow analysis with isogeometric discretization, in calculating the Courant number based on the NURBS mesh local length scale in the flow direction.

1.5 U-duct turbulent flow

The U-duct turbulent flow is a known benchmark problem with sufficient experimental data (see, for example, [188]). The computational challenges encountered include high Reynolds numbers, high curvature and strong flow dependence on the inflow profile. In the U-duct computations we present, a fully-developed flow field in a straight duct with periodicity condition is used as the inflow profile. In the straight-duct computations to obtain the inflow velocity, the periodicity condition is enforced with the ST-SI method. Both the straight-duct and U-duct computations are carried out with quadratic NURBS meshes, which represent the circular arc of the duct exactly in the U-duct computations. We investigate how the results vary with the time-averaging range used in reporting the results, mesh refinement, and the Courant number. The results are compared to experimental data [188] to show how the ST-VMS with ST isogeometric discretization performs in this class of flow problems.

1.6 Outline of the remaining sections

In Sect. 2, we provide the definitions used in the data analysis. The straight-duct computations are presented in Sect. 3, and the U-duct computations in Sect. 4. The concluding remarks are given in Sect. 5. The ST-VMS and ST-SI methods and the stabilization parameters are given in Appendices A and B.

2 Definitions for the data analysis

2.1 Scale separation

We split the velocity scales as

$$\begin{aligned} {\mathbf {u}}&= {\overline{{\mathbf {u}}}} + {\mathbf {u}}^\prime , \end{aligned}$$

(1)

where the overbar indicates the time-averaging over the range ${\mathcal {T}} = (T_1, T_2)$:

$$\begin{aligned} {\overline{f}} = \dfrac{1}{T_2 - T_1} \int _{T_1}^{T_2} f \mathrm{d}t, \end{aligned}$$

(2)

and f can be a vector or scalar. We note that this scale separation is different from the VMS scale separation. It is used only for post processing. We extend that to the ST context as

$$\begin{aligned} \overline{\overline{f}}&= \frac{1}{\int _{{\mathcal {T}}} \int _{\varOmega _t} \mathrm{d}\varOmega \mathrm{d}t} \int _{{\mathcal {T}}} \int _{\varOmega _t} f \mathrm{d}\varOmega \mathrm{d}t \end{aligned}$$

(3)

$$\begin{aligned}&= \frac{1}{\int _Q \mathrm{d}Q} \int _Q f \mathrm{d}Q , \end{aligned}$$

(4)

where $Q = \left\{ \left. {\mathbf {x}}(t) \in \varOmega _t \right| t \in {\mathcal {T}} \right\}$ is the ST domain.

We define the $L_q$ norm of a scalar as

$$\begin{aligned} \left\| f \right\| _{q,{\mathcal {T}}} = \left( \dfrac{1}{T_2 - T_1} \int _{T_1}^{T_2} \left| f \right| ^q \mathrm{d}t \right) ^{\frac{1}{q}}. \end{aligned}$$

(5)

We extend that to the ST context as

$$\begin{aligned} \left\| f \right\| _{q,Q}&= \left( \frac{1}{\int _{{\mathcal {T}}} \int _{\varOmega _t} \mathrm{d}\varOmega \mathrm{d}t} \int _{{\mathcal {T}}} \int _{\varOmega _t} \left| f \right| ^q \mathrm{d}\varOmega \mathrm{d}t \right) ^{\frac{1}{q}}\end{aligned}$$

(6)

$$\begin{aligned}&= \left( \frac{1}{\int _Q \mathrm{d}Q} \int _Q \left| f \right| ^q \mathrm{d}Q \right) ^{\frac{1}{q}}. \end{aligned}$$

(7)

2.2 Nondimensionalization

With $\rho$, U and D being the scales for the density, velocity and length, we define the scaled quantities as

$$\begin{aligned} {\mathbf {u}}^*&= \frac{{\mathbf {u}}}{U}, \end{aligned}$$

(8)

$$\begin{aligned} t^*&= \frac{t U}{D}, \end{aligned}$$

(9)

$$\begin{aligned} p^*&= \frac{p }{\rho U^2}, \end{aligned}$$

(10)

where p is the pressure. For notation convenience, we drop the asterisk, which results in $\rho =1$, for example.

2.3 Wall-related scaling

We define the scaled wall-normal coordinate as

$$\begin{aligned} y^+&= \frac{y u_\tau }{\nu }, \end{aligned}$$

(11)

where y is the coordinate along the wall normal, $\nu = \mu /\rho$ is the kinematic viscosity, and $\mu$ is the viscosity. The friction velocity $u_\tau$ is defined as

$$\begin{aligned} u_\tau = \sqrt{\frac{\left\| {\mathbf {h}}_\mathrm {v}\right\| }{ \rho }}, \end{aligned}$$

(12)

where ${\mathbf {h}}_\mathrm {v}$ is the wall shear stress. The streamwise velocity is scaled near the wall as $u^+ = \frac{u_\mathrm {s}}{u_\tau }$. We note that whether the scaling we do here is done after or before the scaling in Sect. 2.2 does not change the outcome.

3 Straight duct with periodicity condition

3.1 Problem setup

The duct has a $D{\times }D$ cross-section and is 5D long. Figure 1 shows the computational domain and the coordinate basis vectors ${\mathbf {e}}_1$, ${\mathbf {e}}_2$ and ${\mathbf {e}}_\mathrm {s}$. In the Reynolds number definition $\mathrm {Re}=\frac{U D}{\nu }$, U is streamwise velocity averaged in time and over the cross-section. We compute the cases $\mathrm {Re}=4{\times }10^4$ and $10^5$. In the data analysis, we express the velocity components as $u_k = {\mathbf {u}} \cdot {\mathbf {e}}_k$, where $k=1, 2, \mathrm {s}$.

Fig. 1

Straight duct. Computational domain and the coordinate basis vectors. The two end planes are colored red to remind the reader that we are enforcing periodicity condition on those planes. (Color figure online)

3.2 Mesh

Figure 2 shows the mesh used for both $\mathrm {Re} = 4{\times }10^4$ and $10^5$, which is made of $72^3$ quadratic NURBS elements. The mesh is uniform in the streamwise direction. The normal-direction thickness for the first layer of elements near the wall results in $y^+= 0.43$ and 0.95 for $\mathrm {Re}=4{\times }10^4$ and $10^5$. In calculating the $y^+$ values based on Eqs. (11) and (12), $\left\| {\mathbf {h}}_\mathrm {v}\right\|$ is estimated from the pipe friction factor f given [189] as

$$\begin{aligned} \dfrac{1}{f^{0.5}} = 2\log \left( \mathrm {Re}~f^{0.5}\right) -0.8. \end{aligned}$$

(13)

The f values corresponding to $\mathrm {Re} = 4{\times }10^4$ and $10^5$ are $2.2{\times }10^{-2}$ and $1.8{\times }10^{-2}$.

Fig. 2

Straight duct. Control mesh (left). The yellow circles are the control points. The corresponding mesh (middle) and its part in the red frame (right). (Color figure online)

3.3 Boundary conditions

The no-slip conditions on the walls are enforced weakly (see Appendix A.2.2). In each case of Reynolds number, the pressure difference specified across the SI (see Appendix A.2.1) is adjusted until we achieve a Reynolds number close enough to the case Reynolds number. That becomes the actual Reynolds number we compute.

3.4 Computational conditions

We use the ST-VMS, with the stabilization parameters given by Eqs. (22)–(24), (33) and (34). The time-step size $\varDelta t$ is determined from the Courant number $C_{\varDelta t} = \frac{U \varDelta t}{h_\mathrm {s}}$, where $h_\mathrm {s}$ is the “apparent” element length¹ in the streamwise direction. We set $C_{\varDelta t} = 0.322$. The number of nonlinear iterations per time step is 3, and the number of GMRES iterations per nonlinear iteration is 500. We define $T = L/U$, where L is the duct length. After achieving the actual Reynolds number we compute, the duration we compute for data extraction is 20T, which is equivalent to 4474 time steps.

3.5 Results

In computations we carry out to adjust the pressure difference specified across the SI, we achieve the Reynolds number values of $3.96{\times }10^4$ and $9.98{\times }10^4$ at scaled pressure gradient values $5.6{\times }10^{-2}$ and $4.7{\times }10^{-2}$. In both cases, the difference from the case Reynolds number is no more than $1\%$.

Figure 3 shows, for $\mathrm {Re} = 10^5$, $\overline{\overline{u}}_\mathrm {s}$ and $\left\| u'_\mathrm {s}\right\| _{2,Q}$, together with the experimental data from [188]. Figure 4 shows, for both $\mathrm {Re}=4{\times }10^4$ and $10^5$, $\overline{\overline{u^+}}$, together with the DNS data for $\mathrm {Re}=4{\times }10^4$ from [190].

Fig. 3

Straight duct. $\mathrm {Re}=10^5$. $\overline{\overline{u}}_\mathrm {s}$ and $\left\| u_s'\right\| _{2, Q}$ with $Q = \left\{ {\mathbf {x}}=\left. (0, x_2, x_\mathrm {s}) \right| x_\mathrm {s} \in (0, 5D), t \in (8.5 T , 20 T )\right\}$. The experimental data is from [188]

Fig. 4

Straight duct. $\mathrm {Re}=4{\times }10^4$ and $\mathrm {Re}=10^5$. $\overline{\overline{u^+}}$ with $Q = \left\{ \left. {\mathbf {x}}=(0, x_2, x_\mathrm {s}) \right| x_s \in (0, 5D), t \in (8.5 T , 20 T )\right\}$. The DNS data is for $\mathrm {Re} = 4{\times }10^4$ and from [190]

4 U-duct

4.1 Problem setup

Figure 5 shows the computational domain and the coordinate basis vectors. We compute the case $\mathrm {Re}=10^5$. The locations where the flow characteristics are reported are shown in Fig. 6.

Fig. 5

U-duct. Computational domain and the coordinate basis vectors. The red plane is the inlet, and the blue plane is the outlet. (Color figure online)

Fig. 6

U-duct. Locations where the flow characteristics are reported. The red and blue planes are the near-wall ($x_1/D=0.375$) and center ($x_1=0$) planes. The velocity profiles are reported along the red and blue lines. (Color figure online)

Fig. 7

U-duct. Exact representation of the curvature by NURBS. Weight definitions for the mesh. The blue points represent the control point locations where the weight is ${\text {cos}(\pi /4)}$, and the red points are where it is 1. (Color figure online)

4.2 Boundary conditions

The ST-averaged velocity profile, $\overline{\overline{{\mathbf {u}}}}$, from the straight-duct computation is used as the inlet condition (see Sect. 3.5) The no-slip conditions on the walls are enforced weakly. At the outlet, zero-stress condition is used.

Fig. 8

U-duct. Mesh A, B, C, D and E

Table 1 U-duct. Number of control points (nc), number of quadratic NURBS elements in the entire model (ne) and along the bend ($ne_{\mathrm {b}}$), streamwise length of the first element at the inlet ($h_{\mathrm {s}, \mathrm {inlet}}$), minimum streamwise element length ($h_{\mathrm {s}, \mathrm {min}}$), which is in the bend, and $y^+$ at the walls

Fig. 9

U-duct. Flow development as depicted by $\left\| \overline{{\mathbf {u}}}\right\|$ and ${\overline{p}}$ on the center plane, for time-averaging over ${\mathcal {T}} = (0, 10 T)$, (10T, 20T), (20T, 30T), and (30T, 40T)

4.3 Mesh

To define the geometry exactly, we use four quadratic NURBS patches. Figure 7 shows the weight definitions. By using a sequence of knot insertions, we generate five meshes, which we call Mesh A, B, C, D and E. We start with Mesh A. Mesh B, C and D have in all three directions 2, 3 and 4 times the number of elements Mesh A has. Mesh E has along the bend twice the number of elements Mesh D has. Beyond that it has 5 more elements along the lower straight part of the duct, and the element sizes in both the lower and upper straight parts are adjusted along the streamwise direction such that the maximum ratio between two adjacent elements is at most 2. We note that Mesh D and E have the same cross-section as the straight-duct mesh in Sect. 3, and the same streamwise element length at the inlet. Figure 8 shows all five meshes, and Table 1 shows the mesh data for all five.

4.4 Computational conditions

We use the ST-VMS, with the stabilization parameters given by Eqs. (22)–(24), (33) and (34). The time-step size is determined from the Courant number, which is based on the minimum streamwise element length (see Table 1). The number of nonlinear iterations per time step is 3, and the number of GMRES iterations per nonlinear iteration is 500. In calculating T from $T=L/U$, we define $L=(3+0.65\pi +6.5)D$.

4.5 Results

4.5.1 Sequence of computations

First, Mesh A is used in computing for ${\mathcal {T}} = (0, 40T)$. Then, the initial condition for Mesh B computation is obtained by refining the data from Mesh A at $t = 29T$. After computing with Mesh B for a duration of T, the initial condition for Mesh C is obtained by least-squares projection at $t = 30T$. The initial conditions for Mesh D and E are also obtained by least-squares projection, from Mesh C at $t=31T$ and from Mesh D at $t=32T$.

4.5.2 Effect of the time-averaging range

In this study, we use $C_{\varDelta t} = 10$. Figure 9 shows the flow development as depicted by the flow patterns for different time-averaging ranges, all spanning 10T. The differences between the flow patterns for the time ranges beyond ${\mathcal {T}} = (0, 10T)$ are not significant. Figure 10 shows the Fourier decomposition of $u_2$ in various time ranges, all spanning 10T. The patterns in the time ranges beyond ${\mathcal {T}} = (0, 10T)$ are similar. Figure 11 shows the Fourier decomposition of $u_2$ in various time ranges ending at 40T. We see that the lowest frequency of local maximum is around $0.67~T^{-1}$, which corresponds to a period of 1.5T. From that we conclude that an averaging period of 3T is long enough in taking the statistics of the flow field.

Fig. 10

U-duct. Fourier decomposition of $u_\mathrm {2}$ at the blue point on the center plane, in the time ranges indicated in the legend, all spanning 10T . (Color figure online)

Fig. 11

U-duct. Fourier decomposition of $u_\mathrm {2}$ at the blue point on the center plane, in the time ranges indicated in the legend, all ending at 40T. The cyan line marks the lowest frequency of local maximum. (Color figure online)

4.5.3 Effect of the mesh refinement

We compare the data computed with Mesh A to E at $C_{\varDelta t} = 10$. Figure 12 shows the Fourier decomposition of $u_2$ in ${\mathcal {T}} = (33 T, 36 T)$. We see good agreement, though with the representation getting shorter at the higher end of the spectrum as the mesh gets coarser. Figure 13 shows ${\overline{u}}_\mathrm {s}$ on the center plane, time-averaged over ${\mathcal {T}}=(33T, 36T)$, together with the experimental data from [188]. Slight effect of mesh refinement is seen between the results obtained with Mesh A and the finer meshes up to Mesh D, while the results obtained with Mesh B, C and D are in agreement. The results from Mesh E differ from the others slightly around the bend. Section 4.5.4 includes more investigation of the differences between the results obtained with Mesh D and E.

Fig. 12

U-duct. Fourier decomposition of $u_\mathrm {2}$ at the blue point on the center plane, in ${\mathcal {T}} = (33 T, 36 T)$. (Color figure online)

Fig. 13

U-duct. Effect of the mesh refinement. ${\overline{u}}_\mathrm {s}$ on the center plane, at various locations along the duct, time-averaged over ${\mathcal {T}} = (33 T, 36 T)$. The experimental data is from [188]

Fig. 14

U-duct. Isosurfaces corresponding to a positive value of the second invariant of $\pmb {\nabla }{\overline{{\mathbf {u}}}}$, colored by $\left\| \overline{{\mathbf {u}}}\right\|$, time-averaged over ${\mathcal {T}} = (33 T, 36 T)$. Mesh D (left) and E (right). $C_{\varDelta t} = 10$ (top), 5 (middle) and 2.5 (bottom). The yellow contour lines represent the intersection between the isosurfaces and the center plane. (Color figure online)

4.5.4 Effect of the Courant number

We compare the data computed with Mesh D and E at $C_{\varDelta t} = 10$, 5 and 2.5, in ${\mathcal {T}} = (33 T, 36 T)$. Figure 14 shows the isosurfaces corresponding to a positive value of the second invariant of $\pmb {\nabla }{\overline{{\mathbf {u}}}}$. With Mesh D, the effect of the Courant number is not significant. With Mesh E at $C_{\varDelta t} = 10$, the behavior is similar to what we see with Mesh D. On the other hand, with Mesh E at $C_{\varDelta t} = 5$ and 2.5, the recirculation occurs earlier. This observation becomes easier to make when we inspect ${\overline{u}}_\mathrm {s}$. Figures 15 and 16 show that, together with the experimental data from [188], on the center and near-wall planes. Figures 17 and 18 show $\left\| u_\mathrm {s}' \right\| _{2,{\mathcal {T}}}$ on those planes. Overall, the data computed with Mesh E at $C_{\varDelta t} = 5$ and 2.5 is in good agreement with the experimental data.

5 Concluding remarks

We have conducted test and evaluation of the ST-VMS with ST isogeometric discretization in the benchmarking context of the U-duct turbulent flow, which has a number of computational challenges, and there is a good amount of experimental and computational data associated with this benchmark problem. The computational challenges include high Reynolds numbers, high curvature and strong flow dependence on the inflow profile. 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 unsteady flow. The ST isogeometric discretization enables more accurate representation of the duct geometry and increased accuracy in the flow solution. We used the latest stabilization parameters designed in conjunction with the ST-VMS, with the latest element length expressions. The inflow velocity used in the computations comes from a fully-developed flow field in a straight-duct computation with periodicity condition, enforced with the ST-SI method. Both the straight-duct and U-duct computations were carried out with quadratic NURBS meshes, enabling exact representation of the circular arc in the U-duct computations. We investigated how the results vary with the averaging period used in reporting the results, mesh refinement, and the Courant number. We compared the results to experimental data and showed that the ST-VMS with ST isogeometric discretization provides good accuracy in this class of flow problems.

Fig. 15

U-duct. Effect of the Courant number for Mesh D and E. ${\overline{u}}_\mathrm {s}$ on the center plane, at various locations along the duct, time-averaged over ${\mathcal {T}} = (33 T, 36 T)$. The experimental data is from [188]

Fig. 16

U-duct. Effect of the Courant number for Mesh D and E. ${\overline{u}}_\mathrm {s}$ on the near-wall plane, at various locations along the duct, time-averaged over ${\mathcal {T}} = (33 T, 36 T)$. The experimental data is from [188]

Fig. 17

U-duct. Effect of the Courant number for Mesh D and E. $\left\| {u}_\mathrm {s}'\right\| _{2,{\mathcal {T}}}$ on the center plane, at various locations along the duct, time-averaged over ${\mathcal {T}} = (33 T, 36 T)$. The experimental data is from [188]

Fig. 18

U-duct. Effect of the Courant number for Mesh D and E. $\left\| {u}_\mathrm {s}'\right\| _{2,{\mathcal {T}}}$ on the near-wall plane, at various locations along the duct, time-averaged over ${\mathcal {T}} = (33 T, 36 T)$. The experimental data is from [188]

Appendices

ST-VMS and ST-SI

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

1.1 ST-VMS

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 \\&\qquad + \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 \\&\qquad +\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 \\&\qquad + \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) \nonumber \right. \\&\left. \qquad +\pmb {\nabla }q^h \right) \cdot {\mathbf {r}}_{\mathrm {M}}^h \mathrm{d}Q\nonumber \\&\qquad + \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 \\&\qquad -\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 \\&\qquad -\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 = 0, \end{aligned}$$

(14)

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}$$

(15)

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

(16)

are the residuals of the momentum equation and incompressibility constraint. Here, ${\mathbf {f}}$ and ${\mathbf {h}}$ are the body force and the stress specified at the boundary. 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) ^T\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 1

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.

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 Eq. (14) includes added boundary terms corresponding to the SI. 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} - {q}_\mathrm {A}^h {{\mathbf {n}}}_\mathrm {A} \right) \cdot \frac{1}{2} \left( {{\mathbf {u}}}_\mathrm {B}^h - {{\mathbf {u}}}_\mathrm {A}^h \right) \mathrm{d}P \nonumber \\&- \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. - \left( {{\mathcal {F}}}_\mathrm {B}^h - \left| {{\mathcal {F}}}_\mathrm {B}^h \right| \right) {{\mathbf {u}}}_\mathrm {A}^h \right) \mathrm{d}P \nonumber \\&- \int _{\left( P_n\right) _\mathrm {SI}} \rho {{\mathbf {w}}}_\mathrm {A}^h \cdot \frac{1}{2} \left( \left( {{\mathcal {F}}}_\mathrm {A}^h - \left| {{\mathcal {F}}}_\mathrm {A}^h \right| \right) {{\mathbf {u}}}_\mathrm {A}^h \right. \nonumber \\&\quad \left. - \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 {w}}}_\mathrm {A}^h \right) \cdot \left( {{\hat{{\mathbf {n}}}}}_\mathrm {B} \cdot \mu \left( \pmb {\varepsilon }({{\mathbf {u}}}_\mathrm {B}^h) + \pmb {\varepsilon }({{\mathbf {u}}}_\mathrm {A}^h) \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) + \pmb {\varepsilon }\left( {{\mathbf {w}}}_\mathrm {A}^h \right) \right) \cdot \left( {{\mathbf {u}}}_\mathrm {B}^h - {{\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 {w}}}_\mathrm {A}^h \right) \cdot \left( {{\mathbf {u}}}_\mathrm {B}^h - {{\mathbf {u}}}_\mathrm {A}^h \right) \mathrm{d}P, \end{aligned}$$

(17)

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}$$

(18)

$$\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}$$

(19)

$$\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}$$

(20)

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 [135], 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. (17), 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, 106, 139, 140], positive $\varDelta p$ means higher pressure on Side A.

1.2.2 One-side formulation (fluid–solid SI)

On solid surfaces where we prefer weak enforcement of the Dirichlet conditions [191, 192] for the fluid, we use the ST-SI version where the SI is between the fluid and solid domains. This version is obtained (see [106]) 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. (14) 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 \nonumber \\&\quad - \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\nonumber \right. \\&\left. \quad + \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}$$

(21)

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

Stabilization parameters

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 [185]:

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

(22)

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}$$

(23)

and

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

(24)

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}$$

(25)

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}$$

(26)

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

(27)

where

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

(28)

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

(29)

The ST and space-only Jacobian tensors are

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

(30)

and

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

(31)

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 &{} {\mathbf {0}}^T \\ {\mathbf {0}}&{} {\mathbf {D}} \\ \end{bmatrix}. \end{aligned}$$

(32)

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

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

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

(33)

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

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

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

(34)

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}$$

(35)

For more ways of calculating the stabilization parameters in flow computations, see [73, 167, 172,173,174,175,176,177,178,179,180,181,182,183].

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}$$

(36)

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

(37)

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

(38)

$$\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}$$

(39)