A node-numbering-invariant directional length scale for simplex elements

Variational multiscale methods, and their precursors, stabilized methods, have been very popular in flow computations in the past several decades. Stabilization parameters embedded in most of these methods play a significant role. The parameters almost always involve element length scales, most of the time in specific directions, such as the direction of the flow or solution gradient. We require the length scales, including the directional length scales, to have node-numbering invariance for all element types, including simplex elements. We propose a length scale expression meeting that requirement. We analytically evaluate the expression in the context of simplex elements and compared to one of the most widely used length scale expressions and show the levels of noninvariance avoided.


Introduction
In this paper, we introduce a node-numbering-invariant directional length scale for simplex elements used in flow computations with the stabilized and variational multiscale (VMS) methods, discontinuity-capturing (DC) methods, and other special methods that require a directional element length.

Stabilized and VMS methods
Stabilized and VMS methods have been playing a crucial role in computational flow analysis. Two of the earliest and most widely used stabilized methods are the incompressible-flow Streamline-Upwind/Petrov-Galerkin (SUPG) 1,2 and compressible-flow SUPG 3-5 methods. The incompressible-flow Pressure-Stabilizing/Petrov-Galerkin (PSPG) method, 6,7 with its Stokes-flow version introduced in Ref. 8, also has a long history and is 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 they are residual-based, they accomplish stabilization without loss of accuracy. The residual-based VMS (RBVMS), [9][10][11][12] which is now also widely used, subsumes its precursor SUPG/PSPG.

DC methods
When the flow field has a shock or some other discontinuity, stabilized methods are often supplemented with a DC term. Supplementing the SUPG with a DC term goes back to 1986. 13,14 The DC term played a key role in the evolution of the compressible-flow SUPG, which was originally introduced in 1982 in the context of conservation variables. That 1982 method is now called "(SUPG) 82 ". At first, (SUPG) 82 was not used with any DC (shock-capturing) term, and the test computations clearly showed the need for additional stabilization at the shocks. Later, (SUPG) 82 was recast in entropy variables, but also supplemented with a DC term. 15 This of course resulted in better shock profiles. In an ASME paper in 1991, 16 (SUPG) 82 was supplemented with a very similar DC term. It was shown in Refs. 16 and 17 that, with the added DC term, (SUPG) 82 was very comparable in accuracy to (SUPG) 82 recast in entropy variables. The stabilized and DC methods introduced in Ref. 14 for the advection-diffusion-reaction equation accounted for the interaction between the DC and SUPG terms. Taking that interaction into account precludes "compounding" (i.e. augmentation of the SUPG effect by the DC effect when the advection and discontinuity directions are not orthogonal).

ST Slip Interface method
The ST Slip Interface (ST-SI) method was introduced in Ref. 98, in the context of incompressible-flow equations, to retain the desirable moving-mesh features of the ST-VMS and ST-SUPS when we have 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. The starting point in the development of the ST-SI was the version of the ALE-VMS for computations with sliding interfaces. 36,37 Interface terms similar to those in the ALE-VMS version are added to the ST-VMS to account for the compatibility conditions for the velocity and stress at the SI. That accurately connects the two sides of the solution. An ST-SI version where the SI is between fluid and solid domains was also presented in Ref. 98. The SI in this case is a "fluid-solid SI" rather than a standard "fluidfluid SI" and enables weak enforcement of the Dirichlet boundary conditions for the fluid. The ST-SI introduced in Ref. 119 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 vertical-axis wind turbines, 98,48,49 thermo-fluid analysis of disk brakes, 119 flow-driven string dynamics in turbomachinery, 120-122 flow analysis of turbocharger turbines, 123-126 flow around tires with road contact and deformation, 115,127-130 fluid films, 131,130 aerodynamic analysis of ram-air parachutes, 132 and flow analysis of heart valves. 116,117,111,113 In the ST-SI version presented in Ref. 98, 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. 132 The compressible-flow ST-SI methods were introduced in Ref. 133, including the version where the SI is between a thin porous structure and the fluid on its two sides. Compressible-flow porosity models were also introduced in Ref. 133. These, together with the compressible-flow ST SUPG method, 135 extended the ST computational analysis range to compressible-flow aerodynamics of parachutes with fabric A node-numbering-invariant directional length scale for simplex elements 2723 and geometric porosities. That enabled ST computational flow analysis of the Orion spacecraft drogue parachute in the compressible-flow regime. 133,134

ST Isogeometric Analysis
The success of using Isogeometric Analysis (IGA) basis functions in space 136,55,25,36 motivated the integration of the ST methods with isogeometric discretization, which we broadly call "ST-IGA". The ST-IGA was introduced in Ref. 22. Computations with the ST-VMS and ST-IGA were first reported in Ref. 22 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 getting full benefit out of using higher-order basis functions in space. This was demonstrated with the stability and accuracy analysis given in Ref. 22 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 Refs. 22 and 23 and demonstrated in Refs. 99, 100 and 102. 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) 99,100,102 ; the name was given in Ref. 95. 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" 22,99,23,28 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. 137,24,115,[119][120][121][122] 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, 99,100,28,101 bioinspired MAVs 102,103,96,97 and wing-clapping, 104,105 separation aerodynamics of spacecraft, 89 aerodynamics of horizontal-axis [95][96][97]45 and vertical-axis 98,48,49 windturbines, thermo-fluid analysis of ground vehicles and their tires, 24,115 thermo-fluid analysis of disk brakes, 119 flow-driven string dynamics in turbomachinery, [120][121][122] and flow analysis of turbocharger turbines. [123][124][125][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, 123-126 flow-driven string dynamics in turbomachinery, 121,122 ram-air parachutes, 132 spacecraft parachutes, 134 aortas, [111][112][113] heart valves, 116,117,111,113 tires with road contact and deformation, [128][129][130] and fluid films. 131,130 Using IGA basis functions in space is now a key part of some of the newest arterial zero-stress-state estimation methods [138][139][140][141]113 and related shell analysis. 142

Stabilization parameters and element lengths
In all the stabilized and VMS methods discussed in the earlier subsections, an embedded stabilization parameter, known as "τ ", plays a significant role (see Ref. 28). 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 τ s were proposed, starting with those in Refs. 1-5, followed by the the ones introduced in Refs. 13 and 14. In many cases, the element length was seen as an advection length scale, in the flow-velocity direction. The set of τ s introduced in Refs. 3-5 in conjunction with (SUPG) 82 is now called "τ 82 ". The τ definition introduced in Ref. 14, which is for the advective limit and is now called "τ SUGN1 " (and the corresponding element length is now called "h UGN "), automatically yields lower values for higher-order finite element basis functions (see Refs. 143 and 144). Later, other τ definitions that are applicable to higher-order elements were proposed in Ref. 145 in the context of advective-diffusive systems. The τ used in Ref. 16 with (SUPG) 82 was a slightly modified version of τ 82 . Subsequent minor modifications of τ 82 took into account the interaction between the DC and (SUPG) 82 terms in a fashion similar to how it was done in Ref. 14 for the advection-diffusion-reaction equation. Until 2004, all these slightly modified versions of τ 82 were always used with the same DC parameter, which was introduced in the 1991 ASME paper 16 and is now called "δ 91 ". This DC parameter was derived from the one given in Ref. 15 for the entropy variables.
Calculating the τ s based on the element-level matrices and vectors was introduced in Ref. 146 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 146,147 to address the difficulties reported at small time-step sizes. A second element length scale, in the solution-gradient direction and called "h RGN ", was introduced in 2001 Refs. 20 and 148. Recognizing this as a diffusion length scale, a new stabilization parameter for the diffusive limit, "τ SUGN3 ", was introduced in Refs. 20 and 149, to be used together with τ SUGN1 and "τ SUGN2 ", the parameters for the advective and transient limits. For the stabilized ST methods, "τ SUGN12 ", representing both the advective and transient limits, was also introduced in Ref. 20.
New ways of calculating the τ and DC parameter to be used with (SUPG) 82 were introduced in 2004 Refs. 149, 150 and 151. The new τ s, now categorized under the label "τ 04 ", have a matrix structure for viscous flows and reduce to a scalar for inviscid flows. The new DC parameters were of two types: one defined in a style the A node-numbering-invariant directional length scale for simplex elements 2725 DC Directional Dissipation 20,151,152 parameter was defined, and one that is now called "YZβ" DC parameter. The YZβ DC parameter is residual based, and it is simpler than δ 91 . It has options for smoother or sharper computed shocks. A number of 2D and 3D test computations with YZβ DC were reported in Refs. 153-155. These computations showed that in addition to being simpler than δ 91 , the YZβ DC parameter was superior in accuracy. The computations reported in Refs. [153][154][155] were based on the compressible-flow ST SUPG.
Some new options for the stabilization parameters used with the SUPS and VMS were proposed in Refs. 21, 24, 94, 95 and 99. These include a fourth τ component, "τ SUGN4 ", 24 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 24 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 Refs. 72, 156-166 and 167.

Directional element lengths for isogeometric discretization
The stabilization and DC parameters and element lengths discussed in the previous subsection were all originally intended for finite element discretization, but quite often used also for isogeometric discretization. The element lengths and stabilization and DC parameters introduced in Ref. 168 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 Ref. 168 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.

Directional element lengths for simplex elements
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 directional element lengths, to have node-numbering invariance for all element types, including simplex elements. The directional element length expression we introduce here meets that requirement. We analytically evaluate that element length expression in the context of simplex elements and compared to one of the most widely used element length expressions and show the levels of noninvariance avoided.

Outline of the remaining sections
In Sec. 2, we briefly describe the directional element length definition introduced in Ref. 168, which we will start with. In Sec. 3, we describe the parametric spaces we use in calculating the directional element length for simplex elements: the integration parametric space, which is not node-numbering-invariant, and the preferred parametric space, which is. In Sec. 4, we show the level of noninvariance we get with the integration parametric space, and we compare the element lengths based on the integration and preferred parametric spaces. We provide our concluding remarks in Sec. 5.

Overview of the Starting Element Length Definition
We briefly describe the element length definition introduced in Ref. 168, which we will start with. Given x ∈ R n sd and ξ ξ ξ ∈ R n pd , where n sd and n pd are the spatial and parametric dimensions, We define the Jacobian tensor as and the matrix of its components for given basis sets e i (i = 1, . . . , n sd ) and e α (α = 1, . . . , n pd ) as We introduce a transformation tensor D, and based on that we define the following tensor:Q With that, the directional element-length is defined aŝ where r is the unit vector representing the direction, defined in the physical space.
Remark 1. The factor 2 comes from a typical parametric space, which is of square or cubic shape with length 2. Other shapes and different lengths can be accounted for with D.
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. Some of the most widely used element length definitions do not have this invariance. The transformation tensor D was introduced in Ref. 168 to avoid such invariance. It was identified as a "scaling" tensor in Ref. 168, but considering its nondiagonal nature, we now identify it as a "transformation" tensor.
Suppose, we have a preferred parametric spaceξ ξ ξ for calculating the element length. Then,Q = ∂x and We use a different basis set,ê α , for the preferred parametric space.

Simplex Element Parametric Spaces
We consider n pd -dimensional simplex elements (n pd ≥ 2).

Integration parametric space
A typical integration parametric space is defined as ξ ξ ξ ∈ R n pd | ξ α ≥ 0 for α = 1, . . . , n pd and The βth component of the first n pd vertex-position vectors is and the last-vertex position is ξ ξ ξ n pd +1 = 0. Based on these, we define a set of linear shape functions: for α = 1, . . . n pd , and  The derivatives are for β = 1, . . . , n pd , and where (3.8) Figure 1 shows a 2D example.

Preferred parametric space
The integration parametric space, as described in Sec. 3.1, is not ideal for calculating the stabilization parameters. For example, it does not have node-numbering invariance for the directional element length definition given by Eq. (2.6) without D.
As the preferred parametric space, which we want to be free from the nodenumbering order, we select a regular simplex. Here, the regular simplex has vertex distances of 2. The first n pd vertex-position vectors arê To make it a regular simplex, the last-vertex location will be determined following the steps below. We start with defining a centroid based on the first n pd vertices:
which gives 20) and the solution is Because d > 0, we select d = 1 + n pd − 1.

Fig. 2. A triangular element suitable for directional element length calculation.
From the vertex-position vectors, we obtain While the tensor D depends on the orientation of the parametric spaceξ ξ ξ, the following product does not: From that, For n pd = 2, we obtain  We define two vectors: They are the vectors we get when we bring the vector r from the physical space to the two parametric spaces. With that, we havê From Eqs. (4.6) and (4.8), we obtain where • 2 is the spectral norm. From Eqs. (4.7) and (4.9), we obtain Combining Eqs. (4.10) and (4.11), we obtain 12) while noting that we will always have D −1 −1 2 ≤ D 2 . For the Cartesian basis sets, the matrix form of that is (4.13)

Dependence on the node-numbering order
Without D, the element length depends on the node-numbering order because the last vertex n pd + 1 is at an irregular location, and that is the only reason. Therefore, the element length depends only on which node corresponds to the last vertex in A node-numbering-invariant directional length scale for simplex elements 2733 the parametric space. With the affine transformatioñ ξ ξ ξ = R · ξ ξ ξ + b, (4.14) where we can have all the configurations by applying this transformation 1 to n pd times. Applying it n pd + 1 times takes us back to the configuration prior to any transformation. The element length with R is We compare that to the length without R. From Eq. (4.12), First, we prove that R k can be expressed as where m(a) = 1 + ((a − 1) mod (n pd + 1)). We note that, to deal with the algebra in the proof in a more systematic way, we have extended the summation range to n pd + 1, which has no actual effect because ξ ξ ξ n pd +1 = 0.
Proof. Let us define (4.20) Obviously, We now show that Actually, the expressions δ m(a) m(b+k) and δ m(a) m(k) cannot represent ξ ξ ξ n pd +1 · ξ ξ ξ n pd +1 , however, for that a, ξ ξ ξ m(a+1) − ξ ξ ξ 1 = 0. From Eqs. (4.21) and (4.22), We note that To evaluate the spectral norm of R k , we first evaluate (R k ) T · (R k ): − ξ ξ ξ m(b+k) · ξ ξ ξ m(k) + ξ ξ ξ m(k) · ξ ξ ξ m(k) )ξ ξ ξ a ξ ξ ξ b Here, I is the unit tensor, and 1 = n pd α=1 e α . We preclude m(k) = n pd + 1 because that would take us back to the starting parametric space. With that, (4.37) In the matrix version of that, all entries are 1, excluding the diagonal entries, where the value is 2, but not excluding the m(−k)th column. For example, for k = n pd , we get From this example, we can see that the eigenvalues of (R R R k ) T R R R k are independent of the k value, which implies R R R 2 = R R R 2 2 = · · · = R R R n pd 2 . Therefore it will be sufficient to inspect the eigenvalues of (R R R n pd ) T R R R n pd . From Eq. (4.31), R R R n pd = R R R −1 . With that, for notation convenience, we use the characteristic polynomial from det(λI I I − R R R −T R R R −1 ), where I I I is the unit matrix. In deriving that polynomial, we use the special method given in Appendix A, specifically we use c n (Eq. (A.22)) with p = λ − 2, q = −1 and r = λ − 1, obtaining From that, the eigenvalues are and for n pd ≥ 3, λ = 1 is also an eigenvalue with multiplicity of n pd − 2. The maximum eigenvalue is The eigenvector corresponding to λ max , v = v α e α , will be calculated using the matrix and it can be verified that n pd − 1 + (n pd − 1)(n pd + 3) 2 (otherwise).   For n pd = 2, and for n pd = 4, The range expands at both ends with increasing parametric-space dimension.

Comparison of the element lengths based on the integration and preferred parametric spaces
Here, we compare h andĥ. From Eq. (3.32), The maximum eigenvalue is λ max = 2(1 + n pd ), and D 2 = 2(1 + n pd ). (4.59) The eigenvector corresponding to λ max , v = v α e α , will be calculated using the matrix and in that direction,  and for n pd = 4, The range expands at the upper end with increasing parametric-space dimension.

Concluding Remarks
Stabilized and VMS methods are widely used in flow computations. Stabilization parameters embedded in most of these methods almost always involve element lengths, most of the time in specific directions, such as the direction of the flow or solution gradient. We require the element lengths, including the directional element lengths, to have node-numbering invariance for all element types, including simplex elements. We have introduced a directional element length expression meeting that requirement. We accomplished that by using in the element length calculations for simplex elements a preferred parametric space instead of the standard integration parametric space. In the context of simplex elements, we analytically evaluated the element length expressions based on the two parametric spaces. We showed 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. We also showed 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.

A.1. A A A n and B B B n
The first two special matrices are We define a n = det A A A n and b n = det B B B n . To obtain direct expressions for these, we look for recurrence relationships: q q · · · q p 0 p q · · · q q . . . . . . . . . . . .
Note that we used the following relationship:  with a 1 = p and b 1 = q. We eliminate a n and obtain qa n+1 = pb n+1 + nq(p − q)b n . We note that this expression represents b n even for p − q = 0 by letting (p − q) 0 = 1. We substitute this into Eq. (A.6) and obtain (p − q) n q = qa n − nq(p − q) n−1 q. (A. 16) Dividing the equation by q under the condition q = 0, we obtain a n = (p − q) n−1 (p + (n − 1)q). (A.17) We note that even when q = 0, the above equation represents the obvious solution a n = p n , and it generally holds, again by letting (p − q) 0 = 1.

A.2. C C C n
The determinant calculation for the third special matrix will rely on the determinant of A A A n and B B B n , and therefore it is more convenient to provide its definition for n + 1: r q · · · q q q p q . . . . . . . . . = r(p − q) n−1 (p + (n − 1)q) − qn(p − q) n−1 q (A.20) = (p − q) n−1 (r(p + (n − 1)q) − q 2 n).