**Abstract**

To simulate incompressible Navier–Stokes equation, a temporal splitting scheme in time and high-order symmetric interior penalty Galerkin (SIPG) method in space discretization are employed, while the local Lax-Friedrichs flux is applied in the discretization of the nonlinear term. Under a constraint of the Courant–Friedrichs– Lewy (CFL) condition, two benchmark problems in 2D are simulated by the fully discrete SIPG method. One is a lid-driven cavity flow and the other is a circular cylinder flow. For the former, we compute velocity field, pressure contour and vorticity contour. In the latter, while the von Kármán vortex street appears with Reynolds number 50≤ *Re* ≤400, we simulate different dynamical behavior of circular cylinder flows, and numerically estimate the Strouhal numbers comparable to the existing experimental results. The calculations on vortex dominated flows are carried out to investigate the potential application of the SIPG method.

**Keywords:** Navier–stokes equations, von Kármán vortex street, discontinuous Galerkin method, interior penalty

## **1. Introduction**

The Navier–Stokes equations are a concise physics model of low Knudsen number (i.e. non-rarefied) fluid dynamics. Phenomena described with the Navier– Stokes equations include boundary layers, shocks, flow separation, turbulence, and vortices, as well as integrated effects such as lift and drag. Analytical solutions of real flow problems including complex geometries are not available, therefore numerical solutions are necessary. The Navier–Stokes equation has been investigated by many scientists conducting research on numerical schemes for approximation solutions (see [1–8]).

There exist two ways to provide reference data for such problems: One consists in the measurement of quantities of interest in physical experiments and the other is to perform careful numerical studies with highly accurate discretizations. With the prevalence and development of high-performance computers, advanced numerical algorithms are able to be tested for the validation of approaches and codes and for high-order convergence behavior of delicate discretizations.

Among discontinuous Galerkin (DG) methods, primal schemes and mixed methods are distinguished. The former depend on appropriate penalty terms of the discontinuous shape functions, while the latter rely on the mixed methods as the

original second-order or higher-order partial differential equations are written as a system of first order partial differential equations with designed suitable numerical fluxes. Interior penalty discontinuous Galerkin methods (IPDG) are known as the representative of primal schemes whereas local discontinuous Galerkin methods (LDG) belong to the class of mixed methods [9]. About IPDG, symmetric interior penalty Galerkin (SIPG) and non-symmetric interior penalty Galerkin (NIPG) methods were first introduced originally for elliptic problems by Wheeler [10] and Rivière et al. [11]. Then, we presented some numerical analysis and simulations on nonlinear parabolic problems [12]. Recently, some work based on the SIPG and NIPG methods has been successfully applied to the steady-state and transient Navier–Stokes Equations [2, 13, 14], with careful analysis being conducted on optimal error estimates for the velocity.

The physics of Navier–Stokes flows are non-dimensionalized by Mach number *M* and Reynolds number *Re*,

$$M = \frac{\mathbf{u}\_{\Leftrightarrow}}{\mathbf{a}},\tag{1}$$

dynamicist Theodore Kármán (1963; 1994). Vortex streets are ubiquitous in nature and are visibly seen in river currents downstream of obstacles, atmospheric phenomena, and the clouds of Jupiter (e.g. The Great Red Spot). Shed vortices are also the primary driver for the zig-zag motion of bubbles in carbonated drinks. The bubble rising through the drink creates a wake of shed vorticity which impacts the integrated pressures causing side forces and thus side accelerations. The physics of sound generation with an Aeoleans harp operates by alternating vortices creating harmonic surfaces pressure variations leading to radiated acoustic tones. Tones generated by vortex shedding are the so-called Strouhal friction tones. If the diameter of the string, or cylinder immersed in the flow is *D* and the free stream velocity of the flow is **u**<sup>∞</sup> then the

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

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

*St* <sup>¼</sup> *fD* **u**<sup>∞</sup>

Czech physicist who experimented in 1878 with wires experiencing vortex shedding and singing in the wind (Strouhal, 1878; White, 1999). The Strouhal formula provides an experimentally derived shedding frequency for fluid flow. Therefore, we are interested in an investigation of Stouhal numbers of incompressible flow at

fluid flow in 2D is described by the Navier–Stokes equations, which include the equations of continuity and momentum, written in dimensionless form [8] as follows:

*<sup>α</sup>***<sup>u</sup>** <sup>þ</sup> ð Þ <sup>1</sup> � *<sup>α</sup> <sup>∂</sup>***<sup>u</sup>**

**u**j

*<sup>F</sup>*ð Þ¼ **<sup>u</sup> <sup>u</sup>** <sup>⊗</sup> **<sup>u</sup>** <sup>¼</sup> *<sup>u</sup>*<sup>2</sup> *uv*

Here the parameter *α* has the limit values of 0 for the free-slip (no stress) condition (Neumann) and 1 for the no-slip condition (Dirichlet); **u** ¼ ð Þ *u*, *v* is the velocity; *t* is the time; and *p* is the pressure. In general, the external force **f** is not

Using the divergence free constraint, problem (6)–(8) can be rewritten in the

, and *St* is the Strouhal number named after Vincent Strouhal, a

� *ν*Δ**u** þ ð Þ **u** � ∇ **u** þ ∇*p* ¼ **f**, in Ω � ð Þ 0, *T* (6)

*<sup>∂</sup>***<sup>n</sup>** <sup>¼</sup> **<sup>u</sup>**∞*:*

� *ν*Δ**u** þ ∇ � *F* þ ∇*p* ¼ **f**, in Ω � ð Þ 0, *T* (9)

*uv v*<sup>2</sup>

∇ � **u** ¼ 0, in Ω � ð Þ 0, *T* (10)

*<sup>t</sup>*¼<sup>0</sup> <sup>¼</sup> **<sup>u</sup>**0, (11)

*:* (12)

∇ � **u** ¼ 0, in Ω � ð Þ 0, *T* (7)

**u**j*<sup>t</sup>*¼<sup>0</sup> ¼ **u**0, (8)

, (5)

. The dynamics of an incompressible

shedding frequency *f* of the sound is given by the Strouhal formula

where *<sup>f</sup>* <sup>¼</sup> *<sup>T</sup>*�<sup>1</sup>

different Reynolds number.

*∂***u** *∂t*

taken into account in Eq. (6).

following conservative flux form [16]:

with the flux *F* being defined as

**83**

*∂***u** *∂t*

Let Ω be a bounded polygonal domain in <sup>2</sup>

subject to the boundary conditions on ∂Ω :

$$Re = \frac{\rho \mathbf{u}\_{\ast \bullet} D}{\mu},\tag{2}$$

where *ρ* is the density of the fluid, and *μ* is the dynamic viscosity. The kinematic viscosity *ν* is the ratio of *μ* to *ρ*. At low Knudsen numbers, Navier–Stokes surface boundary conditions are effectively no-slip (i.e. zero velocity). Diffusion of momentum from freestream to surface no-slip velocities forms boundary layers decreasing in thickness as Reynolds number increases. Thus, the range of characteristic solution scale increases as the Reynolds number increases. Nonlinear convective terms coupled with the strong velocity gradients in the Navier Stokes equations drive fluid flow at even moderate Reynolds numbers to inherently unsteady behavior. Rotational flow is measured in terms of the vorticity *ω*, defined as the curl of a velocity vector **v**,

$$
\rho = \nabla \times \mathbf{v} \tag{3}
$$

The related concept of circulation Γ is defined as a contour integral of vorticity by

$$
\Gamma = \oint\_{\text{\\$\\$}} \mathbf{v} \cdot \mathbf{ds} = -\iint\_{\text{\\$}} w \cdot \hat{\mathbf{n}} \, \text{d} \mathbf{S} \tag{4}
$$

The concept of a vortex is that of vorticity concentrated along a path [15].

Lid driven cavity flows are geometrically simple boundary conditions testing the convective and viscous portions of the Navier Stokes equation in an enclosed unsteady environment. The cavity flow is characterized by a quiescent flow with the driven upper lid providing energy transfer into the cavity through viscous stresses. Boundary layers along the side and lower surfaces develop as the Reynolds number increases, which tends to shift the vorticity center of rotation towards the center. A presence of the sharp corner at the downstream upper corner increasingly generates small scale flow features as the Reynolds number increases. Full cavity flows remain a strong research topic for acoustics and sensor deployment technologies.

For non-streamlined blunt bodies with a cross-flow, an adverse pressure gradient in the aft body tends to promote flow separation and an unsteady flow field. The velocity field develops into an oscillating separation line on the upper and lower surfaces. This manifests as a series of shed vortices forming and then convecting downstream with the mean flow. The von Kármán vortex street is named after the engineer and fluid

*A Fully Discrete SIPG Method for Solving Two Classes of Vortex Dominated Flows DOI: http://dx.doi.org/10.5772/intechopen.94316*

dynamicist Theodore Kármán (1963; 1994). Vortex streets are ubiquitous in nature and are visibly seen in river currents downstream of obstacles, atmospheric phenomena, and the clouds of Jupiter (e.g. The Great Red Spot). Shed vortices are also the primary driver for the zig-zag motion of bubbles in carbonated drinks. The bubble rising through the drink creates a wake of shed vorticity which impacts the integrated pressures causing side forces and thus side accelerations. The physics of sound generation with an Aeoleans harp operates by alternating vortices creating harmonic surfaces pressure variations leading to radiated acoustic tones. Tones generated by vortex shedding are the so-called Strouhal friction tones. If the diameter of the string, or cylinder immersed in the flow is *D* and the free stream velocity of the flow is **u**<sup>∞</sup> then the shedding frequency *f* of the sound is given by the Strouhal formula

*St* <sup>¼</sup> *fD* **u**<sup>∞</sup> , (5)

where *<sup>f</sup>* <sup>¼</sup> *<sup>T</sup>*�<sup>1</sup> , and *St* is the Strouhal number named after Vincent Strouhal, a Czech physicist who experimented in 1878 with wires experiencing vortex shedding and singing in the wind (Strouhal, 1878; White, 1999). The Strouhal formula provides an experimentally derived shedding frequency for fluid flow. Therefore, we are interested in an investigation of Stouhal numbers of incompressible flow at different Reynolds number.

Let Ω be a bounded polygonal domain in <sup>2</sup> . The dynamics of an incompressible fluid flow in 2D is described by the Navier–Stokes equations, which include the equations of continuity and momentum, written in dimensionless form [8] as follows:

$$\frac{\partial \mathbf{u}}{\partial t} - \nu \Delta \mathbf{u} + (\mathbf{u} \cdot \nabla) \mathbf{u} + \nabla p = \mathbf{f}, \qquad \text{in} \quad \Omega \times (0, T) \tag{6}$$

$$\nabla \cdot \mathbf{u} = \mathbf{0}, \qquad \qquad \text{in} \quad \Omega \times (\mathbf{0}, T) \tag{7}$$

$$\left.\mathbf{u}\right|\_{t=0} = \mathbf{u}\_0,\tag{8}$$

subject to the boundary conditions on ∂Ω :

$$a\mathbf{u} + (1 - a)\frac{\partial \mathbf{u}}{\partial \mathbf{n}} = \mathbf{u}\_{\infty}.$$

Here the parameter *α* has the limit values of 0 for the free-slip (no stress) condition (Neumann) and 1 for the no-slip condition (Dirichlet); **u** ¼ ð Þ *u*, *v* is the velocity; *t* is the time; and *p* is the pressure. In general, the external force **f** is not taken into account in Eq. (6).

Using the divergence free constraint, problem (6)–(8) can be rewritten in the following conservative flux form [16]:

$$\frac{\partial \mathbf{u}}{\partial t} - \nu \Delta \mathbf{u} + \nabla \cdot \mathbf{F} + \nabla p = \mathbf{f}, \qquad \qquad \qquad \text{in} \quad \Omega \times (0, T) \tag{9}$$

$$\nabla \cdot \mathbf{u} = \mathbf{0}, \qquad \qquad \text{in} \quad \Omega \times (\mathbf{0}, T) \tag{10}$$

$$\left.\mathbf{u}\right|\_{t=0} = \mathbf{u}\_0,\tag{11}$$

with the flux *F* being defined as

$$F(\mathbf{u}) = \mathbf{u} \otimes \mathbf{u} = \begin{bmatrix} u^2 & uv \\ uv & v^2 \end{bmatrix}. \tag{12}$$

original second-order or higher-order partial differential equations are written as a system of first order partial differential equations with designed suitable numerical fluxes. Interior penalty discontinuous Galerkin methods (IPDG) are known as the representative of primal schemes whereas local discontinuous Galerkin methods (LDG) belong to the class of mixed methods [9]. About IPDG, symmetric interior penalty Galerkin (SIPG) and non-symmetric interior penalty Galerkin (NIPG) methods were first introduced originally for elliptic problems by Wheeler [10] and Rivière et al. [11]. Then, we presented some numerical analysis and simulations on nonlinear parabolic problems [12]. Recently, some work based on the SIPG and NIPG methods has been successfully applied to the steady-state and transient Navier–Stokes Equations [2, 13, 14], with careful analysis being conducted on

The physics of Navier–Stokes flows are non-dimensionalized by Mach number

**<sup>a</sup>** , (1)

*<sup>μ</sup>* , (2)

*ω* ¼ ∇ � **v** (3)

*ω* � **ndS** ^ (4)

*<sup>M</sup>* <sup>¼</sup> **<sup>u</sup>**<sup>∞</sup>

*Re* <sup>¼</sup> *<sup>ρ</sup>***u**∞*<sup>D</sup>*

where *ρ* is the density of the fluid, and *μ* is the dynamic viscosity. The kinematic viscosity *ν* is the ratio of *μ* to *ρ*. At low Knudsen numbers, Navier–Stokes surface boundary conditions are effectively no-slip (i.e. zero velocity). Diffusion of momentum from freestream to surface no-slip velocities forms boundary layers decreasing in thickness as Reynolds number increases. Thus, the range of characteristic solution scale increases as the Reynolds number increases. Nonlinear convective terms coupled with the strong velocity gradients in the Navier Stokes equations drive fluid flow at even moderate Reynolds numbers to inherently unsteady behavior. Rotational flow is measured in terms of the vorticity *ω*, defined

The related concept of circulation Γ is defined as a contour integral of vorticity by

ðð

**S**

**v** � **ds** ¼ �

The concept of a vortex is that of vorticity concentrated along a path [15]. Lid driven cavity flows are geometrically simple boundary conditions testing the

unsteady environment. The cavity flow is characterized by a quiescent flow with the driven upper lid providing energy transfer into the cavity through viscous stresses. Boundary layers along the side and lower surfaces develop as the Reynolds number increases, which tends to shift the vorticity center of rotation towards the center. A presence of the sharp corner at the downstream upper corner increasingly generates small scale flow features as the Reynolds number increases. Full cavity flows remain a

For non-streamlined blunt bodies with a cross-flow, an adverse pressure gradient in the aft body tends to promote flow separation and an unsteady flow field. The velocity field develops into an oscillating separation line on the upper and lower surfaces. This manifests as a series of shed vortices forming and then convecting downstream with the mean flow. The von Kármán vortex street is named after the engineer and fluid

convective and viscous portions of the Navier Stokes equation in an enclosed

strong research topic for acoustics and sensor deployment technologies.

Γ ¼ ∮ *∂S*

optimal error estimates for the velocity.

*Vortex Dynamics Theories and Applications*

*M* and Reynolds number *Re*,

as the curl of a velocity vector **v**,

**82**

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

$$\mathcal{N}(\mathbf{u}) = \begin{pmatrix} \frac{\partial(u^2)}{\partial \mathbf{x}} + \frac{\partial(uv)}{\partial y} \\ \frac{\partial(uv)}{\partial \mathbf{x}} + \frac{\partial(v^2)}{\partial y} \end{pmatrix} = \nabla \cdot F(\mathbf{u}).$$

Galerkin methods to investigate dynamical behavior of vortex dominated lid-driven

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

The chapter is organized following [17]. In Section 2, a temporal discretization for the Navier–Stokes equation is listed with operator-splitting techniques, and subsequently, the nonlinear term is linearized. Both pressure and velocity field can be solved successively from linear elliptic and Helmholtz-type problems, respectively. In Section 3, a local numerical flux will be given for the nonlinear convection term and an SIPG scheme will be used in spacial discretization for those linear elliptic and Helmholtz-type problems with appropriate boundary conditions, and in Section 4, simulation results are presented for a lid-driven cavity flow up to *Re* ¼ 7500 and a transient flow past a circular cylinder, while numerical investigation on the Strouhal-Reynolds-number has been done, comparable to the experimental values from physics. Finally, Section 5 concludes with a brief summary.

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>*

¼ �*β*0<sup>N</sup> **<sup>u</sup>***<sup>n</sup>* ð Þ� *<sup>β</sup>*1<sup>N</sup> **<sup>u</sup>***<sup>n</sup>*�<sup>1</sup> � *<sup>β</sup>*2<sup>N</sup> **<sup>u</sup>***<sup>n</sup>*�<sup>2</sup> <sup>þ</sup> **<sup>f</sup>**ð Þ *tn*þ<sup>1</sup> ,

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

*Δt*≈ *O* L

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

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

2

, *<sup>α</sup>*<sup>2</sup> <sup>¼</sup> <sup>1</sup> 3

Especially, by using the Euler forward discretization at the first time step (*n* ¼ 0),

<sup>6</sup> , *<sup>α</sup>*<sup>0</sup> <sup>¼</sup> 3, *<sup>α</sup>*<sup>1</sup> ¼ � <sup>3</sup>

*<sup>Δ</sup><sup>t</sup>* ¼ �*β*0<sup>N</sup> **<sup>u</sup>***<sup>n</sup>* ð Þ� *<sup>β</sup>*1<sup>N</sup> **<sup>u</sup>***<sup>n</sup>*�<sup>1</sup> <sup>þ</sup> **<sup>f</sup>**ð Þ *tn*þ<sup>1</sup> (15)

, *β*<sup>0</sup> ¼ 2, *β*<sup>1</sup> ¼ �1*:*

<sup>U</sup>*N*<sup>2</sup> ,

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)–

*<sup>Δ</sup><sup>t</sup>* � *<sup>ν</sup>*Δ**u***<sup>n</sup>*þ<sup>1</sup> <sup>þ</sup> <sup>∇</sup>*pn*þ<sup>1</sup>

<sup>∇</sup> � **<sup>u</sup>***<sup>n</sup>*þ<sup>1</sup> <sup>¼</sup> 0, (14)

and *tn* ¼ *nΔt*. The semi-discrete forms of problem (6)–(8) at time *tn*þ<sup>1</sup> is

*<sup>γ</sup>*0**u***<sup>n</sup>*þ<sup>1</sup> � *<sup>α</sup>*0**u***<sup>n</sup>* � *<sup>α</sup>*1**u***<sup>n</sup>*�<sup>1</sup> � *<sup>α</sup>*2**u***<sup>n</sup>*�<sup>2</sup>

(14) can be decomposed into three stages as follows.

*<sup>γ</sup>*0**u**<sup>~</sup> � *<sup>α</sup>*0**u***<sup>n</sup>* � *<sup>α</sup>*1**u***<sup>n</sup>*�<sup>1</sup> � *<sup>α</sup>*2**u***<sup>n</sup>*�<sup>2</sup>

*<sup>γ</sup>*<sup>0</sup> <sup>¼</sup> <sup>11</sup>

we can get a medium velocity field **u**<sup>1</sup> by

• The first stage

formula can be used

**85**

*Δt* ,

(13)

and cylinder flows.

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

**2. Temporal splitting scheme**

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 triangles, where *h* is the maximum diameter of elements. Let Γ*<sup>I</sup> <sup>h</sup>* be the set of all interior edges of <sup>E</sup>*<sup>h</sup>* and <sup>Γ</sup>*<sup>B</sup> <sup>h</sup>* be the set of all boundary edges. Set <sup>Γ</sup>*<sup>h</sup>* <sup>¼</sup> <sup>Γ</sup>*<sup>I</sup> <sup>h</sup>* ∪ Γ*<sup>B</sup> <sup>h</sup> :* For any nonnegative integer *r* and *s*≥1, the classical Sobolev space on a domain *E* ⊂ <sup>2</sup> is

$$\mathcal{W}^{r,s}(E) = \{ v \in L^s(E) : \forall |m| \le r, \partial^m v \in L^s(E) \}.$$

We define the spaces of discontinuous functions

$$W = \{ \mathbf{v} \in L^2(\Omega)^2 \, : \quad \forall E \in \mathcal{E}\_h, \quad \mathbf{v} \big|\_{E} \in \left( W^{2, 4/3}(E) \right)^2 \},$$

$$M = \{ q \in L^2(\Omega) \, : \quad \forall E \in \mathcal{E}\_h, \; q \big|\_{E} \in \mathcal{W}^{1, 4/3}(E) \}.$$

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

$$\left[\phi\right] = \left(\phi\big|\_{E\_k}\right)\big|\_{\epsilon} - \left(\phi\big|\_{E\_l}\right)\big|\_{\epsilon},$$

$$\left\{\phi\right\} = \frac{1}{2}\left(\left(\phi\_{E\_k}\right)\big|\_{\epsilon} + \left(\phi\big|\_{E\_l}\right)\big|\_{\epsilon}\right).$$

Further, let **v** be a piecewise smooth vector-, or matrix-valued function at **x**∈*e* and denote its jump by

$$[\mathbf{v}] \coloneqq \mathbf{v}^+ \cdot \mathbf{n}\_{E^+} + \mathbf{v}^- \cdot \mathbf{n}\_{E^-},$$

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 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> *i*,*j TijSij*.

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* and *Mh*ð Þ ⊂ *M* for integer *N* ≥1:

$$\mathbf{V}\_{h} = \left\{ \mathbf{v} \in L^{2}(\Omega)^{2} \,:\, \forall E \in \mathcal{E}\_{h}, \quad \mathbf{v}\_{h} \in \left(\mathbb{P}\_{N}(E)\right)^{2} \right\},$$

$$M\_{h} = \left\{ q \in L^{2}(\Omega) \,:\, \forall E \in \mathcal{E}\_{h}, \,\, q \in \mathbb{P}\_{N-1}(E) \right\}.$$

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

*A Fully Discrete SIPG Method for Solving Two Classes of Vortex Dominated Flows DOI: http://dx.doi.org/10.5772/intechopen.94316*

Galerkin methods to investigate dynamical behavior of vortex dominated lid-driven and cylinder flows.

The chapter is organized following [17]. In Section 2, a temporal discretization for the Navier–Stokes equation is listed with operator-splitting techniques, and subsequently, the nonlinear term is linearized. Both pressure and velocity field can be solved successively from linear elliptic and Helmholtz-type problems, respectively. In Section 3, a local numerical flux will be given for the nonlinear convection term and an SIPG scheme will be used in spacial discretization for those linear elliptic and Helmholtz-type problems with appropriate boundary conditions, and in Section 4, simulation results are presented for a lid-driven cavity flow up to *Re* ¼ 7500 and a transient flow past a circular cylinder, while numerical investigation on the Strouhal-Reynolds-number has been done, comparable to the experimental values from physics. Finally, Section 5 concludes with a brief summary.
