**3. The spatial discretizations**

**<sup>u</sup>**<sup>1</sup> � **<sup>u</sup>**<sup>0</sup> *Δt*

, *<sup>α</sup>*<sup>0</sup> <sup>¼</sup> 2, *<sup>α</sup>*<sup>1</sup> ¼ � <sup>1</sup>

*γ*0 **<sup>u</sup>**~~ � **<sup>u</sup>**<sup>~</sup>

*<sup>n</sup>*þ<sup>1</sup> � *<sup>β</sup>*0**<sup>n</sup>** � *<sup>∂</sup>***u***<sup>n</sup>*

*∂t*

*∂t*

�*β*1**<sup>n</sup>** � *<sup>∂</sup>***u***<sup>n</sup>*�<sup>1</sup>

�*β*2**<sup>n</sup>** � *<sup>∂</sup>***u***<sup>n</sup>*�<sup>2</sup>

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

*γ*0

which can be written as a Helmholtz equation for the velocity

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

*νΔt*

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)

**<sup>u</sup>***<sup>n</sup>*þ<sup>1</sup> � **<sup>u</sup>**~~

*<sup>Δ</sup><sup>t</sup>* <sup>¼</sup> *<sup>ν</sup>*Δ**u***<sup>n</sup>*þ<sup>1</sup>

**<sup>u</sup>***<sup>n</sup>*þ<sup>1</sup> <sup>¼</sup> *<sup>γ</sup>*<sup>0</sup> *νΔt*

to update the intermediate velocity **u**~~ by (16).

• The third stage is completed by solving

**86**

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

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

schemes need to be verified with stability analysis.

The pressure projection is as follows

*∂p<sup>n</sup>*þ<sup>1</sup> *<sup>∂</sup>***<sup>n</sup>** <sup>¼</sup> **<sup>f</sup>**

scheme for the time level (*n* ¼ 2) in (15)

*Vortex Dynamics Theories and Applications*

*<sup>γ</sup>*<sup>0</sup> <sup>¼</sup> <sup>3</sup> 2

• The second stage

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

¼ �N **<sup>u</sup>**<sup>0</sup> <sup>þ</sup> **<sup>f</sup>**ð Þ *<sup>t</sup>*<sup>1</sup> ,

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

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

*:* (16)

∇ � **u**~, (17)

≔ *Gn:*

*<sup>n</sup>* <sup>¼</sup> *<sup>n</sup>* � *<sup>Δ</sup>t*. Then we use *<sup>p</sup><sup>n</sup>*þ<sup>1</sup>

**u**~~*:* (18)

which adopts the following coefficients to construct a second-order difference

2

Note that the coefficients in (15) are adjustable, but high-order time discrete

*<sup>Δ</sup><sup>t</sup>* ¼ �∇*pn*þ<sup>1</sup>

*Δt*

<sup>þ</sup> <sup>∇</sup> � *<sup>F</sup>* **un** ð Þ� *<sup>ν</sup>*Δ**u***<sup>n</sup>* 

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

<sup>þ</sup> <sup>∇</sup> � *<sup>F</sup>* **<sup>u</sup><sup>n</sup>**�**<sup>2</sup>** � *<sup>ν</sup>*Δ**u***<sup>n</sup>*�<sup>2</sup> 

,

with a Neumann boundary condition being implemented on the boundaries as

�Δ*pn*þ<sup>1</sup> ¼ � *<sup>γ</sup>*<sup>0</sup>

*∂t*

For spacial approximations, assume that piecewise polynomials of order *N* are employed, then the approximation space can be rewritten as **<sup>V</sup>***<sup>h</sup>* <sup>¼</sup> <sup>⊕</sup>*<sup>K</sup> <sup>k</sup>*¼<sup>1</sup>*N*ð Þ *Ek* 2 . In 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., *<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> -projection of *F*ð Þ **u***<sup>h</sup>* on 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, we have

$$\int\_{\Omega} (\nabla \cdot F) \cdot \mathbf{v}\_{h} d\mathbf{x} = -\sum\_{E\_{k}} \int\_{E\_{k}} (F \cdot \nabla) \cdot \mathbf{v}\_{h} d\mathbf{x} + \sum\_{\epsilon \in \Gamma\_{k}} \int\_{\epsilon} \mathbf{n}\_{\epsilon} \cdot [F \cdot \mathbf{v}\_{h}] d\mathbf{s},\tag{19}$$

where the term ð Þ� *F* � ∇ **v***<sup>h</sup>* equals to *Fij ∂vhi ∂x <sup>j</sup>* , for *i*, *j* ¼ 1, 2 and the indexes *i*, *j* 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 may define *λ*<sup>þ</sup> *<sup>E</sup>*1,*<sup>e</sup>*, *λ*� *<sup>E</sup>*2,*<sup>e</sup>* the largest eigenvalue of the Jacobians *<sup>∂</sup> <sup>∂</sup>***<sup>u</sup>** ð Þ *F* � **n***<sup>e</sup>* � � **u***E*<sup>1</sup> and *∂ <sup>∂</sup>***<sup>u</sup>** ð Þ *F* � **n***<sup>e</sup>* � � **u***E*<sup>2</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 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 integrand in the surface integral as

$$\mathbf{n}\_{\epsilon} \cdot [F \cdot \mathbf{v}\_{h}] = \mathbf{n}\_{\epsilon} \cdot \{F\} \cdot [\mathbf{v}\_{h}] + \frac{\lambda\_{\epsilon}}{2} [\mathbf{u}\_{h}] \cdot [\mathbf{v}\_{h}],$$

with *λ<sup>e</sup>* ¼ max *λ*<sup>þ</sup> *<sup>E</sup>*1,*<sup>e</sup>*, *λ*� *E*2,*e* � �, one can get a DG discretization for the nonlinear term in (19) by the local Lax-Friedrichs flux.

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 continuity on the interior edges by a penalty term, we have for (17)

$$a\left(p\_h^{n+1}, \phi\_h\right) = L\_p(\phi\_h), \ \phi \in M\_h,$$

where

$$a\left(p\_h^{n+1}, \phi\_h\right) = \sum\_{E\_k \in \mathcal{E}\_h} \int\_{E\_k} \nabla p\_h^{n+1} \cdot \nabla \phi\_h d\mathbf{x} - \sum\_{e\_k \in \Gamma\_k} \int\_{e\_k} \left\{\frac{\partial p\_h^{n+1}}{\partial \mathbf{n}}\right\} [\phi\_h] d\mathbf{s}$$

$$- \sum\_{e\_k \in \Gamma\_k} \int\_{e\_k} \left\{\frac{\partial \phi\_h}{\partial \mathbf{n}}\right\} [p\_h^{n+1}] ds + \sum\_{e\_k \in \Gamma\_k} \frac{\sigma\_\epsilon}{|e\_k|^\beta} \int\_{e\_k} \left[p\_h^{n+1}\right] [\phi\_h] ds,$$

$$\text{and} \quad L\_p(\phi\_h) = \sum\_{E\_k \in \mathcal{E}\_k} \int\_{E\_k} \frac{\chi\_0}{\Delta t} \nabla \cdot \bar{\mathbf{u}} \phi\_h d\mathbf{x} + \sum\_{e\_k \in \partial \Omega} \int\_{e\_k} \phi\_h G\_n.$$

*Example 1.* The lid-driven boundary conditions are given by:

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

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

discrete SIPG approach.

**Figure 1.**

**Figure 2.**

**89**

*An initial locally refined mesh.*

*u x*ð Þ¼ , 1 1; *u*ð Þ¼ 1, *y* 0; *u*ð Þ¼ 0, *y* 0; *u x*ð Þ¼ , 0 0; *v x*ð Þ¼ , 1 1; *v*ð Þ¼ 1, *y* 0; *v*ð Þ¼ 0, *y* 0; *v x*ð Þ¼ , 0 0*:*

Here the mesh size of the initial coarse grid is 0.2 and then it is uniformly refined

The boundary condition at the vertex is a jump from zero velocity on the edge to a unit velocity on the upper edge. Nature prevents this singularity with a boundary layer forming along all walls, making the vertex velocity zero. It is reasonable to adopt adaptive meshes for solving those singularity problems. Here, we apply the semi-implicit SIPG method with approximation polynomials of order *N* ¼ 3 in a locally refined mesh in **Figure 1** to solve the incompressible flow. In **Figure 2**, the velocity profiles of ð Þ *u*, *v* through the geometric center of the cavity are plotted with *Re* ¼ 1000, 5000, 7500 taken. From **Figures 3**–**5**, with different Reynolds numbers taken up to 7500, the vorticity field exhibits the expected characteristics of a driven cavity flow consisting of a region of vortical flow centrally located. Energy enters the cavity through the viscous boundary later formed by velocity gradients on the upper driven edge. Convection distributes flow properties throughout the domain.

three times with piecewise discontinuous elements being applied into the fully

*Velocity profiles u*ð Þ , *v through geometric center of the cavity for Re* ¼ 100, 400, 1000, 5000, 7500*.*

In general, *σ<sup>e</sup>* shall be chosen sufficiently large to guarantee coercivity, more accurately, the threshold values of *σ<sup>e</sup>* in [22] are given for *β* ¼ 1 in the above formula, which is referred to an SIPG scheme. Especially, as *β* > 1, the scheme is referred to an over-penalized scheme and the threshold values of *σ<sup>e</sup>* are presented in [23, 24]. Analogously, the SIPG discretization for (18) is given by

$$a\left(\mathbf{u}\_h^{n+1}, \mathbf{v}\_h\right) + \frac{\gamma\_0}{\nu \Delta t} \left(\mathbf{u}\_h^{n+1}, \mathbf{v}\_h\right)\_{\Omega} = L\_{\mathbf{u}}(\mathbf{v}\_h), \quad \forall \mathbf{v}\_h \in \mathbf{V}\_h,$$

where

$$a\left(\mathbf{u}\_{h}^{n+1},\mathbf{v}\_{h}\right) = \sum\_{E\_{k}\in\mathcal{E}\_{k}}\int\_{E\_{k}}\nabla\mathbf{u}\_{h}^{n+1}:\nabla\mathbf{v}\_{h}d\mathbf{x} - \sum\_{e\_{k}\in\Gamma\_{k}}\int\_{e\_{k}}\left\{\nabla\mathbf{u}\_{h}^{n+1}\right\}\mathbf{n}\cdot[\mathbf{v}\_{h}]ds\tag{20}$$

$$-\varepsilon \sum\_{e\_k \in \Gamma\_h} \int\_{e\_k} \{\nabla \mathbf{v}\_h\} \mathbf{n} \cdot \left[\mathbf{u}\_h^{n+1}\right] ds + \sum\_{e\_k \in \Gamma\_h} \frac{\sigma\_\varepsilon}{|\mathbf{c}\_k|^\beta} \int\_{e\_k} \left[\mathbf{u}\_h^{n+1}\right] \cdot [\mathbf{v}\_h] ds,\tag{21}$$

and

$$L\_{\mathbf{u}}(\mathbf{v}\_{h}) = \left(\frac{\mathcal{Y}\_{0}}{\nu \Delta t} \tilde{\mathbf{\tilde{u}}}\_{h}, \mathbf{v}\_{h}\right)\_{\Omega} - \sum\_{c\_{k} \in \partial \Omega} \int\_{c\_{k}} \left(\boldsymbol{\varepsilon} \nabla \mathbf{v}\_{h} \cdot \mathbf{n}\_{\varepsilon} - \frac{\sigma\_{\varepsilon}}{|c\_{k}|^{\beta}} \mathbf{v}\_{h}\right) \mathbf{u}\_{0} \dots$$

where *β* ¼ 1, *ε* ¼ 1*:* Note that the parameter *ε* can be 1, �1, and 0, the scheme (20) becomes SIPG, NIPG and IIPG, respectively. The SIPG scheme exhibits a stiffness matrix with symmetric structure. As a DG method, the SIPG scheme has some attractive advantages of DG methods including high order *hp*� approximation, local mass conservation, robustness and accuracy of DG methods for models with discontinuous coefficients and easy implementation on unstructured grids, while the flexibility of *p*-adaptivity (different orders of polynomials might be used for different elements) in DG methods has become competitive for modeling a wide range of engineering problems.

## **4. Numerical results**

We present a lid-driven flow problem to verify the efficiency and robustness of the interior penalty discontinuous Galerkin method, and then investigate a flow past a cylinder with walls or without a wall, as well as the relationship between the Strouhal number and the Reynold number. Throughout the section, time steps *Δt*≤ 1*E* � 03 are taken.

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

*Example 1.* The lid-driven boundary conditions are given by:

$$\begin{cases} u(\mathbf{x}, \mathbf{1}) = \mathbf{1}; \ u(\mathbf{1}, \mathbf{y}) = \mathbf{0}; \ u(\mathbf{0}, \mathbf{y}) = \mathbf{0}; \ u(\mathbf{x}, \mathbf{0}) = \mathbf{0};\\ v(\mathbf{x}, \mathbf{1}) = \mathbf{1}; \ v(\mathbf{1}, \mathbf{y}) = \mathbf{0}; \ v(\mathbf{0}, \mathbf{y}) = \mathbf{0}; \ v(\mathbf{x}, \mathbf{0}) = \mathbf{0}. \end{cases}$$

Here the mesh size of the initial coarse grid is 0.2 and then it is uniformly refined three times with piecewise discontinuous elements being applied into the fully discrete SIPG approach.

The boundary condition at the vertex is a jump from zero velocity on the edge to a unit velocity on the upper edge. Nature prevents this singularity with a boundary layer forming along all walls, making the vertex velocity zero. It is reasonable to adopt adaptive meshes for solving those singularity problems. Here, we apply the semi-implicit SIPG method with approximation polynomials of order *N* ¼ 3 in a locally refined mesh in **Figure 1** to solve the incompressible flow. In **Figure 2**, the velocity profiles of ð Þ *u*, *v* through the geometric center of the cavity are plotted with *Re* ¼ 1000, 5000, 7500 taken. From **Figures 3**–**5**, with different Reynolds numbers taken up to 7500, the vorticity field exhibits the expected characteristics of a driven cavity flow consisting of a region of vortical flow centrally located. Energy enters the cavity through the viscous boundary later formed by velocity gradients on the upper driven edge. Convection distributes flow properties throughout the domain.

**Figure 1.** *An initial locally refined mesh.*

where

where

and

*a* **u***<sup>n</sup>*þ<sup>1</sup> *<sup>h</sup>* , **v***<sup>h</sup>* � � <sup>¼</sup> <sup>X</sup>

> �*ε* X *ek* ∈Γ*<sup>h</sup>*

*L***u**ð Þ¼ **v***<sup>h</sup>*

range of engineering problems.

**4. Numerical results**

*Δt*≤ 1*E* � 03 are taken.

**88**

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

*Vortex Dynamics Theories and Applications*

*Ek* ∈E*<sup>h</sup>*

ð *ek*

and *Lp <sup>ϕ</sup><sup>h</sup>* ð Þ¼ <sup>X</sup>

� <sup>X</sup> *ek* ∈Γ*<sup>h</sup>*

*a* **u***<sup>n</sup>*þ<sup>1</sup> *<sup>h</sup>* , **v***<sup>h</sup>* � � <sup>þ</sup> *<sup>γ</sup>*<sup>0</sup>

*Ek* ∈E*<sup>h</sup>*

ð *ek* ð *Ek*

*γ*0 *νΔt* ð *Ek*

*∂ϕ<sup>h</sup> ∂***n** � �

[23, 24]. Analogously, the SIPG discretization for (18) is given by

**u***<sup>n</sup>*þ<sup>1</sup> *<sup>h</sup>* , **v***<sup>h</sup>* � �

*<sup>h</sup>* : <sup>∇</sup>**v***hd***<sup>x</sup>** � <sup>X</sup>

� �*ds* <sup>þ</sup> <sup>X</sup>

� <sup>X</sup> *ek* ∈ ∂Ω *ek* ∈Γ*<sup>h</sup>*

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

ð *ek*

where *β* ¼ 1, *ε* ¼ 1*:* Note that the parameter *ε* can be 1, �1, and 0, the scheme (20) becomes SIPG, NIPG and IIPG, respectively. The SIPG scheme exhibits a stiffness matrix with symmetric structure. As a DG method, the SIPG scheme has some attractive advantages of DG methods including high order *hp*� approximation, local mass conservation, robustness and accuracy of DG methods for models with discontinuous coefficients and easy implementation on unstructured grids, while the flexibility of *p*-adaptivity (different orders of polynomials might be used for different elements) in DG methods has become competitive for modeling a wide

We present a lid-driven flow problem to verify the efficiency and robustness of the interior penalty discontinuous Galerkin method, and then investigate a flow past a cylinder with walls or without a wall, as well as the relationship between the Strouhal number and the Reynold number. Throughout the section, time steps

ð *ek*

*σe* j j *ek <sup>β</sup>* ð *ek*

*<sup>ε</sup>*∇**v***<sup>h</sup>* � **<sup>n</sup>***<sup>e</sup>* � *<sup>σ</sup><sup>e</sup>*

!

*νΔt*

∇**u***<sup>n</sup>*þ<sup>1</sup>

f g <sup>∇</sup>**v***<sup>h</sup>* **<sup>n</sup>** � **<sup>u</sup>***<sup>n</sup>*þ<sup>1</sup>

**u**~~*h*, **v***<sup>h</sup>* � � *h*

Ω

∇*pn*þ<sup>1</sup>

*pn*þ<sup>1</sup> *h*

> ð *Ek γ*0 *Δt*

In general, *σ<sup>e</sup>* shall be chosen sufficiently large to guarantee coercivity, more accurately, the threshold values of *σ<sup>e</sup>* in [22] are given for *β* ¼ 1 in the above formula, which is referred to an SIPG scheme. Especially, as *β* > 1, the scheme is referred to an over-penalized scheme and the threshold values of *σ<sup>e</sup>* are presented in

*Ek* ∈E*<sup>h</sup>*

*<sup>h</sup>* � <sup>∇</sup>*ϕhdx* � <sup>X</sup>

� �*ds* <sup>þ</sup> <sup>X</sup>

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

*σe* j j *ek <sup>β</sup>*

<sup>∇</sup> � **<sup>u</sup>**~*ϕhdx* <sup>þ</sup> <sup>X</sup>

<sup>Ω</sup> ¼ *L***u**ð Þ **v***<sup>h</sup>* , ∀**v***<sup>h</sup>* ∈ **V***h*,

∇**u***<sup>n</sup>*þ<sup>1</sup> *h*

> **u***<sup>n</sup>*þ<sup>1</sup> *h*

> > j j *ek <sup>β</sup>* **<sup>v</sup>***<sup>h</sup>*

� �**<sup>n</sup>** � ½ � **<sup>v</sup>***<sup>h</sup> ds* (20)

� � � ½ � **<sup>v</sup>***<sup>h</sup> ds*, (21)

**u**0*:*

ð *ek pn*þ<sup>1</sup> *h* � � *<sup>ϕ</sup><sup>h</sup>* ½ �*ds*,

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

ð *ek*

*∂pn*þ<sup>1</sup> *h ∂***n** � �

*ek* ∈ ∂Ω

ð *ek ϕhGn:*

*ϕ<sup>h</sup>* ½ �*ds*

**Figure 2.** *Velocity profiles u*ð Þ , *v through geometric center of the cavity for Re* ¼ 100, 400, 1000, 5000, 7500*.*

**Figure 3.** *Re =1000, N = 4, mesh* #2*. Left: pressure contour; right: vorticity contour.*

Cylinder flow contains the fundamentals of unsteady fluid dynamics in a simplified geometry. The flow properties and unsteadiness are well defined through years of experimental measurements across a wide range of Reynolds numbers [25], making the cylinder an ideal validation testcase for unsteady numerical fluid

Verification in a numerical domain requires insights from physics for a proper comparison to experimental and theoretical data. In **Figure 7**, we observe that the boundary layers forms along the upper and lower walls. From continuity of mass, the presence of a boundary layer decreasing the flow velocity near the wall requires an increase in the centerline average flow velocity. The cylinders wake provides a similar increase in centerline velocities. This implies a non-intuitive reality that drag can increase velocities within constrained domains. This effect is compensated for in wind tunnel test [26] environments topologically similar to **Figure 7** with a constant mass flow rate and no-slip walls. Drela [15] develops an analysis for 2D

wind tunnels resulting in an effective coefficient of drag of

*Cylinder flow in a channel N* ¼ 3*. Top: with walls; bottom: without a wall.*

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

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

and an effective Reynolds number of

*Cd* <sup>¼</sup> <sup>1</sup> � <sup>1</sup>

*Re* ¼ 1 þ

*Cylinder flow with walls, Re* ¼ 100*. Top: vorticity contour; bottom: pressure contour.*

2 *c <sup>H</sup> Cd* � *<sup>π</sup>* 2 *A H*2

1 4 *c <sup>H</sup> Cd* <sup>þ</sup> *<sup>π</sup>* 6 *A H*2

where the *un* subscript represents the uncorrected value, *H* represents the domain height, *A* represents the cylinder area and *c* represents the cylinder radius. Drelas analysis does not specifically include the boundary layer forming on the

*Cdun* ,

*Re un*,

dynamics simulations.

**Figure 6.**

**Figure 7.**

**91**

**Figure 4.**

*Re =5000, N = 3, the initial mesh refined. Left: pressure contour; right: vorticity contour.*

#### **Figure 5.**

*Re =7500, N = 3, the initial mesh refined. Left: pressure contour; right: vorticity contour.*

Moreover, a video on the dynamical evolution of vorticity isolines ( *Re* ¼ 1000, *N* ¼ 4) can be browsed through a website (Available from: https://youtu.be/UfGWvnoiW58). These numerical simulations are performed for the Navier–Stokes equations which illustrate the effectiveness of the DG method.

*Example 2.* We simulate a channel flow past a circular cylinder with a radius 0*:*05 at the origin 0, 0 ð Þ for *Re* ¼ 100 by the discontinuous Galerkin method in the domain ð Þ� � �1, 3 ð Þ 0*:*5, 0*:*5 . The free stream velocity on the inflow boundary is **<sup>u</sup>**<sup>∞</sup> <sup>¼</sup> ð Þ 1, 0 , while the outflow boundary is *<sup>∂</sup>***<sup>u</sup>** *<sup>∂</sup>***<sup>n</sup>** ¼ 0. To the boundary conditions on the upper and lower sides, we present two different conditions for comparison (see **Figure 6**), which are wall (*u* ¼ 0, *v* ¼ 0) and homogeneous Neumann boundary conditions (*<sup>u</sup>* <sup>¼</sup> 1, *<sup>∂</sup><sup>v</sup> <sup>∂</sup><sup>n</sup>* ¼ 0), respectively. The homogeneous Neumann boundary condition is a special non-reflecting case, where the boundary flux is zero. For reference, the density of the fluid is given by *<sup>ρ</sup>* <sup>¼</sup> <sup>1</sup> *kg* � *<sup>m</sup>*�<sup>3</sup> and a locally refined mesh ( max *h* ¼ 0*:*088) will be used for the simulations.

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

**Figure 6.** *Cylinder flow in a channel N* ¼ 3*. Top: with walls; bottom: without a wall.*

Cylinder flow contains the fundamentals of unsteady fluid dynamics in a simplified geometry. The flow properties and unsteadiness are well defined through years of experimental measurements across a wide range of Reynolds numbers [25], making the cylinder an ideal validation testcase for unsteady numerical fluid dynamics simulations.

Verification in a numerical domain requires insights from physics for a proper comparison to experimental and theoretical data. In **Figure 7**, we observe that the boundary layers forms along the upper and lower walls. From continuity of mass, the presence of a boundary layer decreasing the flow velocity near the wall requires an increase in the centerline average flow velocity. The cylinders wake provides a similar increase in centerline velocities. This implies a non-intuitive reality that drag can increase velocities within constrained domains. This effect is compensated for in wind tunnel test [26] environments topologically similar to **Figure 7** with a constant mass flow rate and no-slip walls. Drela [15] develops an analysis for 2D wind tunnels resulting in an effective coefficient of drag of

$$\mathbf{C}\_d = \left(\mathbf{1} - \frac{\mathbf{1}}{2}\frac{c}{H}\mathbf{C}\_d - \frac{\pi}{2}\frac{A}{H^2}\right)\mathbf{C}\_{d\_{\text{inv}}},$$

and an effective Reynolds number of

$$Re = \left(1 + \frac{1}{4}\frac{c}{H}C\_d + \frac{\pi}{6}\frac{A}{H^2}\right)Re\_{nm},$$

where the *un* subscript represents the uncorrected value, *H* represents the domain height, *A* represents the cylinder area and *c* represents the cylinder radius. Drelas analysis does not specifically include the boundary layer forming on the

**Figure 7.** *Cylinder flow with walls, Re* ¼ 100*. Top: vorticity contour; bottom: pressure contour.*

Moreover, a video on the dynamical evolution of vorticity isolines ( *Re* ¼ 1000, *N* ¼ 4) can be browsed through a website (Available from:

*Re =7500, N = 3, the initial mesh refined. Left: pressure contour; right: vorticity contour.*

*Re =5000, N = 3, the initial mesh refined. Left: pressure contour; right: vorticity contour.*

*Re =1000, N = 4, mesh* #2*. Left: pressure contour; right: vorticity contour.*

*Vortex Dynamics Theories and Applications*

refined mesh ( max *h* ¼ 0*:*088) will be used for the simulations.

**<sup>u</sup>**<sup>∞</sup> <sup>¼</sup> ð Þ 1, 0 , while the outflow boundary is *<sup>∂</sup>***<sup>u</sup>**

boundary conditions (*<sup>u</sup>* <sup>¼</sup> 1, *<sup>∂</sup><sup>v</sup>*

**Figure 3.**

**Figure 4.**

**Figure 5.**

**90**

https://youtu.be/UfGWvnoiW58). These numerical simulations are performed for the Navier–Stokes equations which illustrate the effectiveness of the DG method. *Example 2.* We simulate a channel flow past a circular cylinder with a radius 0*:*05 at the origin 0, 0 ð Þ for *Re* ¼ 100 by the discontinuous Galerkin method in the domain ð Þ� � �1, 3 ð Þ 0*:*5, 0*:*5 . The free stream velocity on the inflow boundary is

the upper and lower sides, we present two different conditions for comparison (see **Figure 6**), which are wall (*u* ¼ 0, *v* ¼ 0) and homogeneous Neumann

boundary condition is a special non-reflecting case, where the boundary flux is zero. For reference, the density of the fluid is given by *<sup>ρ</sup>* <sup>¼</sup> <sup>1</sup> *kg* � *<sup>m</sup>*�<sup>3</sup> and a locally

*<sup>∂</sup>***<sup>n</sup>** ¼ 0. To the boundary conditions on

*<sup>∂</sup><sup>n</sup>* ¼ 0), respectively. The homogeneous Neumann

upper and lower walls. The flow physics associated with wall boundary layer drag differs from cylinder drag in that the wall drag is a distributed effect of monotonically increasing drag with downstream distance rather than a conceptual point source of drag. The wall boundary layer tends to provide a steady acceleration of flow within the interior flow domain (i.e. non-boundary layer portion) leading to an effective buoyancy drag. A secondary feature of the wall boundary layer is that downstream flow features such as vortices are convected at a higher perturbation velocity compared to the initial upstream velocity. For numerical validation of raw experimental data, either the wind tunnel geometry should exactly match the numerical geometry, or the numerical geometry should be corrected using the concepts introduced above to match the actual wind tunnel geometry. Alternatively, the open-air corrected values should be used for validation. The above analysis provides insight into the domain height necessary to reduce volume blockage (*c=H*) and wake buoyancy (*A=H*<sup>2</sup> ) effects.

**Figure 10.**

**Figure 11.**

**Figure 12.**

**93**

*Re* ¼ 100*. Above: velocity in y-direction with walls; below: velocity in y-direction without a wall.*

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

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

*Vorticity contours of flow past a cylinder without a wall at different time.*

*Pressure contours of flow past a cylinder without a wall at different time.*

Alternatively, the flow without walls in **Figure 8** has no interference of the boundary layers along the channel on the up and bottom boundaries, thus the pressure contours expend after flow passing through the cylinder. We also compared the components of the velocity profiles along the *x*� direction in **Figure 9**, and observed that the boundary layers are produced in the top picture rather than in the bottom one. If the effect of the boundary layers disappeared, the velocity in the *x*� direction would reduce dispersively, in other words, the vortex lifespan is less than those produced in the channel with walls. The velocity profiles in the *y*� direction have been given in **Figure 10**.

We localize the domain around the cylinder and refine the mesh, then show the vorticity startup behavior in **Figure 11** as well as the pressure **Figure 12**. Upon startup, two vortices of opposite direction are formed on the upper and lower

#### **Figure 8.**

*Flow past a cylinder in a channel without a wall, Re* ¼ 100*. Above: vorticity contour; below: pressure contour.*

**Figure 9.**

*Re* ¼ 100*. Above: velocity in x-direction with walls; below: velocity in x-direction without a wall.*

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

**Figure 10.**

upper and lower walls. The flow physics associated with wall boundary layer drag differs from cylinder drag in that the wall drag is a distributed effect of monotonically increasing drag with downstream distance rather than a conceptual point source of drag. The wall boundary layer tends to provide a steady acceleration of flow within the interior flow domain (i.e. non-boundary layer portion) leading to an effective buoyancy drag. A secondary feature of the wall boundary layer is that downstream flow features such as vortices are convected at a higher perturbation velocity compared to the initial upstream velocity. For numerical validation of raw experimental data, either the wind tunnel geometry should exactly match the numerical geometry, or the numerical geometry should be corrected using the concepts introduced above to match the actual wind tunnel geometry. Alternatively, the open-air corrected values should be used for validation. The above analysis provides insight into the domain height necessary to reduce volume

) effects.

Alternatively, the flow without walls in **Figure 8** has no interference of the boundary layers along the channel on the up and bottom boundaries, thus the pressure contours expend after flow passing through the cylinder. We also compared the components of the velocity profiles along the *x*� direction in **Figure 9**, and observed that the boundary layers are produced in the top picture rather than in the bottom one. If the effect of the boundary layers disappeared, the velocity in the *x*� direction would reduce dispersively, in other words, the vortex lifespan is less than those produced in the channel with walls. The velocity profiles in the *y*�

We localize the domain around the cylinder and refine the mesh, then show the

Upon startup, two vortices of opposite direction are formed on the upper and lower

*Flow past a cylinder in a channel without a wall, Re* ¼ 100*. Above: vorticity contour; below: pressure contour.*

*Re* ¼ 100*. Above: velocity in x-direction with walls; below: velocity in x-direction without a wall.*

vorticity startup behavior in **Figure 11** as well as the pressure **Figure 12**.

blockage (*c=H*) and wake buoyancy (*A=H*<sup>2</sup>

*Vortex Dynamics Theories and Applications*

direction have been given in **Figure 10**.

**Figure 8.**

**Figure 9.**

**92**

*Re* ¼ 100*. Above: velocity in y-direction with walls; below: velocity in y-direction without a wall.*

**Figure 11.**

*Vorticity contours of flow past a cylinder without a wall at different time.*

**Figure 12.** *Pressure contours of flow past a cylinder without a wall at different time.*

**Figure 13.**

*Drag, lift coefficients and the variation of pressure with time, Re* ¼ 80 *(left), Re* ¼ 140 *(right).*

aft portion of the cylinder. Given a low total simulation time, the flow field resembles the symmetrical low Reynolds number steady flow. As time progresses however, instabilities are magnified and the upper-lower symmetry increases. Given a total time of beyond *t* ¼ 10, an autonomous and phased locked set of street vortices are generated. Surface pressures (**Figure 13** at *Re* ¼ 80, 140) generated reflect the process, including a steady startup portion and the eventual vortex shedding frequency. Validation at *Re* >41 requires sufficient time to obtain the unsteady behavior.

To reduce the effect of the boundary layers along the walls, the coefficients of drag and lift as well as the difference of pressure between the leading edge and the trailing edge on the cylinder shall be computed in a larger domain. Then in a domain Ω ≔ ½ �� � �1, 5 ½ � 1, 1 , higher order DG finite elements have been investigated. Based on the velocity **u**<sup>∞</sup> and the diameter of the cylinder *D* ¼ 0*:*1, we will chose different viscosity coefficients *ν* ¼ 1*e* � 3, 5*e* � 4, 2*:*5*e* � 4 etc. to simulate flow with different Reynolds numbers, that is, the cases *Re* ¼ 100, 200, 400, respectively. Our interest is the drag coefficient *Cd*, the lift coefficient *Cl* on the cylinder and the difference of the pressure between the front and the back of the cylinder

which comes from the classical definition (5). From the evolution of *Cd*, *Cl* and *dp* as in **Figure 13**, we may find a period *Tlift* of the lift coefficients for different values of *Re* to calculate Strouhal number by (22). In **Figure 14**, a comparison of Strouhal numbers between the experimental estimates in Fey etc. [5] and our estimates from (22) indicates a behavioral match between the unsteady onset at approximately *Re* ¼ 50 and the beginning of the transition to turbulence at *Re* ¼ 180. Beyond *Re* ¼ 180, the onset of turbulence changes the flow physics by drastically increasing the energy spectrum of the shed vorticity. The transition appears as a marked decrease in the Strouhal number prior to *Re* ¼ 200. As the present SIPG solver does not include a 3D turbulence model, our estimates follow the laminar results into the actual turbulent region. Moreover, **Figure 14** verifies availability of the periodic of vortex street replaced by the periodic of the lift coefficients. There are many results for flow past a cylinder using Reynolds Averaged Navier–Stokes (RANS) method. For large *Re* up to 2000, the readers are referred to the reference [28] for details, for example, Strouhal number of a

*A comparison of Strouhal-Reynolds-number between our estimate and the linear fit in [27] for* 50 ≤ *Re* ≤400*.*

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

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

A SIPG solver is developed for the incompressible Navier Stokes equations of fluid flow. Two testcases are presented: a lid-driven cavity and a cylinder flow. The DG method produces stable discretizations of the convective operator for high order discretizations on unstructured meshes of simplices, as a requirement for real-world complex geometries. There are still some open problems, such as how the strouhal number of the von Kármán vortex street changes against the Reynolds

The author would like to thank INTECHOPEN for their fully sponsor the publication of this chapter and waive the Open Access Publication Fee completely.

number after oscillations or noises are added to the incident flow.

cylinder flow could experience a slight fall.

**5. Conclusions**

**Figure 14.**

**Acknowledgements**

**95**

$$d\_p = p(t; -0.05, 0) - p(t; 0.05, 0).$$

We use the definition of *Cd* and *Cl* given in [4] as follows:

$$\begin{aligned} \mathbf{C}\_{d} &= \frac{2}{\rho D u\_{\text{max}}^{2}} \int\_{\mathcal{S}} \left( \rho \nu \frac{\partial \mathbf{u}\_{\text{ts}}(t)}{\partial \mathbf{n}} n\_{\text{\textgreater}} - p(t) n\_{\text{\textgreater}} \right) d\mathbf{s}, \\ \text{and} \quad \mathbf{C}\_{l} &= -\frac{2}{\rho D u\_{\text{max}}^{2}} \int\_{\mathcal{S}} \left( \rho \nu \frac{\partial \mathbf{u}\_{\text{ts}}(t)}{\partial \mathbf{n}} n\_{\text{\textgreater}} + p(t) n\_{\text{\textgreater}} \right) d\mathbf{s}, \end{aligned}$$

where **n** ¼ *nx*, *ny* � �*<sup>T</sup>* is the normal vector on the cylinder boundary *S* directing into Ω, **tS** ¼ *ny*, �*nx* � �*<sup>T</sup>* the tangential vector and **utS** the tangential velocity along *<sup>S</sup>*. In the literature, Fey et al. in [27] propose the Strouhal number represented by piecewise linear relationships of the form

$$\text{St} = \text{ST}^\* \, ^\cdot + \frac{m}{\sqrt{Re}}$$

with different values *St* <sup>∗</sup> and *m* in different shedding regimes of the 3D circular cylinder wake. We originally define the periodic *Tlift* of the lift coefficients (see **Figure 13**) by the periodic *T* ≔ <sup>1</sup> *<sup>f</sup>* appearing in the definition of Strouhal number, i.e.,

$$\text{St}(Re) = \frac{D}{T\_{lift}\mathbf{u}\_{\infty}},\tag{22}$$

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

#### **Figure 14.**

aft portion of the cylinder. Given a low total simulation time, the flow field resembles the symmetrical low Reynolds number steady flow. As time progresses however, instabilities are magnified and the upper-lower symmetry increases. Given a total time of beyond *t* ¼ 10, an autonomous and phased locked set of street vortices are generated. Surface pressures (**Figure 13** at *Re* ¼ 80, 140) generated reflect the process, including a steady startup portion and the eventual vortex shedding frequency. Validation at *Re* >41 requires sufficient time to obtain the unsteady behavior.

*Drag, lift coefficients and the variation of pressure with time, Re* ¼ 80 *(left), Re* ¼ 140 *(right).*

To reduce the effect of the boundary layers along the walls, the coefficients of drag and lift as well as the difference of pressure between the leading edge and the trailing edge on the cylinder shall be computed in a larger domain. Then in a domain Ω ≔ ½ �� � �1, 5 ½ � 1, 1 , higher order DG finite elements have been investigated. Based on the velocity **u**<sup>∞</sup> and the diameter of the cylinder *D* ¼ 0*:*1, we will chose different viscosity coefficients *ν* ¼ 1*e* � 3, 5*e* � 4, 2*:*5*e* � 4 etc. to simulate flow with different Reynolds numbers, that is, the cases *Re* ¼ 100, 200, 400, respectively. Our interest is the drag coefficient *Cd*, the lift coefficient *Cl* on the cylinder and the difference of

*dp* ¼ *p t*ð Þ� ; �0*:*05, 0 *p t*ð Þ ; 0*:*05, 0 *:*

*<sup>∂</sup>***utS** ð Þ*<sup>t</sup> ∂***n**

� �

*<sup>∂</sup>***utS** ð Þ*<sup>t</sup> ∂***n**

� �*<sup>T</sup>* is the normal vector on the cylinder boundary *S* directing

� �*<sup>T</sup>* the tangential vector and **utS** the tangential velocity along *<sup>S</sup>*.

*m* ffiffiffiffiffiffiffi *Re* <sup>p</sup>

*ny* � *p t*ð Þ*nx*

� �

*nx* þ *p t*ð Þ*ny*

*<sup>f</sup>* appearing in the definition of Strouhal number, i.e.,

*ds*,

*ds*,

, (22)

the pressure between the front and the back of the cylinder

*Cd* <sup>¼</sup> <sup>2</sup> *ρDu*<sup>2</sup> max

and *Cl* ¼ � <sup>2</sup>

piecewise linear relationships of the form

**Figure 13**) by the periodic *T* ≔ <sup>1</sup>

where **n** ¼ *nx*, *ny*

into Ω, **tS** ¼ *ny*, �*nx*

**94**

**Figure 13.**

*Vortex Dynamics Theories and Applications*

We use the definition of *Cd* and *Cl* given in [4] as follows:

*ρDu*<sup>2</sup> max

ð S *ρν*

> ð S *ρν*

In the literature, Fey et al. in [27] propose the Strouhal number represented by

*St* <sup>¼</sup> *ST*<sup>∗</sup> <sup>þ</sup>

cylinder wake. We originally define the periodic *Tlift* of the lift coefficients (see

*St Re* ð Þ¼ *<sup>D</sup>*

with different values *St* <sup>∗</sup> and *m* in different shedding regimes of the 3D circular

*Tlift***u**<sup>∞</sup>

*A comparison of Strouhal-Reynolds-number between our estimate and the linear fit in [27] for* 50 ≤ *Re* ≤400*.*

which comes from the classical definition (5). From the evolution of *Cd*, *Cl* and *dp* as in **Figure 13**, we may find a period *Tlift* of the lift coefficients for different values of *Re* to calculate Strouhal number by (22). In **Figure 14**, a comparison of Strouhal numbers between the experimental estimates in Fey etc. [5] and our estimates from (22) indicates a behavioral match between the unsteady onset at approximately *Re* ¼ 50 and the beginning of the transition to turbulence at *Re* ¼ 180. Beyond *Re* ¼ 180, the onset of turbulence changes the flow physics by drastically increasing the energy spectrum of the shed vorticity. The transition appears as a marked decrease in the Strouhal number prior to *Re* ¼ 200. As the present SIPG solver does not include a 3D turbulence model, our estimates follow the laminar results into the actual turbulent region. Moreover, **Figure 14** verifies availability of the periodic of vortex street replaced by the periodic of the lift coefficients. There are many results for flow past a cylinder using Reynolds Averaged Navier–Stokes (RANS) method. For large *Re* up to 2000, the readers are referred to the reference [28] for details, for example, Strouhal number of a cylinder flow could experience a slight fall.

## **5. Conclusions**

A SIPG solver is developed for the incompressible Navier Stokes equations of fluid flow. Two testcases are presented: a lid-driven cavity and a cylinder flow. The DG method produces stable discretizations of the convective operator for high order discretizations on unstructured meshes of simplices, as a requirement for real-world complex geometries. There are still some open problems, such as how the strouhal number of the von Kármán vortex street changes against the Reynolds number after oscillations or noises are added to the incident flow.

## **Acknowledgements**

The author would like to thank INTECHOPEN for their fully sponsor the publication of this chapter and waive the Open Access Publication Fee completely.
