# Simulation of Vortex Shedding behind a Flat Plate with Vorticity Based Adaptive Spectral Element Method

Authors: Li-Chieh Hsu; Guo-Jhih Gao
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Vorticity error based adaptive meshes refinement scheme is developed and employed using spectral element method to simulate flow past object problems. In general, it is hard to predict and enhance meshes effectively in a region where  the error is larger in the computational domain by using the conforming mesh method. Employing finer meshes throughout the whole domain leads to lengthy computational time and excessive storage. Therefore, an indicator is used to predict the regions where larger errors exist and mesh refinement is needed. To compare the efficiency of indicators, three kinds of properties are used as mesh refinement indicators, including the synthesis of velocity and pressure estimated error, vorticity estimated error, and estimated error decay rate. Simulations of the cavity flow in Re = 100 and 1000 and the cases of flow past an inclined flat plate in Re = 100 to 1000 are performed with the adaptive mesh method and conforming mesh method. The results show that the adaptive mesh method can provide the same accuracy as that of the conforming mesh method with only 62% of the elements.


---

## Body

## 1. Introduction

Adaptive refinement techniques provide attractive approaches to many classes of problems in computational fluid dynamics. The goal is to produce the best numerical simulations of a complex flow with the least computational effort. Adaptive refinement methods for finite element methods have been proposed by Babuška and Guo [1] and Rachowicz et al. [2]. Basically, there are two major kinds of mesh refinement methods, h-refinement in which the elemental mesh is refined and p-refinement in which the local degree of the polynomial shape function is increased.The spectral element method [3] is a high-order weighted residual technique for numerically solving partial differential equations that combines the generality of h-type finite element methods [4, 5] and the rapid convergence rate of p-type spectral methods [6]. The method breaks up complex computational domains into some macrosubdomains called elements upon which trial functions and solutions are represented by Nth order tensor-product orthogonal polynomial expansions to promise high accuracy. The minimal continuity condition C0 is imposed on the interfaces between elements due to the requirement of variational procedures. Variational projection operators and Gauss numerical quadrature are used in discrete equations which are solved by iteration procedures with tensor-product sum-factorization techniques [7, 8].Although this method is general and flexible enough to handle complex geometries [9], its practical implementation is limited by the arbitrary choice of domain discretization. The nonconforming formulation [10] allows globally unstructured grids whereby elements may match up in an arbitrary manner by means of a mortar space, a set of functions, which mortar the brick-like spectral elements together. The nonconforming spectral element requires only minimum L2 errors on the interfaces between elements as opposed to the C0 continuity required by the conforming spectral element method as shown Figure 1. This means that this method allows a larger element to merge with several smaller elements arbitrarily on a particular face as in Figure 2. Although the nonconforming spectral element method breaks the limitation of the conforming spectral element method, a proper scheme which can precisely indicate the necessary regions to be refined and a strategy to get faster convergence rates is still needed.Figure 1
Conforming spectral element, 4 elements.Figure 2
Nonconforming spectral element, 3 elements.Several adaptive mesh refinement algorithms have been incorporated with the finite element method, finite volume method, and spectral element method to increase computational efficiency. Berrone et al. [11] employed the adaptive finite element method and finite volume method, separately, to simulate the flow past a rectangular cylinder. The Galerkin least square finite element method is used in space discretization. A semi-implicit Euler method is applied in time discretization. Two posterior error indicators are used to control the mesh distribution and time length of each step, respectively. The space error indicator is related to the stress-jump on the internal edge and residuals in the Neumann boundary conditions. The time error indicator is based on the square magnitude of the solution gradient in a time interval for deciding the time length of a step. This study shows that the adaptive finite element method can solve problems accurately and make good agreements with solutions obtained by the finite volume method with respect to Strouhal number, drag coefficient, and recirculation length. Galvão et al. [12] used the hp-type adaptive least-squares spectral element method to solve hyperbolic partial differential equations. For nonconforming neighboring elements, the mortar element method is applied to patch the interface between elements. In their study, the integral condition was applied to ensure weak continuity across the nonconforming interface. The vertex condition ensures that the solutions agree at the end points of the mortar. Dupont and Lin [13] used the discontinuous triangular spectral element method to simulate a shallow water problem to compare their results with those obtained by second-order and fourth-order finite difference methods as well as the finite element method. The Legendre polynomials are chosen as the basis functions for solution expression inside the elements. Fourth-order Runge-Kutta method is used for time integration. The results of the spectral element method show its superiority to finite element methods in energy conservation with no dissipation and higher accuracy within the same CPU time. However, though the adaptive spectral element method is better than the spectral element and finite difference methods, the accuracy-to-cost convergence is not as good as that achieved by the spectral element method. A tree-based adaptive mesh refinement (AMR) algorithm for the high-order discontinuous Galerkin method on quadrilateral grids with nonconforming elements is proposed by Kopera and Giraldo [14] to solve the atmospheric simulations which require highly expensive large computational cost. In that study, the rising thermal bubble case and the general features of the solution are captured correctly by AMR simulations with 15 times speed-up.For the spectral element method, there are two major mesh refinement approaches, includingh-type and p-type refinements. In h-refinement, computational grids are enhanced by increasing number (K) of elements. The p-refinement is based on increasing the order (N) of polynomial of the basis function within an element. Mavriplis and Hsu [15] developed an adaptive h-refinement algorithm for incompressible flow simulation with the spectral element method. Feng and Mavriplis [16] further used p-type and h-type with adaptive mesh refinement and mesh coarsening techniques to capture a 2D thin flame front, successfully. Later, Feng et al. [17] continued to extend this approach to the 3D mortar element method for nonconforming mesh refinement. Henderson [18] developed another nonconforming element method, namely, the piecewise mortar element method by incorporating dynamic meshes for 2D incompressible flow simulation. Further, the nonconforming element method has also been extended to be three-dimensional scheme with the Fourier expansion technique for turbulence flow simulation. Recently, the mesh refinement is also favorable in high-order computation. Premasuthan et al. [19] applied the spectral difference (SD) method and local mesh refinement on high-order computation for compressible fluid flows with discontinuities such as supersonic flow past a bump, transonic flow past NACA0012 airfoil, and supersonic flow past a cylinder. In their study, the mortar element method is applied on the communication of the interfaces between nonconforming elements by projection matrices which are responsible project subdomains solutions to mortar space and project the computed fluxes from mortar back to subdomains.For mesh adaptation, several physical properties can be used as the indicators of mesh refinement or coarsening. For example, the discontinuity of density and velocity are most commonly used in compressible flow. The derivatives of velocity or pressure are usually used in smooth solutions such as for incompressible flow. Oh et al. [20] solved Navier-Stokes equations with the finite difference method to simulate vortex shedding behind a circular cylinder (Re=200) as well as the transonic flow past an airfoil by applying adaptive meshes with refinement and coarsening functions for which the error indicator was composed of the truncation error of the magnitude of the density gradient and density change of the time rate. Rausch et al. [21] proposed using the absolute value of the substantial derivative of density as an enrichment indicator. Hwang and Kuo [22] introduced a new mesh-enrichment indicator which is based on the gradient of the substantial derivative of density.In this study, an error indicator which is based on error estimation, a posteriori error, is used to decide which region needs to be refined. In the past, estimated errors of the primary properties such as pressure, velocity, and even the decay rate of the sequence of coefficients in the Legendre series expansion of the estimated error were selected as the mesh refinement indicators. However, the mesh refinement quality and solution accuracy are not stable enough when using these indicators, especially for cases of external flow. Therefore, in this study, we present a vorticity error based nonconforming adaptive spectral element method in which meshes are automatically refined with time in order to capture the underresolved regions of the domain and follow regions requiring high resolution. To compare the accuracy and efficiency, the estimated error of velocity, pressure, vorticity, and the decay rate of the sequence of coefficients in the Legendre series expansion of estimated error are selected as the error indicators to perform mesh refinement.This nonconforming method is applied in the simulation of cavity flow, flow past an inclined flat plate, and flow past a cylinder for comparison with the results of the conforming method.


## 2. Numerical Method

Our objective is to simulate complex flows in the laminar regimes by solving the incompressible Navier-Stokes equations:(1)∂u⃑∂t+u⃑·∇u⃑=-1ρ∇p+ν∇2u⃑,∇·u⃑=0,
where u⃑ is the velocity vector, p is the pressure, ρ is the density, and ν is the kinematic viscosity. A time splitting scheme [23] is introduced to treat the nonlinear convective terms explicitly while the pressure and velocity diffusion terms can be solved implicitly as follows:
(2)u⃑^-u⃑nΔt=∑r=02βr(-u⃑·∇u⃑)n-r,(3)u⃑^^-u⃑^Δt=1ρ∇pn+1,(4)u⃑n+1-u⃑^^Δt=ν∇2u⃑n+1.
The scheme uses intermediate time step values of velocity, u⃑^ and u⃑^^, where ∧ and ∧^ are between the nth and (n+1)st time steps. The nonlinear terms are advanced by a third order Adams-Bashforth scheme with coefficients βr. The pressure term can be treated as an elliptic problem separately by taking divergence on both sides of (3):
(5)∇2pn+1=ρΔt∇·u⃑^.
Hence, the pressure and diffusion terms are modeled as Helmholtz problems which are discretized by the spectral element method as shown below. A general model of the Helmholtz equation is shown as follows:
(6)(∇2-λ2)ϕ=g,on⁡Ω.
The variational form of (6) with an arbitrary trial function ψ=δϕ is derived as
(7)-∫Ω∇ϕ·∇ψdx-λ2∫Ωϕψdx=∫Ωgψdx,∀ψ.
The computational domain is broken up into macroelements as shown in Figure 1. All variables can be approximated by discretization as follows:
(8)ϕhk=∑p=0N∑q=0Mϕpqkhp(r)hq(s),
where the h functions are the Lagrangian interpolants based on the orthogonal set of Legendre polynomials of high degree N (or M) and r and s are the local coordinates on each element k. By performing Gauss-Lobatto quadrature and summing contributions from adjacent elements, the global matrix equation is obtained as follows:
(9)(A-λ2B)ϕ=Bg,
where A is the discrete Laplacian operator and B is the mass matrix. Here, the global matrix equation is solved by the preconditioned conjugate gradient. Details can be found in a previous study [24].


### 2.1. Nonconforming Element Method

There are two methods for constructing nonconforming discretizations for elliptic problems. The first is the iterative method originally proposed by Funaro et al. [25] for collocation schemes to patch subdomains with a C1 condition across the interfaces. This algorithm has been modified and incorporated within the variational framework of spectral elements [26, 27]. The second method is based on the standard spectral element method, namely, the nonconforming spectral element or mortar element method, as introduced by Maday et al. [10, 28].The nonconforming spectral element provides a technique which only requires minimumL2 errors on the interfaces between elements rather than C0 continuity, as required by the conforming spectral element method. This means that the nonconforming spectral element method allows a larger element to merge with several smaller elements arbitrarily on a particular face. Nonconforming discretizations greatly improve the flexibility of the spectral element approach with regard to automatic mesh generation and local mesh refinement while preserving the convergence properties of spectral element discretization. Standard spectral element discretization and interpolations are employed within local subdomains, whereas a mortar space facilitates communication at the interfaces between nonconforming subdomains. With projection operators, the mortar space collects the information from local subdomains and sums them together at each collocation point in the mortar space. The composite solutions are fed back to local subdomains by inverse projection operators.There are some accuracy losses once the meshes become nonconforming. To solve this problem, the consistency error is kept small through the combination of theL2 condition and vertex condition. The L2 condition, (12), minimizes the jump of functions in the interior of the boundaries of the elements. The vertex condition, (11), ensures exact continuity at cross points. The integral matching enables exponential convergence. Mortar data are transferred into the elemental edge data by a projection operator matrix.Mortar space,Wh, is defined as
(10)Wh={ϕ∈C0(S),∀p=1,…,M,ϕ|γpWh=ϕ|γp∈PN(γp),ϕ|∂Ω=0},
where PN(γp) is the space of all polynomials on γp of degree ≤N. The nonconforming space for solutions on the elemental edges is given as Xh=v∈L2(Ω), ∀k=1,…,K, v|Ωk∈PN(Ωk), such that ∃∅∈Wh for which
(11)∀q=1,…,V,∀k=1,…,K,suchthatvqisavertexofΩk,v|Ωk(vq)=ϕ(vq),(12)∀l=1,…,4,∀k=1,…,K,∀ψ∈PN-2(Γk,l),∫Γk,l(v|Ωk-ϕ)ψds=0,
where M is the total number of mortar elements, V is the total number of vertices, and K is the total number of elements. Mortar function ϕ|γp is simply restricted to the edge of an element. Therefore, it can be expressed by one-dimensional Gauss-Lobatto-Legendre (GLL) interpolants:
(13)ϕ|γp=∑j=0NϕjPhjN(s~),∀p∈{1,…,M}.
The discrete mortar function and local elemental boundary solution and basis are applied to the L2 condition, (12). The L2 and vertex conditions can be represented in discrete forms by choosing the basis function, ψ, as follows:
(14)ψ|Γk,l=∑q=1N-1βqηqN-2(s^),
where
(15)ηqN-2(z)=(-1)N-qLN′(z)εq-z,z∈|-1,1|,q∈{1,…,N-1}.
Mortar data are transferred into the elemental edge data by a projection operator matrix Q which is derived with the basis function as (14) and (15). For instance, a standard Poisson equation, the conforming discretized global matrix equation, is
(16)Au=Bf.
The corresponding system in the nonconforming case is
(17)QTAQu=QTBf.Q is a global projection operator matrix which projects the solutions from mortars to the interior points of the local edge of elements, whereas QT is also a global projection operator matrix which projects solutions from the local edge of elements to mortars. The detailed derivation can be found in the previous works [10, 28]. Similar work about mortar method on treatment of sub-domain refinement also can be found in the study of Kopriva [29] who incorporated mortar method with the spectral difference method [30] to solve compressible flow.


## 2.1. Nonconforming Element Method

There are two methods for constructing nonconforming discretizations for elliptic problems. The first is the iterative method originally proposed by Funaro et al. [25] for collocation schemes to patch subdomains with a C1 condition across the interfaces. This algorithm has been modified and incorporated within the variational framework of spectral elements [26, 27]. The second method is based on the standard spectral element method, namely, the nonconforming spectral element or mortar element method, as introduced by Maday et al. [10, 28].The nonconforming spectral element provides a technique which only requires minimumL2 errors on the interfaces between elements rather than C0 continuity, as required by the conforming spectral element method. This means that the nonconforming spectral element method allows a larger element to merge with several smaller elements arbitrarily on a particular face. Nonconforming discretizations greatly improve the flexibility of the spectral element approach with regard to automatic mesh generation and local mesh refinement while preserving the convergence properties of spectral element discretization. Standard spectral element discretization and interpolations are employed within local subdomains, whereas a mortar space facilitates communication at the interfaces between nonconforming subdomains. With projection operators, the mortar space collects the information from local subdomains and sums them together at each collocation point in the mortar space. The composite solutions are fed back to local subdomains by inverse projection operators.There are some accuracy losses once the meshes become nonconforming. To solve this problem, the consistency error is kept small through the combination of theL2 condition and vertex condition. The L2 condition, (12), minimizes the jump of functions in the interior of the boundaries of the elements. The vertex condition, (11), ensures exact continuity at cross points. The integral matching enables exponential convergence. Mortar data are transferred into the elemental edge data by a projection operator matrix.Mortar space,Wh, is defined as
(10)Wh={ϕ∈C0(S),∀p=1,…,M,ϕ|γpWh=ϕ|γp∈PN(γp),ϕ|∂Ω=0},
where PN(γp) is the space of all polynomials on γp of degree ≤N. The nonconforming space for solutions on the elemental edges is given as Xh=v∈L2(Ω), ∀k=1,…,K, v|Ωk∈PN(Ωk), such that ∃∅∈Wh for which
(11)∀q=1,…,V,∀k=1,…,K,suchthatvqisavertexofΩk,v|Ωk(vq)=ϕ(vq),(12)∀l=1,…,4,∀k=1,…,K,∀ψ∈PN-2(Γk,l),∫Γk,l(v|Ωk-ϕ)ψds=0,
where M is the total number of mortar elements, V is the total number of vertices, and K is the total number of elements. Mortar function ϕ|γp is simply restricted to the edge of an element. Therefore, it can be expressed by one-dimensional Gauss-Lobatto-Legendre (GLL) interpolants:
(13)ϕ|γp=∑j=0NϕjPhjN(s~),∀p∈{1,…,M}.
The discrete mortar function and local elemental boundary solution and basis are applied to the L2 condition, (12). The L2 and vertex conditions can be represented in discrete forms by choosing the basis function, ψ, as follows:
(14)ψ|Γk,l=∑q=1N-1βqηqN-2(s^),
where
(15)ηqN-2(z)=(-1)N-qLN′(z)εq-z,z∈|-1,1|,q∈{1,…,N-1}.
Mortar data are transferred into the elemental edge data by a projection operator matrix Q which is derived with the basis function as (14) and (15). For instance, a standard Poisson equation, the conforming discretized global matrix equation, is
(16)Au=Bf.
The corresponding system in the nonconforming case is
(17)QTAQu=QTBf.Q is a global projection operator matrix which projects the solutions from mortars to the interior points of the local edge of elements, whereas QT is also a global projection operator matrix which projects solutions from the local edge of elements to mortars. The detailed derivation can be found in the previous works [10, 28]. Similar work about mortar method on treatment of sub-domain refinement also can be found in the study of Kopriva [29] who incorporated mortar method with the spectral difference method [30] to solve compressible flow.


## 3. Error Estimator

The numerical solution can be represented by the sum of finite terms of Legendre polynomials. Because the solution is expressed by finite terms of polynomials, truncation error inL2 norm can be represented as (18) for one dimension and (19) for two dimensions, respectively:
(18)∥u~-uh∥=(∑n=N+1∞a~n2(2n+1)/2)1/2,(19)∥u~-uh∥=(∑n=N+1∞∑m=M+1∞a~nm2(2n+1)(2m+1)/4)1/2,
where u~ is the exact solution, uh is a numerical solution, and an is a Legendre polynomials coefficient. In addition, the numerical solution often diverges from the best polynomial fit. This is referred to as quadrature error. The an can be obtained exactly by GLL quadrature, for n≦N-1, whereas it is just approximated, for all n≧N. Therefore, the approximation error is shown as follows:
(20)∥uh-∏hu~∥=(2(bN)22N+1)1/2.bN can be approximated by aN. The total estimated error is the sum of the truncation error, (18), and quadrature error, (20), as shown below:
(21)εestk=((aNk)2(2N+1)/2+∫N+1∞(ank(n))2(2n+1)/2dn)1/2.
For n=N+1, ∞, there is no solved solution, uh, for calculation of the Legendre polynomials coefficient, an. Therefore, it is assumed that an follows a model which decays exponentially as the polynomial order, N, increases:
(22)an=Ce-σn.
In each dimension, the known last four coefficient points, an, n=N-3, N, are used with the least squares best fit method to extrapolate the coefficient, C, and the decay rate σ. The decay rate, σ, indicates the quality of resolution. Generally, the result is regarded as having poor resolution for σ<1 or good resolution for σ>1. In the adaptive process, this information is used to choose refinement options: σ<1, increasing the number of elements or, for σ>1, increasing the order of the polynomial N. In this study, estimated errors of velocity, pressure, and vorticity as well as estimated error coefficient decay rate σ are used as the mesh refinement indicators to compare the accuracy and efficiency. Among them, once the estimated error coefficient decay rate σ is used as an indicator, only the elements with σ<0.8 are chosen for mesh refinement.Error estimators are based on a posteriori estimates [31] of the L2 errors in the spectral representation of the solution for each element to determine which regions should be refined. Once the new grid has been defined after mesh refinement, the current (old time) step solution must be interpolated onto the new topology. The coordinates of the new elements and their Gauss-Lobatto collocation points are determined by the old mesh with spectral distribution. The solutions on the new Gauss-Lobatto collocation points inside the new elements are evaluated by the Lagrangian interpolation from the old mesh. In short, the solution strategy is as follows: compute an initial solution with a suitable initial mesh; estimate errors in the solution locally for each element; refine or coarsen the mesh according to the error estimators; interpolate old mesh solutions onto the new elements; and resume the numerical solution process.


## 4. Results and Discussion

### 4.1. Validation

The case of laminar incompressible flow past a square cavity is commonly used to validate numerical methods. There are discontinuities in velocity on the upper corners, whereu=0 derives from the no-slip condition on the vertical walls and u=1 results due to the imposed driving flow on the upper boundary. Hence, there are singularities on these two corners. The cavity flow regime from Reynolds number 100 to 1000 has been investigated by Ghia et al. [32]. Ghia’s results are used for comparison with our results of conforming meshes and adaptive meshes. The conforming approach with 64 equal-sized elements and an order of shape function of 11 is used for simulation. The nonconforming approach is performed with the initial mesh which has four equal-sized elements with shape function order, N=11.Three mesh refinement indicators which are the synthesis of estimated error of velocity and pressure, estimated error of vorticity, and estimated error decay rate (σ) are employed for the sake of comparison. At each refinement step, the element with the largest error is applied to conduct mesh refinement. After mesh adaptation, the mesh becomes 34 nonconforming elements. The meshes and results are shown in Figure 3 for Re=100 and Figure 4 for Re=1000, respectively. There are two obvious secondary vortices in the bottom corners, especially in the case of Re=1000. The profile of y-directional velocity, v, along the horizontal line through geometric center and the profile of x-directional velocity, u, along the vertical line through the geometric center are used for comparison with the result of multigrid method [32]. The results from all conforming and nonconforming approaches show good agreement with those in [32], as shown in Figures 5 and 6 for Re=100 and Figures 7 and 8 for Re=1000, respectively. From the above cases, we can see that these adaptive nonconforming schemes not only capture the physical phenomena precisely with 53% elements of that of the conforming method but also attain the same accuracy as achieved by the conforming meshes approach and other numerical methods. The global L2 error decays exponentially as the element number is increased which fits the property of h-type mesh refinement as shown in Figure 9.Figure 3
Re
=
100 streamlines, nonconforming approach by vorticity error base.Figure 4
Re
=
1000 streamlines, nonconforming approach by vorticity error base.Figure 5
Re
=
100, v-profile along central line of x-span in square cavity.Figure 6
Re
=
100, u-profile along central line of x-span in square cavity.Figure 7
Re
=
1000, v-profile along central line of x-span in square cavity.Figure 8
Re
=
1000, u-profile along central line of x-span in square cavity.Figure 9
Global error decay in adaptive mesh refinement, cavity flow,Re=1000.


### 4.2. Flow Past a Flat Plate

Flow past a thin airfoil is an interesting topic for science and engineering because it abounds with rich physical phenomena such as shear layers, leading and trailing edge vorticity, laminar separation, and vortex shedding. It is also important for microaerial vehicles (MAVs) design.Lam and Leung [33] investigated experimentally the shedding of vortices from an inclined flat plate. The results showed that the trailing edge vortex in the separated wake possesses more intense vorticity levels but a smaller vortex size than the leading edge vortex. Breuer and Jovičić [34] conducted a large eddy simulation (LES) for the flow past an inclined flat plate in Re=20000, at an angle of attack of 18°. The findings revealed Kelvin-Helmholtz instability in the free shear layer behind the leading edge. A highly asymmetric wake with vortices of unequal strength was observed in the wake region.However, there have been few studies of flow past a flat plate in hundreds of Reynolds number. Taira et al. [35] studied the flow past an inclined thin flat plate with immersed boundary method. The numerical results of flow Re=100 past a three-dimensional flat plate with an aspect ratio=2 agree with experimental lift and drag coefficients in various angles of attack. The effects of aspect ratio and planform geometry on the leading-edge vortex were studied with Re=300 flow. The result showed that the elliptic and semicircular plates generate a wake with less distinct transition between leading edge and tip vortices than that of rectangular profile. Hence, the flow control for flat plate in low Reynolds numbers has aroused the attentions of several researchers. For instance, Taira et al. [36] simulated a synthetic jet type actuator to modify the dynamics of wake vortices and the corresponding forces exerted on wings. Brunton and Rowley [37] further used the immersed boundary method, a direct numerical method, and proper orthogonal decomposition/Galerkin projection to simulate stationary flat plate in Re=100, 300 flow with angles of attack, α=30° and α=45°, and the full dynamic wake structures of a pitching or plunging flat plate.In our study, for constructing quadrilateral spectral elements on the surface of a flat plate, an extremely thin two-dimensional inclined flat plate with maximum thickness to chord ratio, 1 × 10−5, is used to approach a theoretical straight-line-like flat plate profile. Hence, there are very sharp leading edge and trailing edge which are prone to triggering some critical physical phenomena such as leading edge vortex and trailing edge vortex. A rectangular global computational domain is selected as in Figure 10. The chord of the plate is 1 unit (D) with an angle of attack 30°. The domain is extended from the center of the plate to upstream boundary with distance of 20 D and to downstream boundary with distance of 30 D. The distances from upper and lower boundaries to the plate center are both 20 D. The uniform flow boundary conditions (u=U0, v=0) are employed in the inflow and upper and lower boundaries so as to set different flow conditions, Re=100, 140, and 200. The outflow boundary condition is employed in the most downstream boundary. There are 395 macrospectral elements (K), employed for conforming approach. The order (N) of the basis function, Gauss-Lobatto-Legendre polynomial, is 11.Figure 10
Flat plateα=30°, K=395, conforming meshes.In the nonconforming approach, the initial coarse mesh with 95 elements used for simulation leads to unsatisfactory result. This result is used as the initial input of mesh adaption process. Figure11 shows the results of the mesh refinement after carrying out mesh adaptation with the indicator of vorticity estimated error 50 times. The mesh refinement regions are concentrated in the area around the flat plate, wake region, and upwind region. The refinement regions do not vary much for these three Reynolds numbers. The results of adaptive nonconforming meshes based on vorticity error are shown in Figures 13, 15, and 17. These results not only show good agreement with those counterparts done by conforming meshes as shown in Figures 12, 14, and 16 but also can define clearer vorticity contours.Figure 11
Flat plateα=30°, K=245, nonconforming meshes refinement started from initial meshes (K=95) by vorticity error indicator.Figure 12
Conforming approach with 395 elements, vorticity,Re=100, α=30°.Figure 13
Nonconforming approach by vorticity error indicator,Re=100, α=30°, vorticity with 245 elements.Figure 14
Conforming approach with 395 elements, vorticity,Re=140, α=30°.Figure 15
Nonconforming approach by vorticity error indicator,Re=140, α=30°, vorticity with 247 elements.Figure 16
Conforming approach with 395 elements, vorticity,Re=200, α=30°.Figure 17
Nonconforming approach by vorticity error indicator,Re=200, α=30°, vorticity with 249 elements.Among those mesh refinement regions, more intensive mesh refinements are executed around the flat plate, especially in the leading edge and the trailing edge where several tiny elements are employed, automatically, as shown in Figures19, 21, and 23. The result shows that vorticity, the first derivative property, can be solved smoothly by the adaptive nonconforming method. This mesh distribution indicates that the vorticity error occurs in the stronger shear layer area where the vorticity gradient is also larger. In addition, since the sharp leading edge causes more difficulties in computation, it attracts more mesh refinement activities. The results from adaptive nonconforming meshes show good agreements with those from conforming meshes as shown in Figures 18, 20, and 22. However, the element number, 247, required by adaptive nonconforming approach is much less than the 395 elements used in the conforming approach. Hence, only 62% of the element numbers are needed when adaptive mesh refinement is employed. This reveals that adaptive mesh approach can compute more efficiently than the conforming approach does, if meshes are employed in the right locations by a proper mesh refinement indicator.Figure 18
Conforming approach with 395 elements, closer look vorticity contours,Re=100, α=30°.Figure 19
Nonconforming approach with 247 elements, closer look vorticity contours,Re=100, α=30°.Figure 20
Conforming approach with 395 elements, closer look vorticity contours,Re=140, α=30°.Figure 21
Nonconforming approach with 247 elements, closer look vorticity contours,Re=140, α=30°.Figure 22
Conforming approach with 395 elements, closer look vorticity contours,Re=200, α=30°.Figure 23
Nonconforming approach with 249 elements, closer look vorticity contours,Re=200, α=30°.The clockwise leading edge vortex first forms and grows. Later, it is cut off by the strong developing trailing vortex which rolls up in a counterclockwise manner. The trailing edge is also pushed away from the plate by a new growing leading vortex again. This mechanism causes the leading and trailing edge vortices to detach from the upper surface of the plate, repeatedly, with alternative shedding as the example ofRe=200 shown in Figure 24.Re
=
200 streamlines during vortex formation process.
(a)
(b)
(c)
(d)
(e)
(f)Before steady vortex shedding is achieved, the transient stage, the coherent structure can be found as shown in Figures25, 26, and 27 where a strong vorticity gradient is generated on the surface of the plate. Those results indicate that the nonconforming approach also can capture the same transient phenomena as the conforming approach does. Furthermore, if the view is extended more widely as in Figure 28, more obvious coherent structures can be found in Re=200 compared to those in lower Reynolds numbers such as Re=140.Re
=
100, α=30° transient vorticity contours: (a) conforming approach, (b) adaptive nonconforming approach.
(a)
(b)Re
=
140, α=30° transient vorticity contours: (a) conforming approach, (b) adaptive nonconforming approach.
(a)
(b)Re
=
200, α=30° transient vorticity contours: (a) conforming approach, (b) adaptive nonconforming approach.
(a)
(b)α
=
30
° transient vorticity contours: (a) Re=140, (b) Re=200.
(a)
(b)


### 4.3. Comparison

The Strouhal number and surface forces can be used to examine the accuracy of nonconforming computation with conforming one. Based on Strouhal number comparison, as shown in Table1, the vorticity error based method (k=249) only deviates from the conforming meshes result (k=395) by 0.8~1% for various Reynolds numbers. However, the synthesis of velocity and pressure error based method leads to more deviation, ranging from 0.41 to 4.46%. Nevertheless, the estimated error decay rate (σ) based mesh refinement generates the largest deviation among three refinement methods.Table 1
The comparison of results from adaptive nonconforming and conforming methods.
Methods
Conforming
Vorticity (deviation)
V and P (deviation)
σ
<
0.8 (deviation)
St.
Re = 100
0.2434
0.2454 (0.78%)
0.2424 (0.41%)
0.2570 (5.59%)
Re = 140
0.2870
0.2910 (1.08%)
0.2742 (4.96%)
0.2712 (5.51%)
Re = 200
0.3222
0.3194 (0.8%)
0.3125 (1.37%)
0.3086 (4.22%)
Cd
Re = 100
0.2915
0.2907 (0.275%)
0.2814 (3.589%)
0.3144 (7.86%)
Re = 140
0.3389
0.3417 (0.82%)
0.3227 (5.02%)
0.3323 (1.95%)
Re = 200
0.4072
0.4095 (0.56%)
0.3745 (8.73%)
0.3394 (16.65%)
Cl
Re = 100
0.2139
0.2192 (2.417%)
0.2289 (7.01%)
0.2432 (2.48%)
Re = 140
0.3417
0.3432 (0.437%)
0.3268 (4.36%)
0.3144 (7.98%)
Re = 200
0.4842
0.4989 (2.946%)
0.4351 (10.14%)
0.3595 (25.75%)
Notes:
Conforming: Conforming mesh (k=395).
Vorticity: Vorticity error based method (k=249).
V and P: Synthesis of velocity and pressure error based method.
σ
<
0.8: estimated error decay rate (σ) based method which the error decay rate is lower than 0.8.
Deviation: The comparison of the results of the adaptive meshes method and those of the conforming meshes.
St.: Strouhal number.From the drag coefficient history diagrams, the periodical oscillation curve of vorticity error based mesh refinement almost matches with that of conforming meshes in amplitude and period as well as phase angle as shown in Figure29 for Re=140 and Figure 30 for Re=200, respectively. As comparison of the mean drag coefficient, the results from vorticity error based refinement deviate from those from conforming meshes by only 0.27~0.82%, as shown in Table 1. In contrast to the vorticity based method, the mean drag coefficients from the synthesis error of velocity and pressure based method generate larger deviation of 3.59~8.73% from those with conforming approach.Figure 29
Drag coefficient history comparison with conforming and nonconforming approaches,Re=140.Figure 30
Drag coefficient history comparison with conforming and nonconforming approaches,Re=200.Comparing the mean lift coefficient in the nonconforming approaches and conforming method, as shown in Table1, the deviation is about 0.44% for Re=140 and 2.95% for Re=200, when the vorticity error based method is applied. However the deviation is augmented, when the synthesis of velocity and pressure error based method is employed. As seen in the lift coefficient history diagrams shown in Figure 31, for Re=140, the result of nonconforming meshes when using the vorticity error based method matches well with that of conforming meshes. However, for Re=200, there is a slight difference in amplitude between the conforming and nonconforming approaches as shown in Figure 32. Among these three refinement methods, the estimated error decay rate (σ) based mesh refinement generates the greatest deviation with unstable accuracy. However, the vorticity error based method demonstrates its reliability and accuracy with less than 1% deviation from results of conforming method for most of cases.Figure 31
Lift coefficient history comparison with conforming and nonconforming approaches,Re=140.Figure 32
Lift coefficient history comparison with conforming and nonconforming approaches,Re=200.


### 4.4. Computational Efficiency

The globalL2 error is the square root of the sum of each elemental contribution of estimated error which is composed of truncation and quadrature errors as shown below:
(23)εestglobal=(∑k=1K((aNk)2(2N+1)/2+∫N+1∞(ank(n))2(2n+1)/2dn))1/2.
The global estimated L2 error spans the whole computational domain rather than each individual element. Therefore, this global index can indicate the quality and efficiency of mesh adaptation.Within the computational domain, the element which holds the largest error is selected to be refined. Each kind of error estimation is employed to compare the efficiency with the global error convergence rate and the CPU time spent within a specific physical time. As shown in Figure33, for Re=140, the global L2 error is kept in about 10-3 when conforming meshes with 395 elements. For nonconforming meshes, before the mesh refinement process, a large global error results from initial coarse meshes (K=95). The coarse meshes and their solution are used as the input for further mesh evolution. During mesh adaptation process, the global L2 error is declining from initial accuracy, about 10-2, because more elements are deployed as time goes on. Among these three mesh refinement methods, the vorticity based method has the fastest global L2 error decay rate. The sharply reduced error is smaller than that of conforming method after 0.05 physical time unit and it decreases further to be 10-5 before 0.2 physical time unit. The global L2 error of the synthesis of velocity and pressure error based method also decreases to lower than that of conforming method after 0.17 physical time unit but is still higher than that of the vorticity based method. The worst result comes from the estimated error decay rate (σ) method leading an unsatisfactory global L2 error and decay rate. From the point of view of the CPU time spent, as shown in Figure 34, all the adaptive nonconforming methods take much less CPU time than does the conforming approach. For example, to reach 0.15 physical time, the CPU time spent by the vorticity based method is much less than any other method. The synthesis of velocity and pressure method and estimated error decay rate (σ) method both cost almost the same CPU time but still less than that of conforming method. Overall, the CPU time spent by adaptive nonconforming methods is only 20 to 50% of that of conforming method. The global L2 error decay rate and cost of CPU time in Re=200 has the same pattern as those of Re=140. Hence, among all methods, the vorticity based method has the fastest global L2 error decay rate and achieves minimum error with the least computational cost. From the above analyses, the vorticity error based method and synthesis of pressure and velocity error based method are more competitive than the error decay rate method.Figure 33
Global error decay time rate in simulated time scale.Figure 34
Instant CPU time lapse in simulated time scale.


### 4.5. Mesh Refinement on a Curved Element

Other than exploring cases in which mesh refinement is conducted in elements with linear edges, this study also addresses the feasibility of the implementation of mesh refinement on curved elements to explore the potential of application on more complex geometries. The case of flow past a cylinder is used for testing the algorithm of mesh adaptivity on curved elements such as elements around cylinder surfaces. Due to its higher accuracy and efficiency, the estimated error of vorticity is chosen as the mesh refinement indicator. The initial mesh is composed of 95 conforming elements as shown in Figure35. After processes of mesh refinement, there are vigorous mesh refinement activities in the front of the cylinder and on the surface as well as in the wake regions with 253 total element numbers, as shown in Figure 36. The smoothness of vorticity contours as shown in Figure 37 indicates that the accuracy of the results implemented by nonconforming elements is sufficient and this mesh refinement algorithm is applicable for curvy elements. For time average pressure coefficient on the cylinder surface, the result of adaptive nonconforming method is exactly the same as that of the conforming method. Both results agree well with those of Park et al. [38] as shown in Figure 38. For the lift coefficient history, as shown in Figure 39, both of the conforming and adaptive nonconforming methods generate exactly the same results, with both revealing a Strouhal number. The Strouhal number obtained by this adaptive nonconforming mesh method is 0.1816 which deviate by 0.78% from the experimental result, 0.1802, of Williamson [39]. The above results show that the adaptive nonconforming mesh method can also obtain the same accuracy as conforming mesh method does even for curved elements. However, with adaptive nonconforming mesh method, the element number is reduced by 35% compared to the conforming mesh method.Figure 35
Initial mesh is composed by 95 conforming elements.Figure 36
Nonconforming meshes with 253 elements after mesh refinement,Re=140.Figure 37
Vorticity contours, flow past cylinder,Re=140.Figure 38
Surface pressure coefficients of cylinder,Re=140, by conforming (388 elements) and nonconforming (253 elements) meshes.Figure 39
Lift coefficients of cylinder,Re=140, by conforming (388 elements) and nonconforming (253 elements) meshes.


### 4.6. High Reynolds Number Flow

This adaptive nonconforming scheme is also implemented in higher Reynolds number flow to validate its feasibility. In the case of flow past an inclined flat plate in high Reynolds number, the global mesh distribution is similar to that in lower Reynolds number. Figure40 shows that adaptive nonconforming scheme can describe the detailed flow structure well as conforming method does. The flow structures on the surface of plate are slightly different from each other because there is minor time lag of snapshot between two methods. However, the resolution of wake structure obtained by adaptive mesh scheme may be improved by allocating nonconforming meshes in wake region efficiently without increasing computational load. Currently, the total element numbers is limited to less than 260, 66% of conforming meshes, to make sense for choosing adaptive nonconforming scheme. Hence, for considering computational cost, how to improve the precision of mesh distribution in high Reynolds number, especially in high shear stress regions and large vorticity gradient regions, is important for the next step of this study.Re
=
1000, α=30° instantaneous streamlines (a) conforming approach, k=395 (b) adaptive nonconforming approach, k=252.
(a)
(b)


## 4.1. Validation

The case of laminar incompressible flow past a square cavity is commonly used to validate numerical methods. There are discontinuities in velocity on the upper corners, whereu=0 derives from the no-slip condition on the vertical walls and u=1 results due to the imposed driving flow on the upper boundary. Hence, there are singularities on these two corners. The cavity flow regime from Reynolds number 100 to 1000 has been investigated by Ghia et al. [32]. Ghia’s results are used for comparison with our results of conforming meshes and adaptive meshes. The conforming approach with 64 equal-sized elements and an order of shape function of 11 is used for simulation. The nonconforming approach is performed with the initial mesh which has four equal-sized elements with shape function order, N=11.Three mesh refinement indicators which are the synthesis of estimated error of velocity and pressure, estimated error of vorticity, and estimated error decay rate (σ) are employed for the sake of comparison. At each refinement step, the element with the largest error is applied to conduct mesh refinement. After mesh adaptation, the mesh becomes 34 nonconforming elements. The meshes and results are shown in Figure 3 for Re=100 and Figure 4 for Re=1000, respectively. There are two obvious secondary vortices in the bottom corners, especially in the case of Re=1000. The profile of y-directional velocity, v, along the horizontal line through geometric center and the profile of x-directional velocity, u, along the vertical line through the geometric center are used for comparison with the result of multigrid method [32]. The results from all conforming and nonconforming approaches show good agreement with those in [32], as shown in Figures 5 and 6 for Re=100 and Figures 7 and 8 for Re=1000, respectively. From the above cases, we can see that these adaptive nonconforming schemes not only capture the physical phenomena precisely with 53% elements of that of the conforming method but also attain the same accuracy as achieved by the conforming meshes approach and other numerical methods. The global L2 error decays exponentially as the element number is increased which fits the property of h-type mesh refinement as shown in Figure 9.Figure 3
Re
=
100 streamlines, nonconforming approach by vorticity error base.Figure 4
Re
=
1000 streamlines, nonconforming approach by vorticity error base.Figure 5
Re
=
100, v-profile along central line of x-span in square cavity.Figure 6
Re
=
100, u-profile along central line of x-span in square cavity.Figure 7
Re
=
1000, v-profile along central line of x-span in square cavity.Figure 8
Re
=
1000, u-profile along central line of x-span in square cavity.Figure 9
Global error decay in adaptive mesh refinement, cavity flow,Re=1000.


## 4.2. Flow Past a Flat Plate

Flow past a thin airfoil is an interesting topic for science and engineering because it abounds with rich physical phenomena such as shear layers, leading and trailing edge vorticity, laminar separation, and vortex shedding. It is also important for microaerial vehicles (MAVs) design.Lam and Leung [33] investigated experimentally the shedding of vortices from an inclined flat plate. The results showed that the trailing edge vortex in the separated wake possesses more intense vorticity levels but a smaller vortex size than the leading edge vortex. Breuer and Jovičić [34] conducted a large eddy simulation (LES) for the flow past an inclined flat plate in Re=20000, at an angle of attack of 18°. The findings revealed Kelvin-Helmholtz instability in the free shear layer behind the leading edge. A highly asymmetric wake with vortices of unequal strength was observed in the wake region.However, there have been few studies of flow past a flat plate in hundreds of Reynolds number. Taira et al. [35] studied the flow past an inclined thin flat plate with immersed boundary method. The numerical results of flow Re=100 past a three-dimensional flat plate with an aspect ratio=2 agree with experimental lift and drag coefficients in various angles of attack. The effects of aspect ratio and planform geometry on the leading-edge vortex were studied with Re=300 flow. The result showed that the elliptic and semicircular plates generate a wake with less distinct transition between leading edge and tip vortices than that of rectangular profile. Hence, the flow control for flat plate in low Reynolds numbers has aroused the attentions of several researchers. For instance, Taira et al. [36] simulated a synthetic jet type actuator to modify the dynamics of wake vortices and the corresponding forces exerted on wings. Brunton and Rowley [37] further used the immersed boundary method, a direct numerical method, and proper orthogonal decomposition/Galerkin projection to simulate stationary flat plate in Re=100, 300 flow with angles of attack, α=30° and α=45°, and the full dynamic wake structures of a pitching or plunging flat plate.In our study, for constructing quadrilateral spectral elements on the surface of a flat plate, an extremely thin two-dimensional inclined flat plate with maximum thickness to chord ratio, 1 × 10−5, is used to approach a theoretical straight-line-like flat plate profile. Hence, there are very sharp leading edge and trailing edge which are prone to triggering some critical physical phenomena such as leading edge vortex and trailing edge vortex. A rectangular global computational domain is selected as in Figure 10. The chord of the plate is 1 unit (D) with an angle of attack 30°. The domain is extended from the center of the plate to upstream boundary with distance of 20 D and to downstream boundary with distance of 30 D. The distances from upper and lower boundaries to the plate center are both 20 D. The uniform flow boundary conditions (u=U0, v=0) are employed in the inflow and upper and lower boundaries so as to set different flow conditions, Re=100, 140, and 200. The outflow boundary condition is employed in the most downstream boundary. There are 395 macrospectral elements (K), employed for conforming approach. The order (N) of the basis function, Gauss-Lobatto-Legendre polynomial, is 11.Figure 10
Flat plateα=30°, K=395, conforming meshes.In the nonconforming approach, the initial coarse mesh with 95 elements used for simulation leads to unsatisfactory result. This result is used as the initial input of mesh adaption process. Figure11 shows the results of the mesh refinement after carrying out mesh adaptation with the indicator of vorticity estimated error 50 times. The mesh refinement regions are concentrated in the area around the flat plate, wake region, and upwind region. The refinement regions do not vary much for these three Reynolds numbers. The results of adaptive nonconforming meshes based on vorticity error are shown in Figures 13, 15, and 17. These results not only show good agreement with those counterparts done by conforming meshes as shown in Figures 12, 14, and 16 but also can define clearer vorticity contours.Figure 11
Flat plateα=30°, K=245, nonconforming meshes refinement started from initial meshes (K=95) by vorticity error indicator.Figure 12
Conforming approach with 395 elements, vorticity,Re=100, α=30°.Figure 13
Nonconforming approach by vorticity error indicator,Re=100, α=30°, vorticity with 245 elements.Figure 14
Conforming approach with 395 elements, vorticity,Re=140, α=30°.Figure 15
Nonconforming approach by vorticity error indicator,Re=140, α=30°, vorticity with 247 elements.Figure 16
Conforming approach with 395 elements, vorticity,Re=200, α=30°.Figure 17
Nonconforming approach by vorticity error indicator,Re=200, α=30°, vorticity with 249 elements.Among those mesh refinement regions, more intensive mesh refinements are executed around the flat plate, especially in the leading edge and the trailing edge where several tiny elements are employed, automatically, as shown in Figures19, 21, and 23. The result shows that vorticity, the first derivative property, can be solved smoothly by the adaptive nonconforming method. This mesh distribution indicates that the vorticity error occurs in the stronger shear layer area where the vorticity gradient is also larger. In addition, since the sharp leading edge causes more difficulties in computation, it attracts more mesh refinement activities. The results from adaptive nonconforming meshes show good agreements with those from conforming meshes as shown in Figures 18, 20, and 22. However, the element number, 247, required by adaptive nonconforming approach is much less than the 395 elements used in the conforming approach. Hence, only 62% of the element numbers are needed when adaptive mesh refinement is employed. This reveals that adaptive mesh approach can compute more efficiently than the conforming approach does, if meshes are employed in the right locations by a proper mesh refinement indicator.Figure 18
Conforming approach with 395 elements, closer look vorticity contours,Re=100, α=30°.Figure 19
Nonconforming approach with 247 elements, closer look vorticity contours,Re=100, α=30°.Figure 20
Conforming approach with 395 elements, closer look vorticity contours,Re=140, α=30°.Figure 21
Nonconforming approach with 247 elements, closer look vorticity contours,Re=140, α=30°.Figure 22
Conforming approach with 395 elements, closer look vorticity contours,Re=200, α=30°.Figure 23
Nonconforming approach with 249 elements, closer look vorticity contours,Re=200, α=30°.The clockwise leading edge vortex first forms and grows. Later, it is cut off by the strong developing trailing vortex which rolls up in a counterclockwise manner. The trailing edge is also pushed away from the plate by a new growing leading vortex again. This mechanism causes the leading and trailing edge vortices to detach from the upper surface of the plate, repeatedly, with alternative shedding as the example ofRe=200 shown in Figure 24.Re
=
200 streamlines during vortex formation process.
(a)
(b)
(c)
(d)
(e)
(f)Before steady vortex shedding is achieved, the transient stage, the coherent structure can be found as shown in Figures25, 26, and 27 where a strong vorticity gradient is generated on the surface of the plate. Those results indicate that the nonconforming approach also can capture the same transient phenomena as the conforming approach does. Furthermore, if the view is extended more widely as in Figure 28, more obvious coherent structures can be found in Re=200 compared to those in lower Reynolds numbers such as Re=140.Re
=
100, α=30° transient vorticity contours: (a) conforming approach, (b) adaptive nonconforming approach.
(a)
(b)Re
=
140, α=30° transient vorticity contours: (a) conforming approach, (b) adaptive nonconforming approach.
(a)
(b)Re
=
200, α=30° transient vorticity contours: (a) conforming approach, (b) adaptive nonconforming approach.
(a)
(b)α
=
30
° transient vorticity contours: (a) Re=140, (b) Re=200.
(a)
(b)


## 4.3. Comparison

The Strouhal number and surface forces can be used to examine the accuracy of nonconforming computation with conforming one. Based on Strouhal number comparison, as shown in Table1, the vorticity error based method (k=249) only deviates from the conforming meshes result (k=395) by 0.8~1% for various Reynolds numbers. However, the synthesis of velocity and pressure error based method leads to more deviation, ranging from 0.41 to 4.46%. Nevertheless, the estimated error decay rate (σ) based mesh refinement generates the largest deviation among three refinement methods.Table 1
The comparison of results from adaptive nonconforming and conforming methods.
Methods
Conforming
Vorticity (deviation)
V and P (deviation)
σ
<
0.8 (deviation)
St.
Re = 100
0.2434
0.2454 (0.78%)
0.2424 (0.41%)
0.2570 (5.59%)
Re = 140
0.2870
0.2910 (1.08%)
0.2742 (4.96%)
0.2712 (5.51%)
Re = 200
0.3222
0.3194 (0.8%)
0.3125 (1.37%)
0.3086 (4.22%)
Cd
Re = 100
0.2915
0.2907 (0.275%)
0.2814 (3.589%)
0.3144 (7.86%)
Re = 140
0.3389
0.3417 (0.82%)
0.3227 (5.02%)
0.3323 (1.95%)
Re = 200
0.4072
0.4095 (0.56%)
0.3745 (8.73%)
0.3394 (16.65%)
Cl
Re = 100
0.2139
0.2192 (2.417%)
0.2289 (7.01%)
0.2432 (2.48%)
Re = 140
0.3417
0.3432 (0.437%)
0.3268 (4.36%)
0.3144 (7.98%)
Re = 200
0.4842
0.4989 (2.946%)
0.4351 (10.14%)
0.3595 (25.75%)
Notes:
Conforming: Conforming mesh (k=395).
Vorticity: Vorticity error based method (k=249).
V and P: Synthesis of velocity and pressure error based method.
σ
<
0.8: estimated error decay rate (σ) based method which the error decay rate is lower than 0.8.
Deviation: The comparison of the results of the adaptive meshes method and those of the conforming meshes.
St.: Strouhal number.From the drag coefficient history diagrams, the periodical oscillation curve of vorticity error based mesh refinement almost matches with that of conforming meshes in amplitude and period as well as phase angle as shown in Figure29 for Re=140 and Figure 30 for Re=200, respectively. As comparison of the mean drag coefficient, the results from vorticity error based refinement deviate from those from conforming meshes by only 0.27~0.82%, as shown in Table 1. In contrast to the vorticity based method, the mean drag coefficients from the synthesis error of velocity and pressure based method generate larger deviation of 3.59~8.73% from those with conforming approach.Figure 29
Drag coefficient history comparison with conforming and nonconforming approaches,Re=140.Figure 30
Drag coefficient history comparison with conforming and nonconforming approaches,Re=200.Comparing the mean lift coefficient in the nonconforming approaches and conforming method, as shown in Table1, the deviation is about 0.44% for Re=140 and 2.95% for Re=200, when the vorticity error based method is applied. However the deviation is augmented, when the synthesis of velocity and pressure error based method is employed. As seen in the lift coefficient history diagrams shown in Figure 31, for Re=140, the result of nonconforming meshes when using the vorticity error based method matches well with that of conforming meshes. However, for Re=200, there is a slight difference in amplitude between the conforming and nonconforming approaches as shown in Figure 32. Among these three refinement methods, the estimated error decay rate (σ) based mesh refinement generates the greatest deviation with unstable accuracy. However, the vorticity error based method demonstrates its reliability and accuracy with less than 1% deviation from results of conforming method for most of cases.Figure 31
Lift coefficient history comparison with conforming and nonconforming approaches,Re=140.Figure 32
Lift coefficient history comparison with conforming and nonconforming approaches,Re=200.


## 4.4. Computational Efficiency

The globalL2 error is the square root of the sum of each elemental contribution of estimated error which is composed of truncation and quadrature errors as shown below:
(23)εestglobal=(∑k=1K((aNk)2(2N+1)/2+∫N+1∞(ank(n))2(2n+1)/2dn))1/2.
The global estimated L2 error spans the whole computational domain rather than each individual element. Therefore, this global index can indicate the quality and efficiency of mesh adaptation.Within the computational domain, the element which holds the largest error is selected to be refined. Each kind of error estimation is employed to compare the efficiency with the global error convergence rate and the CPU time spent within a specific physical time. As shown in Figure33, for Re=140, the global L2 error is kept in about 10-3 when conforming meshes with 395 elements. For nonconforming meshes, before the mesh refinement process, a large global error results from initial coarse meshes (K=95). The coarse meshes and their solution are used as the input for further mesh evolution. During mesh adaptation process, the global L2 error is declining from initial accuracy, about 10-2, because more elements are deployed as time goes on. Among these three mesh refinement methods, the vorticity based method has the fastest global L2 error decay rate. The sharply reduced error is smaller than that of conforming method after 0.05 physical time unit and it decreases further to be 10-5 before 0.2 physical time unit. The global L2 error of the synthesis of velocity and pressure error based method also decreases to lower than that of conforming method after 0.17 physical time unit but is still higher than that of the vorticity based method. The worst result comes from the estimated error decay rate (σ) method leading an unsatisfactory global L2 error and decay rate. From the point of view of the CPU time spent, as shown in Figure 34, all the adaptive nonconforming methods take much less CPU time than does the conforming approach. For example, to reach 0.15 physical time, the CPU time spent by the vorticity based method is much less than any other method. The synthesis of velocity and pressure method and estimated error decay rate (σ) method both cost almost the same CPU time but still less than that of conforming method. Overall, the CPU time spent by adaptive nonconforming methods is only 20 to 50% of that of conforming method. The global L2 error decay rate and cost of CPU time in Re=200 has the same pattern as those of Re=140. Hence, among all methods, the vorticity based method has the fastest global L2 error decay rate and achieves minimum error with the least computational cost. From the above analyses, the vorticity error based method and synthesis of pressure and velocity error based method are more competitive than the error decay rate method.Figure 33
Global error decay time rate in simulated time scale.Figure 34
Instant CPU time lapse in simulated time scale.


## 4.5. Mesh Refinement on a Curved Element

Other than exploring cases in which mesh refinement is conducted in elements with linear edges, this study also addresses the feasibility of the implementation of mesh refinement on curved elements to explore the potential of application on more complex geometries. The case of flow past a cylinder is used for testing the algorithm of mesh adaptivity on curved elements such as elements around cylinder surfaces. Due to its higher accuracy and efficiency, the estimated error of vorticity is chosen as the mesh refinement indicator. The initial mesh is composed of 95 conforming elements as shown in Figure35. After processes of mesh refinement, there are vigorous mesh refinement activities in the front of the cylinder and on the surface as well as in the wake regions with 253 total element numbers, as shown in Figure 36. The smoothness of vorticity contours as shown in Figure 37 indicates that the accuracy of the results implemented by nonconforming elements is sufficient and this mesh refinement algorithm is applicable for curvy elements. For time average pressure coefficient on the cylinder surface, the result of adaptive nonconforming method is exactly the same as that of the conforming method. Both results agree well with those of Park et al. [38] as shown in Figure 38. For the lift coefficient history, as shown in Figure 39, both of the conforming and adaptive nonconforming methods generate exactly the same results, with both revealing a Strouhal number. The Strouhal number obtained by this adaptive nonconforming mesh method is 0.1816 which deviate by 0.78% from the experimental result, 0.1802, of Williamson [39]. The above results show that the adaptive nonconforming mesh method can also obtain the same accuracy as conforming mesh method does even for curved elements. However, with adaptive nonconforming mesh method, the element number is reduced by 35% compared to the conforming mesh method.Figure 35
Initial mesh is composed by 95 conforming elements.Figure 36
Nonconforming meshes with 253 elements after mesh refinement,Re=140.Figure 37
Vorticity contours, flow past cylinder,Re=140.Figure 38
Surface pressure coefficients of cylinder,Re=140, by conforming (388 elements) and nonconforming (253 elements) meshes.Figure 39
Lift coefficients of cylinder,Re=140, by conforming (388 elements) and nonconforming (253 elements) meshes.


## 4.6. High Reynolds Number Flow

This adaptive nonconforming scheme is also implemented in higher Reynolds number flow to validate its feasibility. In the case of flow past an inclined flat plate in high Reynolds number, the global mesh distribution is similar to that in lower Reynolds number. Figure40 shows that adaptive nonconforming scheme can describe the detailed flow structure well as conforming method does. The flow structures on the surface of plate are slightly different from each other because there is minor time lag of snapshot between two methods. However, the resolution of wake structure obtained by adaptive mesh scheme may be improved by allocating nonconforming meshes in wake region efficiently without increasing computational load. Currently, the total element numbers is limited to less than 260, 66% of conforming meshes, to make sense for choosing adaptive nonconforming scheme. Hence, for considering computational cost, how to improve the precision of mesh distribution in high Reynolds number, especially in high shear stress regions and large vorticity gradient regions, is important for the next step of this study.Re
=
1000, α=30° instantaneous streamlines (a) conforming approach, k=395 (b) adaptive nonconforming approach, k=252.
(a)
(b)


## 5. Conclusions

Three kinds of properties, the synthesis estimated error of velocity and pressure, the estimated error coefficient decay rate (σ), and vorticity estimated error are chosen as the indicators for mesh refinement. Among those, the vorticity error based method performs with the fastest convergence and highest accuracy with least computational cost. The convergence and accuracy of the estimated error decay rate (σ) method are unsatisfactory. From the comparison of Strouhal numbers, mean drag, and lift coefficients from different methods, the results of the vorticity error based adaptive nonconforming method are the most consistent with those from conforming method. The drag and lift history diagram of the vorticity error based adaptive nonconforming method match well with ones using the conforming method. This study shows that the adaptive nonconforming method guided by vorticity error is a useful technique which performs the same or even better than the conforming method with respect to accuracy and computational cost. The results show that this adaptive mesh method can be applied on the mesh with linear or curved boundaries. Therefore, this mesh adaptation scheme can be implemented in the simulations of flow with more complex geometric boundaries in the future. In addition, how to allocate nonconforming meshes precisely on the most needed region in high Reynolds number external flow to capture detailed physical phenomena without increasing computational load is important for further investigation.



---

*Source: 959615-2014-09-25.xml*