**2. Temporal splitting scheme**

and **u** ⊗ **v** ¼ *uiv <sup>j</sup>*, *i*, *j* ¼ 1, 2*:* Indeed, it holds

*Vortex Dynamics Theories and Applications*

N ð Þ¼ **u**

interior edges of <sup>E</sup>*<sup>h</sup>* and <sup>Γ</sup>*<sup>B</sup>*

and denote its jump by

and *Mh*ð Þ ⊂ *M* for integer *N* ≥1:

**84**

*Wr*,*<sup>s</sup>*

*<sup>V</sup>* ¼ f**v**<sup>∈</sup> *<sup>L</sup>*<sup>2</sup>

*<sup>M</sup>* ¼ f*q*∈*L*<sup>2</sup>

*<sup>∂</sup> <sup>u</sup>*<sup>2</sup> ð Þ *∂x* þ *<sup>∂</sup>*ð Þ *uv ∂y*

1

*<sup>h</sup>* be the set of all boundary edges. Set <sup>Γ</sup>*<sup>h</sup>* <sup>¼</sup> <sup>Γ</sup>*<sup>I</sup>*

ð Þ *<sup>E</sup>* : <sup>∀</sup>j*m*j<sup>≤</sup> *<sup>r</sup>*, *<sup>∂</sup>mv*<sup>∈</sup> *Ls* f g ð Þ *<sup>E</sup> :*

� �

> � �

*<sup>E</sup>* <sup>∈</sup> *<sup>W</sup>*2,4*=*<sup>3</sup>

*<sup>E</sup>* <sup>∈</sup>*W*1,4*=*<sup>3</sup>

*:*

ð Þ *E* � �<sup>2</sup>

ð Þg *E :*

*i*,*j TijSij*.

,

g,

A locally conservative DG discretization will be employed for the Navier–Stokes Eq. (9)–(11). We denote by E*<sup>h</sup>* a shape-regular triangulation of the domain Ω into

any nonnegative integer *r* and *s*≥1, the classical Sobolev space on a domain *E* ⊂ <sup>2</sup> is

ð Þ <sup>Ω</sup> <sup>2</sup> : <sup>∀</sup>*E*∈E*h*, **<sup>v</sup>**

The jump and average of a function *ϕ* on an edge *e* are defined by:

½ �¼ *<sup>ϕ</sup> <sup>ϕ</sup> Ek* <sup>Þ</sup> �

<sup>2</sup> *<sup>ϕ</sup>Ek* � �� � *<sup>e</sup>* þ *ϕ El* Þ � � � � *e*

f g*<sup>ϕ</sup>* <sup>¼</sup> <sup>1</sup>

tensor product of two tensors **<sup>T</sup>** and **<sup>S</sup>** is defined as **<sup>T</sup>** : **<sup>S</sup>** <sup>¼</sup> <sup>P</sup>

**<sup>V</sup>***<sup>h</sup>* <sup>¼</sup> **<sup>v</sup>**∈*L*<sup>2</sup>

*Mh* <sup>¼</sup> *<sup>q</sup>*∈*L*<sup>2</sup>

ð Þ Ω : ∀*E* ∈E*h*, *q*

� � � *<sup>e</sup>* � *ϕ El* Þ � � � � *e* ,

� �

� � � �

Further, let **v** be a piecewise smooth vector-, or matrix-valued function at **x**∈*e*

½ � **v** ≔ **v**<sup>þ</sup> � **n***<sup>E</sup>*<sup>þ</sup> þ **v**� � **n***<sup>E</sup>*� ,

where *e* is shared by two elements *E*<sup>þ</sup> and *E*�, and an outward unit normal vector **n***<sup>E</sup>*<sup>þ</sup> (or **n***<sup>E</sup>*� ) is associated with the edge *e* of an element *E*<sup>þ</sup> (or *E*�). The

Let *N*ð Þ *E* be the set of polynomials on an element *E* with degree no more than *N:* Based on the triangulation, we introduce two approximate subspaces **V***h*ð Þ ⊂*V*

ð Þ <sup>Ω</sup> <sup>2</sup> : <sup>∀</sup>*E*<sup>∈</sup> <sup>E</sup>*h*, **<sup>v</sup>***<sup>h</sup>* <sup>∈</sup> ð Þ *N*ð Þ *<sup>E</sup>* <sup>2</sup> n o

ð Þ <sup>Ω</sup> : <sup>∀</sup>*E*<sup>∈</sup> <sup>E</sup>*h*, *<sup>q</sup>* <sup>∈</sup>*<sup>N</sup>*�<sup>1</sup>ð Þ *<sup>E</sup>* � �*:*

We mainly cite the content of [17], in which was motivated by the work of Girault, Rivière and Wheeler in a series of papers [2, 14]. Some projection methods [6, 18] have been developed to overcome the incompressibility constraints ∇ � **u** ¼ 0. An implementation of the operator-splitting idea for discontinuous Galerkin elements was developed in [2]. We appreciate the advantages of the discontinuous Galerkin methods, such as local mass conservation, high order of approximation, robustness and stability. In this work, we will make use of the underlying physical nature of incompressible flows in the literature and extend the interior penalty discontinuous

CCCA <sup>¼</sup> <sup>∇</sup> � *<sup>F</sup>*ð Þ **<sup>u</sup>** *:*

*<sup>h</sup>* be the set of all

*<sup>h</sup>* ∪ Γ*<sup>B</sup> <sup>h</sup> :* For

0

BBB@

triangles, where *h* is the maximum diameter of elements. Let Γ*<sup>I</sup>*

ð Þ¼ *<sup>E</sup> <sup>v</sup>*<sup>∈</sup> *Ls*

We define the spaces of discontinuous functions

*<sup>∂</sup>*ð Þ *uv ∂x* þ *<sup>∂</sup> <sup>v</sup>*<sup>2</sup> ð Þ *∂y*

> We consider here a third-order time-accurate discretization method at each time step by using the previous known velocity vectors. Let *<sup>Δ</sup><sup>t</sup>* be the time step, *<sup>M</sup>* <sup>¼</sup> *<sup>T</sup> Δt* , and *tn* ¼ *nΔt*. The semi-discrete forms of problem (6)–(8) at time *tn*þ<sup>1</sup> is

$$\frac{\gamma\_0 \mathbf{u}^{n+1} - a\_0 \mathbf{u}^n - a\_1 \mathbf{u}^{n-1} - a\_2 \mathbf{u}^{n-2}}{\Delta t} - \nu \Delta \mathbf{u}^{n+1} + \nabla p^{n+1} \tag{13}$$

$$\mathbf{f} = -\beta\_0 \mathcal{N}(\mathbf{u}^n) - \beta\_1 \mathcal{N}\left(\mathbf{u}^{n-1}\right) - \beta\_2 \mathcal{N}\left(\mathbf{u}^{n-2}\right) + \mathbf{f}(t\_{n+1}),$$

$$\nabla \cdot \mathbf{u}^{n+1} = \mathbf{0}, \tag{14}$$

which has a timestep constraint based on the CFL condition (see [19]):

$$
\Delta t \approx O\left(\frac{\mathcal{L}}{l \mathcal{M}^2}\right),
$$

where L is an integral length scale (e.g. the mesh size) and U is a characteristic velocity. Because the semi-discrete system (13)–(14) is linearized, thus, a timesplitting scheme can be applied naturally, i.e., the semi-discretization in time (13)– (14) can be decomposed into three stages as follows.

• The first stage

When **u***<sup>n</sup>* and **u***<sup>n</sup>*�<sup>1</sup> (*n*≥ 1) are known, the following linearized third-order formula can be used

$$\frac{\gamma\_0 \tilde{\mathbf{u}} - a\_0 \mathbf{u}^n - a\_1 \mathbf{u}^{n-1} - a\_2 \mathbf{u}^{n-2}}{\Delta t} = -\beta\_0 \mathcal{N}(\mathbf{u}^n) - \beta\_1 \mathcal{N}(\mathbf{u}^{n-1}) + \mathbf{f}(t\_{n+1}) \tag{15}$$

with the following coefficients for the subsequent time levels (*n*≥2)

$$y\_0 = \frac{11}{6}, \quad a\_0 = 3,\\ a\_1 = -\frac{3}{2}, \quad a\_2 = \frac{1}{3}, \quad \beta\_0 = 2, \quad \beta\_1 = -1.$$

Especially, by using the Euler forward discretization at the first time step (*n* ¼ 0), we can get a medium velocity field **u**<sup>1</sup> by

$$\frac{\mathbf{u}^1 - \mathbf{u}^0}{\Delta t} = -\mathcal{N}(\mathbf{u}^0) + \mathbf{f}(t\_1),$$

and **u**<sup>2</sup> by

$$\frac{\gamma\_0 \mathbf{u}^2 - a\_0 \mathbf{u}^1 - a\_1 \mathbf{u}^0}{\Delta t} = -\beta\_0 \mathcal{N} \left( \mathbf{u}^1 \right) - \beta\_1 \mathcal{N} \left( \mathbf{u}^0 \right) + \mathbf{f}(t\_2),$$

which adopts the following coefficients to construct a second-order difference scheme for the time level (*n* ¼ 2) in (15)

$$y\_0 = \frac{3}{2}, \quad a\_0 = 2,\\ a\_1 = -\frac{1}{2}, \quad a\_2 = 0, \ \beta\_0 = 2, \ \beta\_1 = -1.$$

Note that the coefficients in (15) are adjustable, but high-order time discrete schemes need to be verified with stability analysis.

• The second stage

The pressure projection is as follows

$$
\gamma\_0 \frac{\tilde{\tilde{\mathbf{u}}} - \tilde{\mathbf{u}}}{\Delta t} = -\nabla p^{n+1}.\tag{16}
$$

are obviously a type of elliptic and Helmholtz problems at each time step as *n*≥2. We decouple the incompressibility condition and the nonlinearity, then the pressure and velocity semi-discretizations (17)–(18) will be formulated by the interior penalty discontinuous Galerkin methods in spacial discretizations in the

*A Fully Discrete SIPG Method for Solving Two Classes of Vortex Dominated Flows*

For spacial approximations, assume that piecewise polynomials of order *N* are

the approximating polynomial space for the velocity or pressure restricted to each element, a high-order nodal basis can be chosen, consisting of Lagrange interpolating polynomials defined on a reference simplex introduced in [20, 21]. We let **u** be approximated by **u***<sup>h</sup>* ∈ **V***<sup>h</sup>* and adopt a suitable approximation for the term *F*, i.e.,

each element of T *<sup>h</sup>*. Multiplying the nonlinear term by a test function **v***<sup>h</sup>* ∈ **V***h*, integrating over the computational domain, and applying integration by parts,

ð Þ� *<sup>F</sup>* � <sup>∇</sup> **<sup>v</sup>***hdx* <sup>þ</sup> <sup>X</sup>

*∂vhi ∂x <sup>j</sup>*

, respectively, where **u***<sup>E</sup>*<sup>1</sup> and **u***<sup>E</sup>*<sup>2</sup> are the mean values of *uh* on the

correspond to the components of the related vectors. On each edge *e*∈*∂E*<sup>1</sup> ∩ *∂E*<sup>2</sup> shared by two elements, to ensure the flux Jacobian of purely real eigenvalues, we

*<sup>E</sup>*2,*<sup>e</sup>* the largest eigenvalue of the Jacobians *<sup>∂</sup>*

elements *E*<sup>1</sup> and *E*2, respectively. The global Lax-Friedrichs flux is generally more dissipative than the local Lax-Friedrichs flux, therefore, we primarily consider the local flux on each edge. Although the Lax-Friedrichs flux is perhaps the simplest numerical flux and often the most efficient flux, it is not the most accurate scheme. A remedy of the problem is to employ high-order finite elements. By replacing the

**n***<sup>e</sup>* � ½ �¼ *F* � **v***<sup>h</sup>* **n***<sup>e</sup>* � f g*F* � ½ �þ **v***<sup>h</sup>*

*e*∈Γ*<sup>h</sup>*

*λe*

� �, one can get a DG discretization for the nonlinear

For the pressure correction step (17) and the viscous correction step (18), we use the SIPG method to approximate the correction steps. Choosing the orthonormal Legendre basis and the Legendre-Gauss-Lobatto quadrature points gives a wellconditioned Vandermonde matrix and the resulting interpolation well behaved, which greatly simplifies the formulas. The *C*<sup>0</sup> continuity condition of the basis in the discontinuous Galerkin formulation is not required. Enforcing a weak continu-

� � <sup>¼</sup> *Lp <sup>ϕ</sup><sup>h</sup>* ð Þ, *<sup>ϕ</sup>*<sup>∈</sup> *Mh*,

<sup>2</sup> ½ �� **<sup>u</sup>***<sup>h</sup>* ½ � **<sup>v</sup>***<sup>h</sup>* ,

ð *e* *<sup>k</sup>*¼<sup>1</sup>*N*ð Þ *Ek*


**n***<sup>e</sup>* � ½ � *F* � **v***<sup>h</sup> ds*, (19)

� � **u***E*<sup>1</sup> and

, for *i*, *j* ¼ 1, 2 and the indexes *i*, *j*

*<sup>∂</sup>***<sup>u</sup>** ð Þ *F* � **n***<sup>e</sup>*

2 . In

employed, then the approximation space can be rewritten as **<sup>V</sup>***<sup>h</sup>* <sup>¼</sup> <sup>⊕</sup>*<sup>K</sup>*

*<sup>F</sup>*ð Þ **<sup>u</sup>** <sup>≈</sup> *<sup>F</sup>*ð Þ **<sup>u</sup>***<sup>h</sup>* , where *<sup>F</sup>*ð Þ **<sup>u</sup>***<sup>h</sup>* also can be represented as the *<sup>L</sup>*<sup>2</sup>

*Ek*

ð *Ek*

next section.

we have

ð Ω

may define *λ*<sup>þ</sup>

� � **u***E*<sup>2</sup>

*∂ <sup>∂</sup>***<sup>u</sup>** ð Þ *F* � **n***<sup>e</sup>*

**87**

**3. The spatial discretizations**

*DOI: http://dx.doi.org/10.5772/intechopen.94316*

ð Þ� <sup>∇</sup> � *<sup>F</sup>* **<sup>v</sup>***hd***<sup>x</sup>** ¼ �<sup>X</sup>

*<sup>E</sup>*1,*<sup>e</sup>*, *λ*�

integrand in the surface integral as

*<sup>E</sup>*1,*<sup>e</sup>*, *λ*� *E*2,*e*

ity on the interior edges by a penalty term, we have for (17)

*a pn*þ<sup>1</sup> *<sup>h</sup>* , *ϕ<sup>h</sup>*

term in (19) by the local Lax-Friedrichs flux.

with *λ<sup>e</sup>* ¼ max *λ*<sup>þ</sup>

where the term ð Þ� *F* � ∇ **v***<sup>h</sup>* equals to *Fij*

To seek *pn*þ<sup>1</sup> such that <sup>∇</sup> � **<sup>u</sup>**~~ <sup>¼</sup> 0, we solve the system

$$-\Delta p^{n+1} = -\frac{\chi\_0}{\Delta t} \nabla \cdot \bar{\mathbf{u}},\tag{17}$$

with a Neumann boundary condition being implemented on the boundaries as

$$\begin{split} \frac{\partial p^{n+1}}{\partial \mathbf{n}} &= \mathbf{f}^{n+1} - \beta\_0 \mathbf{n} \cdot \left( \frac{\partial \mathbf{u}^n}{\partial t} + \nabla \cdot F(\mathbf{u}^n) - \nu \Delta \mathbf{u}^n \right) \\ &- \beta\_1 \mathbf{n} \cdot \left( \frac{\partial \mathbf{u}^{n-1}}{\partial t} + \nabla \cdot F(\mathbf{u}^{n-1}) - \nu \Delta \mathbf{u}^{n-1} \right) \\ &- \beta\_2 \mathbf{n} \cdot \left( \frac{\partial \mathbf{u}^{n-2}}{\partial t} + \nabla \cdot F(\mathbf{u}^{n-2}) - \nu \Delta \mathbf{u}^{n-2} \right) \coloneqq G\_n. \end{split}$$

One can compute the vorticity *<sup>ω</sup><sup>n</sup>* <sup>¼</sup> <sup>∇</sup> � **<sup>u</sup>***<sup>n</sup>* at time *<sup>t</sup> <sup>n</sup>* <sup>¼</sup> *<sup>n</sup>* � *<sup>Δ</sup>t*. Then we use *<sup>p</sup><sup>n</sup>*þ<sup>1</sup> to update the intermediate velocity **u**~~ by (16).

• The third stage is completed by solving

$$
\gamma\_0 \frac{\mathbf{u}^{n+1} - \tilde{\tilde{\mathbf{u}}}}{\Delta t} = \nu \Delta \mathbf{u}^{n+1},
$$

which can be written as a Helmholtz equation for the velocity

$$-\Delta \mathbf{u}^{n+1} + \frac{\mathcal{Y}\_0}{\nu \Delta t} \mathbf{u}^{n+1} = \frac{\mathcal{Y}\_0}{\nu \Delta t} \tilde{\mathbf{u}}.\tag{18}$$

From the three stages given above, we notice that (15) in the semi-discrete systems is presented in a linearized and explicit process, moreover, (17) and (18) are obviously a type of elliptic and Helmholtz problems at each time step as *n*≥2. We decouple the incompressibility condition and the nonlinearity, then the pressure and velocity semi-discretizations (17)–(18) will be formulated by the interior penalty discontinuous Galerkin methods in spacial discretizations in the next section.
