Research article

A High-Order Symmetric Interior Penalty Discontinuous Galerkin Scheme to Simulate Vortex Dominated Incompressible Fluid Flow

  • Received: 16 April 2016 Accepted: 25 April 2016 Published: 27 April 2016
  • A high-order Symmetric Interior Penalty discontinuous Galerkin (SIPG) method has been used for solving the incompressible Navier-Stokes equation. We apply the temporal splitting scheme in time and the SIPG discretization in space with the local Lax-Friedrichs flux for the discretization of nonlinear terms. A fully discrete semi-implicit splitting scheme has been presented and high-order discontinuous Galerkin (DG) finite elements are available. Under a constraint of the CFL condition, two benchmark problems in 2D are investigated: one is a lid-driven cavity flow to verify the high-order discontinuous Galerkin method is accurate and robust; the other is a flow past a circular cylinder, for which we mainly check the Strouhal numbers with the von Kármán vortex street, and also simulate the boundary layers with walls and corresponding dynamical behavior with Neumann conditions on the top and bottom boundaries, respectively. We predict the Strouhal number for the range of Reynolds number 50 ≤ Re ≤ 400, making a comparison between the predicted values by our numerical method and the referenced values from physical experiments.

    Citation: Song Lunji. A High-Order Symmetric Interior Penalty Discontinuous Galerkin Schemeto Simulate Vortex Dominated Incompressible Fluid Flow[J]. AIMS Mathematics, 2016, 1(1): 43-63. doi: 10.3934/Math.2016.1.43

    Related Papers:

    [1] Cagnur Corekli . The SIPG method of Dirichlet boundary optimal control problems with weakly imposed boundary conditions. AIMS Mathematics, 2022, 7(4): 6711-6742. doi: 10.3934/math.2022375
    [2] Martin Stoll, Hamdullah Yücel . Symmetric interior penalty Galerkin method for fractional-in-space phase-field equations. AIMS Mathematics, 2018, 3(1): 66-95. doi: 10.3934/Math.2018.1.66
    [3] M. Chandru, T. Prabha, V. Shanthi, H. Ramos . An almost second order uniformly convergent method for a two-parameter singularly perturbed problem with a discontinuous convection coefficient and source term. AIMS Mathematics, 2024, 9(9): 24998-25027. doi: 10.3934/math.20241219
    [4] Lingling Sun, Hai Bi, Yidu Yang . A posteriori error estimates of mixed discontinuous Galerkin method for a class of Stokes eigenvalue problems. AIMS Mathematics, 2023, 8(9): 21270-21297. doi: 10.3934/math.20231084
    [5] Yanhua Gu . High-order numerical method for the fractional Korteweg-de Vries equation using the discontinuous Galerkin method. AIMS Mathematics, 2025, 10(1): 1367-1383. doi: 10.3934/math.2025063
    [6] Xiaoxia Wang, Jinping Jiang . The pullback attractor for the 2D g-Navier-Stokes equation with nonlinear damping and time delay. AIMS Mathematics, 2023, 8(11): 26650-26664. doi: 10.3934/math.20231363
    [7] Kaifang Liu, Lunji Song . A family of interior-penalized weak Galerkin methods for second-order elliptic equations. AIMS Mathematics, 2021, 6(1): 500-517. doi: 10.3934/math.2021030
    [8] Ruby, Vembu Shanthi, Higinio Ramos . A numerical approach to approximate the solution of a quasilinear singularly perturbed parabolic convection diffusion problem having a non-smooth source term. AIMS Mathematics, 2025, 10(3): 6827-6852. doi: 10.3934/math.2025313
    [9] Şuayip Toprakseven, Seza Dinibutun . Error estimations of a weak Galerkin finite element method for a linear system of $ \ell\geq 2 $ coupled singularly perturbed reaction-diffusion equations in the energy and balanced norms. AIMS Mathematics, 2023, 8(7): 15427-15465. doi: 10.3934/math.2023788
    [10] Liangkun Xu, Hai Bi . A multigrid discretization scheme of discontinuous Galerkin method for the Steklov-Lamé eigenproblem. AIMS Mathematics, 2023, 8(6): 14207-14231. doi: 10.3934/math.2023727
  • A high-order Symmetric Interior Penalty discontinuous Galerkin (SIPG) method has been used for solving the incompressible Navier-Stokes equation. We apply the temporal splitting scheme in time and the SIPG discretization in space with the local Lax-Friedrichs flux for the discretization of nonlinear terms. A fully discrete semi-implicit splitting scheme has been presented and high-order discontinuous Galerkin (DG) finite elements are available. Under a constraint of the CFL condition, two benchmark problems in 2D are investigated: one is a lid-driven cavity flow to verify the high-order discontinuous Galerkin method is accurate and robust; the other is a flow past a circular cylinder, for which we mainly check the Strouhal numbers with the von Kármán vortex street, and also simulate the boundary layers with walls and corresponding dynamical behavior with Neumann conditions on the top and bottom boundaries, respectively. We predict the Strouhal number for the range of Reynolds number 50 ≤ Re ≤ 400, making a comparison between the predicted values by our numerical method and the referenced values from physical experiments.


    1. Introduction

    As a fundamental equation of fluid dynamics, the Navier-Stokes equations have been investigated by many scientists conducting research on numerical schemes for their numerical solutions [7,8,10,17,18,19,22,23]. Analytical solutions of real flow problems including complex geometries are not available, nor likely in the foreseeable future. There are 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. For example, Symmetric Interior Penalty Galerkin (SIPG) and Non-symmetric Interior Penalty Galerkin (NIPG) methods were first introduced originally for elliptic problems by Wheeler [24] and Riviˊere et al. [16]. Recently, some work based on the SIPG and NIPG methods has been successfully applied to the steady-state and transient Navier-Stokes equations [8,9,15], for which optimal error estimates for the velocity have been derived.

    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. The physics of Navier-Stokes flows are non-dimensionalized by Mach number M and Reynolds number Re,

    M=ua, (1)
    Re=ρuDμ, (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 couped 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,
    w=×v (3)
    The related concept of circulation Γ is defined as a contour integral of vorticity by
    Γ=Svds=SωˆndS (4)
    The concept of a vortex is that of vorticity concentrated along a path [3].

    Lid driven cavity flows are geometrically simple boundary conditions testing the convective and viscous portions of the Navier Stokes equation in a 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ˊarmˊan vortex street is named after the engineer and fluid dynamicist Theodore Kˊarmˊan (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 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 Aeolean's 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 then the shedding frequency f of the sound is given by the Strouhal formula

    St=fDu, (5)
    where f=T1, 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 a 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 R2. 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 [23] as follows:

    utνΔu+(u)u+p=f,inΩ×(0,T) (6)
    u=0,inΩ×(0,T) (7)
    u|t=0=u0, (8)
    subject to the boundary conditions on Ω:
    αu+(1α)un=u.
    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:

    utνΔu+F+p=f,inΩ×(0,T) (9)
    u=0,inΩ×(0,T) (10)
    u|t=0=u0, (11)
    with the flux F being defined as
    F(u)=uu=[u2uvuvv2]. (12)
    and uv=uivj,i,j=1,2. Indeed, it holds
    N(u)=((u2)x+(uv)y(uv)x+(v2)y)=F(u).

    A locally conservative DG discretization will be employed for the Navier-Stokes equation (9)-(11). We denote by Eh a shape-regular triangulation of the domain ˉΩ into triangles, where h is the maximum diameter of elements. Let ΓIh be the set of all interior edges of Eh and ΓBh be the set of all boundary edges. Set Γh=ΓIhΓBh. For any nonnegative integer r and s1, the classical Sobolev space on a domain ER2 is

    Wr,s(E)={vLs(E)||m|r,mvLs(E)}.
    We define the spaces of discontinuous functions
    V={vL2(Ω)2:EEh, v|E(W2,4/3(E))2},M={qL2(Ω):EEh, q|EW1,4/3(E)}.
    The jump and average of a function ϕ on an edge e are defined by:
    [ϕ]=(ϕ|Ek)|e(ϕ|El)|e,{ϕ}=12((ϕ|Ek)|e+(ϕ|El)|e)).
    Further, let v be a piecewise smooth vector-, or matrix-valued function at xe and denote its jump by
    [v]:=v+nE++vnE,
    where e is shared by two elements E+ and E, and an outward unit normal vector nE+ (or nE) is associated with the edge e of an element E+ (or E). The tensor product of two tensors T and S is defined as T:S=i,jTijSij.

    Let PN(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 Vh(V) and Mh(M) for integer N1:

    Vh={vL2(Ω)2:EEh, vh(PN(E))2},Mh={qL2(Ω):EEh, qPN1(E)},
    The work was motivated by the work of Girault, Riviˊere and Wheeler on discontinuous finite elements for incompressible flows presented in a series of papers [8,15]. Some projection methods [2,19] have been developed to overcome the incompressibility constraints u=0. An implementation of the operator-splitting idea for discontinuous Galerkin elements was developed in [8]. 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 Galerkin methods to investigate dynamical behavior of vortex dominated lid-driven and cylinder flows.

    The paper is organized as follows. In Section 2, a temporal discretization for the Navier-Stokes equation is presented by operator-splitting techniques, and subsequently, the nonlinearity is linearized. Both pressure and velocity 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 interior penalty discontinuous Galerkin scheme will be used in spacial discretization for those linear equations, and in Section 4, simulations of a lid-driven cavity flow up to Re=7500, and a transient flow past a circular cylinder are presented, while a numerical investigation on the Strouhal-Reynolds-number has been done, comparable to the experimental results. Finally, Section 5 concludes with a brief summary.

    2. Temporal splitting scheme

    We consider here a third-order time-accurate discretization method at each time step by using the previous known velocity vectors. Let t be the time step, M=Tt, and tn=nt. The semi-discrete forms of problem (6)-(8) at time tn+1 is

    γ0un+1α0unα1un1α2un2tνΔun+1+pn+1=β0N(un)β1N(un1)β2N(un2)+f(tn+1), (13)
    un+1=0, (14)
    which has a timestep constraint based on the CFL condition (see [14]):
    tO(LUN2),
    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 time-splitting 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 un and un1 (n1) are known, the following linearized third-order formula can be used

    γ0˜uα0unα1un1α2un2t=β0N(un)β1N(un1)+f(tn+1) (15)
    with the following coefficients for the subsequent time levels (n2)
    γ0=116, α0=3, α1=32, α2=13, β0=2, β1=1.
    Especially, by using the Euler forward discretization at the first time step (n=0), we can get a medium velocity field u1 by
    u1u0t=N(u0)+f(t1),
    and u2 by
    γ0u2α0u1α1u0t=β0N(u1)β1N(u0)+f(t2),
    which adopts the following coefficients to construct a second-order difference scheme for the time level (n=2) in (15)
    γ0=32, α0=2, α1=12, α2=0, β0=2, β1=1.
    • The second stage

    The pressure projection is as follows

    γ0˜˜u˜ut=pn+1. (16)
    To seek pn+1 such that ˜˜u=0, we solve the system
    Δpn+1=γ0t˜u, (17)
    with a Neumann boundary condition being implemented on the boundaries as
    pn+1n=fn+1β0n(unt+F(un)νΔun)β1n(un1t+F(un1)νΔun1)β2n(un2t+F(un2)νΔun2):=Gn.

    One can compute the vorticity ωn=×un at time tn=nt. Then we use pn+1 to update the intermediate velocity ˜˜u by (16).

    • The third stage is completed by solving

    γ0un+1˜˜ut=νΔun+1,

    which can be written as a Helmholtz equation for the velocity

    Δun+1+γ0νtun+1=γ0νt˜˜u.
    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 n2. 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.

    3. The spatial discretizations

    For spacial approximations, assume that piecewise polynomials of order N are employed, then the approximation space can be rewritten as Vh=Kk=1PN(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 [11,12]. We let u be approximated by uhVh and adopt a suitable approximation for the term F, i.e., F(u)F(uh), where F(uh) also can be represented as the L2-projection of F(uh) on each element of Th. Multiplying the nonlinear term by a test function vhVh, integrating over the computational domain, and applying integration by parts, we have

    Ω(F)vhdx=EkEk(F)vhdx+eΓhene[Fvh]ds, (19)
    where the term (F)vh equals to Fijvhixj, for i,j=1,2 and the indexes i,j correspond to the components of the related vectors. On each edge eE1E2 shared by two elements, to ensure the flux Jacobian of purely real eigenvalues, we may define λ+E1,e, λE2,e the largest eigenvalue of the Jacobians u(Fne)|ˉuE1 and u(Fne)|ˉuE2, respectively, where ˉuE1 and ˉuE2 are the mean values of uh on the elements E1 and E2, 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
    ne[Fvh]=ne{F}[vh]+λe2[uh][vh],
    with λe=max(λ+E1,e,λE2,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 well-conditioned Vandermonde matrix and the resulting interpolation well behaved, which greatly simplifies the formulas. The C0 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(pn+1h,ϕh)=Lp(ϕh), ϕMh,
    where
    a(pn+1h,ϕh)=EkEhEkpn+1hϕhdxekΓhek{pn+1hn}[ϕh]dsekΓhek{ϕhn}[pn+1h]ds+ekΓhσe|ek|βek[pn+1h][ϕh]ds,andLp(ϕh)=EkEhEkγ0t˜uϕhdx+ekΩekϕhGn.
    In general, σe shall be chosen sufficiently large to guarantee coercivity, more accurately, the threshold values of σe in ntee coercivity, more accurat 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 σe are presented in [20,21]. Analogously, the SIPG discretization for (18) is given by
    a(un+1h,vh)+γ0νt(un+1h,vh)Ω=Lu(vh),vhVh,
    where
    a(un+1h,vh)=EkEhEkun+1h:vhdxekΓhek{un+1h}n[vh]dsekΓhek{vh}n[un+1h]ds+ekΓhσe|ek|βek[un+1h][vh]ds,
    and
    Lu(vh)=(γ0νt˜˜uh,vh)ΩekΩek(vhneσe|ek|βvh)u0.
    where β=1. As a DG method, these SIPG schemes have 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 polynomialsmight 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 t1E03 are taken.

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

    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 three times with piecewise discontinuous elements being applied into the fully discrete SIPG approach.

    Employing the uniform mesh (see the left profile in Figure1) and approximation polynomials of order 3, we illustrate the velocity vector, pressure and vorticity profiles in Figure2, observing that the main characterization of solutions can be captured by fine meshes well, except a few very small oscillations occurred close to the upper boundary along the lid for computing velocity in ydirection and pressure. Oscillations on the upper lid propagates from a velocity singularity that exists at the corners. 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 adaptable to adopt an adaptive meshes for solving those singularity problems. Here, we apply the semi-implicit SIPG method in a locally refined mesh (see the right profile in Figure1) to solve the incompressible flow. In Figure3, the velocity profiles of (u,v) through the geometric center of the cavity are plotted with Re=100,400,1000,5000,7500 taken. From Figures 4-8, 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. 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.

    Figure 1. #1: A uniform mesh (left); #2: An initial locally refined mesh (right).
    Figure 2. The profiles of velocity components (u,v), pressure and vorticity for Re=2000 by using mesh #1 and N=4.
    Figure 3. Velocity profiles (u,v) through geometric center of the cavity for Re=100,400,1000,5000,7500.
    Figure 4. Re=100 and N=3, mesh #2. Left: pressure contour; Right: vorticity contour.
    Figure 5. Re=400 and N=3, mesh #2. Left: pressure contour; Right: vorticity contour.
    Figure 6. Re=1000, N=4, mesh #2. Left: pressure contour; Right: vorticity contour.
    Figure 7. Re=5000, N=3, mesh #2 with a refinement once. Left: pressure contour; Right: vorticity contour.
    Figure 8. Re=7500, N=3, a refined mesh of #2. Left: pressure contour; Right: vorticity contour.

    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 inflow boundary condition is (u,v)=(1,0), while the outflow boundary is un=0. To the boundary conditions on the upper and lower sides, we present two different conditions for comparison (see Figure9), which are wall (u=0,v=0) and homogeneous Neumann boundary conditions (u=1,vn=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 ρ=1 kgm3 and a locally refined mesh (maxh=0.088) will be used for the simulations.

    Figure 9. Flow past a cylinder and N=3. Top: A channel with walls; Bottom: A channel without a wall. The free stream velocity on the inflow boundary is u=(1,0).

    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 [13], 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 Fig.10, 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 cylinder's 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 [1] environments topologically similar to Fig.10 with a constant mass flow rate and no-slip walls. Drela [3] develops an analysis for 2D wind tunnels resulting in an effective coefficient of drag of

    Cd=(112cHCdπ2AH2)Cdun,
    and an effective Reynolds number of
    Re=(1+14cHCd+π6AH2)Reun,
    where the un subscript represents the uncorrected value, H represents the domain height, A represents the cylinder area and c represents the cylinder radius. Drela's analysis does not specifically include the boundary layer forming on the 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/H2) effects.

    Figure 10. Flow past a cylinder in a channel with walls, Re=100. Top: Vorticity contour; Bottom: Pressure contour.

    Alternatively, the flow without walls in Figure11 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 Figure12, 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 xdirection would reduce dispersively, in other words, the vortex lifespan is less than those produced in the channel with walls. The velocity profiles in the ydirection have been given in Figure13.

    Figure 11. Flow past a cylinder in a channel without a wall, Re=100. Top: Vorticity contour; Bottom: Pressure contour.
    Figure 12. Flow past a cylinder in a channel, Re=100. Top: Velocity in x-direction with walls; Bottom: Velocity in x-direction without a wall.
    Figure 13. Flow past a cylinder in a channel, Re=100. Top: Velocity in y-direction with walls; Bottom: Velocity in y-direction without a wall.

    We localize the domain around the cylinder and refine the mesh, then show the vorticity startup behavior in Figure14 as well as the pressure Figure15. Upon startup, two vortices of opposite direction are formed on the upper and lower 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 (Figure16 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.

    Figure 14. Vorticity contours of flow past a cylinder without a wall at different time.
    Figure 15. Pressure contours of flow past a cylinder without a wall at different time.
    Figure 16. Drag coefficients, 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 and the diameter of the cylinder D=0.1, we will chose different viscosity coefficients ν=1e3,5e4,2.5e4 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

    dp=p(t;0.05,0)p(t;0.05,0).
    We use the definition of Cd and Cl given in [17] as follows:
    Cd=2ρDu2maxS(ρνutS(t)nnyp(t)nx)ds,andCl=2ρDu2maxS(ρνutS(t)nnx+p(t)ny)ds,
    where n=(nx,ny)T is the normal vector on the cylinder boundary S directing into Ω, tS=(ny,nx)T the tangential vector and utS the tangential velocity along S. In the literature, Fey et al. in [6] propose the Strouhal number represented by piecewise linear relationships of the form
    St=ST+mRe
    with different values St and m in different shedding regimes of the 3D circular cylinder wake. We apply the periodic Tlift of the lift coefficients (see Figure16) to express the periodic T:=1f appearing in the definition of Strouhal number, i.e.,
    St(Re)=DTliftu, (20)
    which comes from the classical definition (5). From the evolution of Cd, Cl and dp as in Figure16, we may find a period Tlift of the lift coefficients for different values of Re to calculate Strouhal number by (20). In Figure17, a comparison of Strouhal numbers between the experimental estimates in Fey etc. [6] and our estimates from (20) 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.

    Figure 17. A comparison of Strouhal-Reynolds-number between our estimate and the linear fit in [6] for 50Re400.

    5. Conclusion

    This paper developed a Symmetric Interior Penalty discontinuous Galerkin numerical solver for the incompressible Navier Stokes equations of fluid flow, with the temporal splitting technique applied in decoupling the diffusion and pressure terms and the local Lax-Friedrichs flux for the DG discretization of the nonlinear convection term. Two testcases are presented: a lid-driven cavity and a cylinder flow. The SIPG 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.

    Acknowledgments

    The authors would like to thank the anonymous referee and editor very much for their valuable comments and suggestions, which greatly help us improve the presentation of this article. The work of the first author was partially supported by the Natural Science Foundation of Gansu Province, China (Grant 145RJZA046), and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase).

    Conflict of Interest

    We declare no conflicts of interest in this paper.

    [1] J. B. Barlow, W. H. Rae, and A. Pope (1999) Low-Speed Wind Tunnel Testing, John Wiley.
    [2] A. J. Chorin (1969) On the convergence of discrete approximations to the Navier-Stokes equations, Math. Comp., 23.
    [3] M. Drela (2014) Flight Vehicle Aerodynamics.MIT Press, Boston .
    [4] Y. Epshteyn, B. Rivi`ere (2007) Estimation of penalty parameters for symmetric interior penalty Galerkin methods,.J. Comput. Appl. Math., 206 .
    [5] U. Fey, M. K¨onig, and H. Eckelmann (1998) A new Strouhal-Reynolds-number relationship for the circular cylinder in the range 47 ≤ Re ≤ 2 × 105,.Physics of Fluids 1547-1549.
    [6] U.Ghia, K. N. Ghia, C. T. Shin (1982) High-Re solutions for incompressible flow using the Navier-Stokes equations and a multigrid method..J. Comput. Phys., 48: 387-411.
    [7] V. Girault, B. Rivi`ere, and M. F. Wheeler (2005) A splitting method using discontinuous Galerkin for the transient incompressible Navier-Stokes equations,.ESAIM: Mathematical Modelling and Numerical Analysis, 39: 1115-1147.
    [8] V. Girault, B. Rivi`ere, and M. F. (2005) Wheeler, A discontinuous Galerkin method with non-overlapping domain decomposition for the Stokes and Navier-Stokes problems,.Math. Comp 53-84.
    [9] O. Goyon (1996) High-Reynolds number solutions of Navier-Stokes equations using incremental unknowns,.Comput. Method. Appl. M.130 319-335.
    [10] J. S. Hesthaven (1998) From electrostatics to almost optimal nodal sets for polynomial interpolation in a simplex,.SIAM J. Numer. Anal. 35: 655-676.
    [11] J. S. Hesthaven, C. H. Teng (2000) Stable spectral methods on tetrahedral elements, SIAM.J. Sci. Comput., 2352-2380.
    [12] S. F Hoerner (1965) Fluid-Dynamic Drag, Hoerner Fluid Dynamics.Bakersfield .
    [13] G. Karniadakis, S. J. Sherwin (2005) Spectral/hp element methods for CFD, Oxford University Press.New York .
    [14] S. Kaya, B. Rivi`ere (2005) A discontinuous subgrid eddy viscosity method for the time-dependent Navier-Stokes equations, SIAM J.Numer. Anal 43: 1572-1595.
    [15] B. Rivi`ere, M. F. Wheeler, and V. Girault (1999) Improved energy estimates for interior penalty, constrained and discontinuous Galerkin methods for elliptic problems..Part I, Comput. Geosci., 337-360.
    [16] M. Sch¨afer, S. Turek (1996) The benchmark problem ‘flow around a cylinder’, In Flow Simulation with High-Performance Computers II, Hirschel, E.H.(ed.). Notes on Numerical Fluid Mechanics, vol. 52.Vieweg, Braunschweig, 547-566.
    [17] J. Shen (1991) Hopf bifurcation of the unsteady regularized driven cavity flow,.J. Comput. Phys 95: 228-245.
    [18] J. Shen (1996) On error estimates of the projection methods for the Navier-Stokes equations: Secondorder schemes,.Math. Comp 65: 1039-1065.
    [19] L. Song, Z. Zhang (2015) Polynomial preserving recovery of an over-penalized symmetric interior penalty Galerkin method for elliptic problems,.Discrete Contin. Dyn. Syst. – Ser. B 20: 1405-1426.
    [20] L. Song, Z. Zhang (2015) Superconvergence property of an over-penalized discontinuous Galerkin finite element gradient recovery method,.J. Comput. Phys 299: 1004-1020.
    [21] R. Temam (2001) Navier-Stokes Equations: Theory and Numerical Analysis, AMS Chelsea publishing.Providence .
    [22] R. Temam (1995) Navier-Stokes Equations and Nonlinear Functional Analysis, Volume 66 of CBMSNSF Regional Conference Series in Applied Mathematics.SIAM, Philadelphia, second edition .
    [23] M. F. Wheeler (1978) An elliptic collocation-finite element method with interior penalties,.SIAM J. Numer. Anal., 15: 152-161.
  • This article has been cited by:

    1. Lunji Song, 2020, Chapter 4, 978-1-83962-616-6, 10.5772/intechopen.94316
  • Reader Comments
  • © 2016 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(7208) PDF downloads(1797) Cited by(1)

Article outline

Figures and Tables

Figures(17)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog