A Grid-Insensitive LDA Method on Triangular Grids Solving the System of Euler Equations

Abstract

The performance of the classic upwind-type residual distribution (RD) methods on skewed triangular grids are rigorously investigated in this paper. Based on an improved signals distribution, an improved second order RD method based on the LDA approach is proposed to faithfully replicate the flow physics on skewed triangular grids. It will be mathematically and numerically shown that the improved LDA method is found to have minimal accuracy variations when grids are skewed compared to classic RD and cell vertex finite volume methods on scalar equations and system of Euler equations.

Appendices

Appendix 1: General Equation for First Order Finite Volume

According to Fig. 32, the scaled normal vector for each edge and the upwind first order value for the edge are shown in Table 4. Therefore, the line integration of the first order finite volume will be,

$$\begin{aligned} \begin{array}{rl} \oint \varvec{F}\cdot \hat{n}dS= &{}u_0\left( \frac{ak}{2}+\frac{bh}{6}+\frac{ak}{2}+ \frac{bh}{3}+\frac{bh}{3}+\frac{bh}{6}+ \frac{ak}{2}-\frac{bh}{6}+\frac{ak}{2}-\frac{bh}{6}\right) \\ &{}+u_3\left( -\frac{ak}{2}+\frac{bh}{6}-\frac{ak}{2}+ \frac{bh}{6}\right) +u_4\left( -ak-\frac{bh}{3}\right) +u_5\left( -\frac{2bh}{3}\right) \end{array} \end{aligned}$$

(66)

Thus,

$$\begin{aligned} \text {TE}=\frac{3 a k (2 u_0-u_3-u_4)+b h (2 u_0+u_3-u_4-2 u_5)}{6 h k} \end{aligned}$$

(67)

Table 4 The scaled normal of the cell vertex finite volume cell; and, the upwind edge value ($\frac{b}{a}<\frac{k}{h}$)

Fig. 32

The isosceles grid topology for two kinds (finite volume cells are shaded)

Appendix 2: Finite Volume Method Characteristics

The second order finite volume which is used in this paper is based on the following characteristics:

• The node center (cell vertex) is used in order to maintain same computational points as residual distribution approach.

• The gradients are calculated based on Least Square approach for each point [17].

• The solver is also same as the residual distribution to ensure a true comparison. Hence, the first order explicit is used [17].

• The local time step is implemented based on the specific local time step [24].

• The upwind discretization for the scalar equation is constructed by Roe’s flux. Also, for the Euler equations in this study we use the Roe’s flux [25].

Appendix 3: Lax–Friedrichs

In the Lax–Friedrichs method the sub-residuals should be constructed to ensure the first order accuracy and also $L_\infty$ stability.

$$\begin{aligned} \phi _i^\tau =\bar{\phi }^\tau +\alpha \left( u_i-\bar{u}\right) ,\quad \alpha >\max ||K_i|| \end{aligned}$$

(68)

These choices guaranty that the scheme is $L_\infty$ stable [22]. Note that, the $\bar{\phi }^\tau$ and $\bar{u}$ are the arithmetic average of $\phi ^\tau _i$ and $u_i$ for a cell.

Appendix 4: Non-positivity of Weighted N scheme

Lemma 7.1

The weighted N-scheme approach does NOT satisfy positivity for a skewed element.

Proof

Using Eq. 25 for the N-scheme to construct weighted N leads us to,

$$\begin{aligned} \begin{array}{l} \phi ^\text {N(W)}_i= \left( \frac{1}{3}- \frac{\vec {N}_i\cdot \vec {N}_i}{3\min \left( \vec {N}_i\cdot \vec {N}_p\right) }\right) \phi _i^\text {N} +\left( \frac{1}{3}- \frac{\vec {N}_j\cdot \vec {N}_i}{3\min \left( \vec {N}_j \cdot \vec {N}_p\right) }\right) \phi _j^\text {N}\\ \qquad \qquad \,\, +\left( \frac{1}{3}- \frac{\vec {N}_k\cdot \vec {N}_i}{3\min \left( \vec {N}_k\cdot \vec {N}_p\right) }\right) \phi _k^\text {N}\\ \phi ^\tau _j= \left( \frac{1}{3}- \frac{\vec {N}_i\cdot \vec {N}_j}{3\min \left( \vec {N}_i \cdot \vec {N}_p\right) }\right) \phi _i^\text {N} +\left( \frac{1}{3}- \frac{\vec {N}_j\cdot \vec {N}_j}{3\min \left( \vec {N}_j \cdot \vec {N}_p\right) }\right) \phi _j^\text {N}\\ \quad \qquad +\left( \frac{1}{3}- \frac{\vec {N}_k\cdot \vec {N}_j}{3\min \left( \vec {N}_k \cdot \vec {N}_p\right) }\right) \phi _k^\text {N}\\ \phi ^\tau _k= \left( \frac{1}{3}- \frac{\vec {N}_i\cdot \vec {N}_k}{3\min \left( \vec {N}_i \cdot \vec {N}_p\right) }\right) \phi _i^\text {N} +\left( \frac{1}{3}- \frac{\vec {N}_j\cdot \vec {N}_k}{3\min \left( \vec {N}_j \cdot \vec {N}_p\right) }\right) \phi _j^\text {N}\\ \quad \qquad +\left( \frac{1}{3}- \frac{\vec {N}_k\cdot \vec {N}_k}{3\min \left( \vec {N}_k \cdot \vec {N}_p\right) }\right) \phi _k^\text {N} \end{array} \end{aligned}$$

(69)

where,

$$\begin{aligned} \phi _{i,j,k}^\text {N}=k_{i,j,k}^+(u_{i,j,k}-\hat{u}),\quad \hat{u}=\frac{\sum _pk_p^-u_p}{\sum _pk_p^-} \end{aligned}$$

(70)

to satisfy local positivity for a signal ($\phi _i$) all the coefficients of the main nodes $u_i$ should be positive and all the coefficients for other points within the same element should be negative. For instance in $\phi _i$ all the coefficients of $u_i$ must be positive. The other coefficients need to be negative. While, $K_j^+,K_k^+$ are positive; also, all the values of w are positive; the positivity is lost according to the second and third term. Therefore, we shall abandon the weighted N-scheme approach and will not discuss it further in the paper. However, the ’weighted’ idea will be applied on the LDA. $\square$

Appendix 5: Grid Topology

1.1 Isotropic Grid

The isotropic grid is an unstructured-type grid with right triangular elements. In order to control the skewness one could use the ratio of two right angle edges as shown in Fig. 33. The relation between the stretching parameter ($s=\frac{k}{h}$) and the skewness is shown in the next equation, although its derivation is omitted here for brevity. The details can be found in [12].

$$\begin{aligned} Q=\frac{2}{\pi }\arctan \left( \frac{s}{2}\right) \end{aligned}$$

(71)

Note that the minimum value for s will be 1. Substituting $s=1$, the minimum value for the skewness will be $Q\simeq 0.3$ which is for right triangles. In other words, the best condition for isotropic grid would be having skewness of 0.3.

The grid skewness will be determined by setting the stretching parameter. In this grid type there are two kinds of points, one with eight neighboring elements and one with four neighbors. Note that, for the finite volume cell vertex the median points of the neighboring cells and midpoints of the edges are used to demonstrate the cell. Thus, the outer of a finite volume cell is depended on the arrangement of the neighboring cells. For the first and second type of isotropic grids, the finite volume cell vertex element is shaded in Fig. 33. It is much easier to handle the isotropic grids in terms of achieving a uniformly skewed grids over the whole computational domain and also when doing the mathematical analyses. As such, the isotropic grids are used for scalar problems where rigorous mathematical analyses would be performed.

Fig. 33

The isotropic grid element and topology for RD and FV methods (finite volume cells are shaded)

1.2 Anisotropic Grid

The Delaunay triangulation is used to generate a fully unstructured or anisotropic grid over a cylinder. After generating the grid it would be randomized the in a way that different quality of the grids in terms of skewness could be built. This is where achieving a uniform grid skewness over the complete domain would not be possible but rather the grid skewness would have a range (or distribution). It should be mentioned that each randomization construct a different skewness distribution. Since the elements are not necessarily right triangular elements (unlike isotropic), the skewness ranges from 0 to 1.0.

Fig. 34

The randomize grid element area with radius of the minimum distance of the each point from the surrounding edges

According to Fig. 34, each point will move in fully randomize direction with a finite maximum distance (R) which avoids grid overlapping.

• The randomization percentage: The maximum distance that a point can move from its original place is R which we can be controlled in terms of percentage defined as $\alpha \times R$. A suitable value for $\alpha$ is chosen to implement grid irregularity. Larger values of $\alpha$ denote a higher percentage grid randomization.

• Randomization number: To build a much more realistic unstructured grid one could perform the whole process (n) times, to build even more randomized grid.

The two options above might be written as $(\alpha ,n)$. It should be mentioned that in this study, we are using three different combination of randomization to cover the possibilities in the engineering problems which are $(20\%,2)$, $(50\%,5)$ and $(90\%,9)$. For simplicity, we are calling these three randomization grids as G1, G2 and G3. The anisotropic grids will only be used when solving the system of Euler equations.

1.3 Showing Zero Weight for One Characteristic Projection

Lemma 7.2

At least one of the ratios in Eq. 27 for a specific edge number e is zero but the summation is always one.

Proof

The minimum value of $\vec {N}_i\cdot \vec {N}_p$ could be found by $p=i,j,k$. Consider $p=i$ then,

$$\begin{aligned} \frac{A_i}{A^\tau }=\frac{1}{3}-\frac{\vec {N}_i \cdot \vec {N}_i}{3\left( \vec {N}_i\cdot \vec {N}_i\right) }=0 \end{aligned}$$

(72)

If $p=j$, then,

$$\begin{aligned} \frac{A_j}{A^\tau }=\frac{1}{3}- \frac{\vec {N}_i\cdot \vec {N}_j}{3\left( \vec {N}_i\cdot \vec {N}_j\right) }=0 \end{aligned}$$

(73)

And, if $p=k$, then,

$$\begin{aligned} \frac{A_k}{A^\tau }=\frac{1}{3} -\frac{\vec {N}_i\cdot \vec {N}_k}{3\left( \vec {N}_i\cdot \vec {N}_k\right) }=0 \end{aligned}$$

(74)

Moreover, it is obvious that,

$$\begin{aligned} \sum _p\frac{A_p}{A^\tau }=\frac{A_i}{A^\tau }+\frac{A_j}{A^\tau }+\frac{A_k}{A^\tau }=1,\quad p=i,j,k. \end{aligned}$$

(75)

$\square$