c
TU Delft, The Netherlands, 2006
THE SAMPLING BASED DYNAMIC PROCEDURE FOR
NUMERICAL DISCRETIZATION ENHANCEMENT
Dieter Fauconnier∗, Chris De Langhe and Erik Dick
Ghent University, Faculty of Engineering, St.-Pietersnieuwstraat 41, B-9000 Gent, Belgium
e-mail∗: dieter.fauconnier@ugent.be
web page: http://www.floheacom.ugent.be/
Key words: Dynamic Procedure, High-order accuracy, Richardson Extrapolation Abstract. In CFD computations, discretization errors or truncation errors should be small providing an acceptable level of accuracy . In this paper an alternative use is made of the recently proposed LES formalism based on sampling operators. It is shown that the sampling-based dynamic procedure, in combination with an appropriate truncation error model, can be used as a technique to increase the numerical accuracy of a discretization. The technique is resemblant to the well-known Richardson extrapolation. The procedure is tested on a 1D convection-diffusion equation and a 2D lid-driven cavity at Re = 400, using a finite difference method. Very promising results are obtained.
1 INTRODUCTION
All numerical simulations (laminar, DNS, LES,...) are liable to numerical discretization errors, due to the finite representation of the derivative operators on the computational grid. The reliability and accuracy of a simulation depend partially on the ability to control these discretization errors. It is known that these errors can be quite large for high wavenumbers appearing e.g. in steep gradients or small vortices. Especially in LES, the smallest resolved scales of the spectrum contain a significant amount of energy [1, 2, 3]. Therefore, these errors should be small enough providing an acceptable level of accuracy. However, the complexity of implementing high-order methods, to improve the numerical accuracy, or the computational cost of fine meshes are prohibitive for CFD simulations in complex geometries.
was proposed that, by relying on a so-called generalized dynamic procedure, succeeded in accounting for the subgrid scales. In this paper, we investigate the ability of this sampling based dynamic procedure, in combination with an appropriate model for the truncation error to obtain higher-order numerical accuracy. Two such possible models are presented. Further, we show that Richardson extrapolation is a special formulation of this procedure. First we introduce the sampling formalism in a finite difference context. Then, the generalized dynamic procedure for truncation error modelling is explained and analyzed, and its relation to Richardson extrapolation is discussed. We propose two different models for the truncation error of the Navier-Stokes equations, including a Smagorinsky-like model. Finally, in order to evaluate the numerical qualities of the proposed method, without any turbulence modeling ambiguities, we test the concept on a 1D convection-diffusion equation and a 2D laminar lid-driven cavity at Re = 400.
2 THE SAMPLING FORMALISM
Consider the continuity equation and the Navier-Stokes equations1, for the vector field
~
u(~x, t) and the pressure field p (~x, t) of an incompressible fluid (ρ = 1) in IRn, n ∈ {1, 2, 3}
∂ui ∂xi = 0 (1) ∂ui ∂t + uj ∂ui ∂xj = −∂p ∂xi + ν∂ 2u i ∂x2 j (2) Projecting the Navier-Stokes equations from a continuum physical domain Ω ⊂ IRn to a corresponding discrete physical domain Ω∆1
= {~x1, . . . , ~xN}, ~xl ∈ Ω, with N the number
of grid points, representing the grid with grid spacing ∆1, requires the definition of an
appropriate projection operator.
Hence, we define the sampling operator S∆1
, which operates between Ω and Ω∆1
. This sampling operator S∆1 is idempotent, and commutative with the product of the non-linear
terms. S∆1 ◦ S∆1 = S∆1 (3) S∆1 ◦ [φ · ψ] = S∆1 ◦ [φ] · S∆1 ◦ [ψ] (4) However, S∆1
does not commute with spatial derivatives as they cannot be evaluated at an infinitesimal interval ǫ.
dφ
dx 6= lim∆→ǫ>0
∆φ ∆x We use the notation S∆1 ◦ u
i = ui and S∆1 ◦ ∂ = δ.
Applying now S∆1 to the continuity equation (1) and the Navier-Stokes equations (2) 1
gives δui δxi = Π∆1 (5) ∂ui ∂t + uj δui δxj = −δp δxi + νδ 2u i δx2 j + Σ∆1 i (6)
leading to truncation errors in both continuity and momentum equations due to the non-commutativity of the operator S∆1
with the spatial derivatives. These truncation terms have the basic form
Π∆1 = δui δxi − ∂ui ∂xi (7) Σ∆1 i = uj δui δxj − ∂ui ∂xj + δp δxi − ∂p ∂xi − ν δ 2u i δx2 j − ∂2ui ∂x2 j ! (8) Although both errors vanish in the limit N → ∞, they are generally not negligible for finite grids. Given a discretization scheme, the exact form for the truncation errors (7) and (8) can be obtained from Taylor series expansion, provided that the field is sufficiently smooth on the grid. We choose the discretization scheme a priori to be 2nd
order central difference approximation for both 1st and 2nd order partial derivatives. The
grid Ω∆1 = {~x
1, . . . , ~xN}, ~xl∈ Ω has a uniform spacing ∆xi in each spatial direction. The
finite difference approximations of the derivatives of a scalar ui with respect to xj in a
node xj = xlj are δui δxj = ui x l+1 j − ui xl−1j 2∆xj (9) δ2ui δx2 j = ui x l+1 j − 2ui xlj + ui xl−1j ∆x2 j (10) Then the exact Taylor series expansion of the continuous derivative can be written as2
∂ui ∂xj = δui δxj − 1 6∆x 2 j ∂3u i ∂x3 j − 1 120∆x 4 j ∂5u i ∂x5 j − O ∆x6 (11) ∂2u i ∂x2 j = δ 2u i δx2 j − 1 12∆x 2 j ∂4u i ∂x4 j − 1 360∆x 4 j ∂6u i ∂x6 j − O ∆x6 (12)
Substituting the previous expressions in the continuity equation (1) and the Navier-Stokes equations (2), which is equivalent to applying the sample operator S∆1, leads to the exact
2
expressions for the truncation terms Π∆1(7)and Σ∆1 i (8) Π∆1 = 1 6∆x 2 i ∂3u i ∂x3i + 1 120∆x 4 i ∂5u i ∂x5i + O ∆x 6 (13) Σ∆1 i = uj " 1 6∆x 2 j ∂3u i ∂x3 j + 1 120∆x 4 j ∂5u i ∂x5 j + O ∆x6 # + " 1 6∆x 2 i ∂3p ∂x3i + 1 120∆x 4 i ∂5p ∂x5i + O ∆x 6 # −ν " 1 12∆x 2 j ∂4u i ∂x4j + 1 360∆x 4 j ∂6u i ∂x6j + O ∆x6 # (14) Our objective is to obtain a model that increases the accuracy of the second order ap-proximations, by using information from two different grid resolutions, in the philosophy of the generalized sampling based dynamic procedure.
3 THE GENERALIZED DYNAMIC PROCEDURE FOR TRUNCATION
ERROR MODELLING
3.1 Concept
The original dynamic procedure, based on the Germano [7] identity can be extended to a more general approach in the sampling context. Note that Jeanmart and Winck-elmans [8], already suggested the use of a sampling operator in the dynamic procedure. The same terminology as in Section (2) is used. Projection is done of the Navier-Stokes equations from a continuum physical domain Ω ⊂ IRnto a corresponding discrete physical
domain Ω∆1 = {~x
1, . . . , ~xN1}, ~xl ∈ Ω, and to Ω
∆2 = {~x
1, . . . , ~xN2}, ~xl ∈ Ω with N2 < N1
(in practice, Ω∆2 ⊂ Ω∆1
). This corresponds with the sampling operators S∆1
and S∆2
, projecting respectively Ω → Ω∆1 and Ω → Ω∆2. S∆2 also projects Ω∆1 → Ω∆2, since
S∆2
◦ S∆1
= S∆2
(15) We keep the notation S∆1 ◦ u
i = ui and introduce S∆2 ◦ ui = eui = eui. We also keep the
same notation for the discrete derivative-operator S∆2 ◦ ∂ = δ. Applying the operator
S∆1 on the continuity equation (1) and Navier-Stokes equations (2) leads to
0 = C∆1 (ui) + Π∆1 = − δui δxi + Π∆1 (16) ∂ui ∂t = N ∆1 i (ui) + Σ∆i 1 = −uj δui δxj − δp δxi + νδ 2u i δx2j + Σ ∆1 i (17) Π∆1and Σ∆1
i are truncation errors that are known theoretically, but cannot be evaluated
numerically. C∆2 and N∆2
Subsequently, applying S∆2 to the continuous set (1) and (2) gives 0 = C∆2 (eui) + Π∆2 = − δuei δxi + Π∆2 (18) ∂uei ∂t = N ∆2 i (eui) + Σ∆i 2 = −euj δuei δxj − δpe δxi + νδ 2eu i δx2 j + Σ∆2 i (19)
Ideally, the latter set should also be obtained by applying the sampling operator S∆2
to the first set of equations (16-17) giving
0 = S∆2 ◦ C∆1 (ui) + S∆2 ◦ Π∆1 (20) ∂uei ∂t = S ∆2 ◦ N∆1 i (ui) + S∆2 ◦ Σ∆i 1 (21)
Consistency between formulations (18-19) and (20-21) imposes the following relations: S∆2 ◦ C∆1 (ui) − C∆2(eui) = Π∆2 − S∆2 ◦ Π∆1 (22) S∆2 ◦ N∆1 i (ui) − Ni∆2(eui) = Σ∆i 2 − S∆ 2 ◦ Σ∆1 i (23)
These explicitly express the commutation errors made by the projection Ω∆1 → Ω∆2.
The left-hand sides of (22) and (23) are scalar level and vector level equivalents of the Germano identity, respectively. They can be determined in terms of the resolved velocity ui since S∆2 ◦ ui = eui = eui and play the role of the Leonard scalar (22) or the Leonard
vector in (23). Suppose models are adopted for both truncation errors Π∆1
and Σ∆1
i which
have the basic forms
Π∆1 = Cπmπ,∆1 (24) Σ∆1 i = C σ im σ,∆1 i (25)
and analogously for the test-level ∆2. Using a similar terminology as in the classic dynamic
procedure, the expressions (22) and (23) can be written as
Lπ = CπMπ (26) Lσ i = C σ iM σ i (27)
The Leonard terms are thus resemblant to the expressions (7-8). Explicitely written, they are given as Lπ = gδui δxi − δeui δxi (32) Lσ i = euj gδui δxj − δeui δxj ! + gδp δxi − δep δxi ! − ν δg 2u i δx2 j − δ 2eu i δx2 j ! (33) Similar expressions are obtained for the model terms. In the approach of Knaepen et al. [5], a single field Cσ was used to model the error of the separate momentum equations leading
to a least square optimization procedure as a compromise between the 3 independent conditions (27). Here we propose separate fields Cσ
i for the separate equations, from
which the optimal parameter can be determined for every independent condition (27). It could be necessary to filter the quantities before calculating the scalar fields Cπ and
Ciσ, to remove undesirable high frequency pollution or singularities, leading to instability
of the algorithm. We suggest to use a least square method leading to3
Cπ = hL πMπi hMπMπi (34) Ciσ = hL σ iMσii hMσ iMσii (35) where a smoothing filter is applied, such as a local moving average or a global average. For the latter averaging, constant fields are obtained. Cπ and Cσ
i can be positive or negative,
dependent on the adopted discretization scheme. Clipping may be necessary, if excessive values jeopardize the algoritm’s stability. Finally, Cπ and Cσ
i, given on the coarse grid,
are interpolated to the fine grid using a piecewise Cubic Hermite interpolation. A fully embedded test grid is applied.
3.2 Analysis
We use the notation δnu
δxn
∆
= δnu
δxn (x, ∆). Consider the exact Taylor series expansion of
the n-th order derivative, for a k-th order central discretization scheme ∂nu ∂xn (x) = δnu δxn ∆ + ck∆k ∂k+nu ∂xk+n + O ∆ k+2 (36) ∂nu ∂xn (x) = δnue δxn 2∆ + ck(2∆)k ∂k+nu ∂xk+n + O ∆ k+2 (37)
We assume now, that the leading order truncation term is an adequate model for the com-plete truncation error. Discretization of this term and applying the generalized dynamic
3
procedure leads to a coefficient ck ck = δnue δxn 2∆ − δ nu δxn ∆ ∆k δk+nu δxk+n ∆ − (2∆)k δ k+nue δxk+n 2∆ (38)
Substitution of which in (36) results finally in
∂nu ∂xn(x) ≈ 2kδ nu δxn ∆ δk+nue δxk+n 2∆ − δ nue δxn 2∆ δk+nu δxk+n ∆ 2kδk+neu δxk+n 2∆ − δ k+nu δxk+n ∆ (39)
This expression is closely related to a generalized form of Richardson extrapolation. Under the assumption δk+nue δxk+n 2∆ ≡ δ k+nu δxk+n ∆ ≡ ∂ k+nu ∂xk+n (40)
the generalized Richardson extrapolation formula is obtained ∂nu ∂xn(x) ≈ 2kδ nu δxn(x, ∆) − δn e u δxn (x, 2∆) 2k− 1 + O ∆ k+2 (41)
which is fully equivalent with a (k + 2)th order accurate central scheme for the nth
deriv-ative.
For further analysis of expression (39), and the implications of approximation (40), a Fourier analysis is performed for the 2nd order accurate gradient. We adopt the Taylor
series expansion ∂u ∂x(x) ≈ δu δx ∆ + c∆2 δ 3u δx3 ∆ (42) ∂u ∂x(x) ≈ δeu δx 2∆ + c (2∆)2 f δ 3u δx3 2∆ + (1 − f ) δ 3u δx3 ∆! (43) resulting with the dynamic procedure in
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 f=1/5 f=1 f=1/2 f=1/3 f=1/4 f=1/10 κ/κ max κ′ / κmax
Figure 1: Modified wavenumbers for ∂u∂x: (△), spectral; (. . .), 2 nd
order central; (- . -), 4th
order central; (- - -), 6th order central; (—), Dynamic Procedure.
The blending factor f is introduced in order to switch between the dynamic procedure (39) (f = 1) and Richardson extrapolation (41) (f = 0), and to investigate intermediate behaviour.
In Fourier space, the gradient can be written as F ∂u ∂x = iκ′F (u) (45)
with κ′ the modified wavenumber. Modified wavenumbers are given in figure 1 for a
2nd, 4th (equals Richardson Extrapolation f = 0) and 6thorder central scheme, and also
for the expression (44). The latter is given at values of f = 1,12,13,14,15,101. The full dynamic procedure f = 1 leads to singularities for certain wavenumbers, and is for this reason useless. We noticed from a semilog error plot that f = 15 is the value for which the dynamic procedure obtains best results. For the lower wavenumbers the accuracy is slightly better than 6th order, for the higher wavenumbers it is slightly less.
4 MODELLING OF THE TRUNCATION ERRORS
as the basic modeling ingredients, and rewrite (13) and (14) as Π∆1 = 1 6+ O ∆ 2 ∆x2i ∂3u i ∂x3 i (46) Σ∆1 i = uj 1 6 + O ∆ 2 ∆x2j ∂3u i ∂x3j + 1 6 + O ∆ 2 ∆x2i ∂3p ∂x3i −ν 2 1 6+ O ∆ 2 ∆x2j∂ 4u i ∂x4 j (47) We now make the assumption that to a certain degree of accuracy we can merge the different series between the brackets in (47) into the vector field Cσ
i . Similarly denoting
the terms between brackets in (46) by the scalar field Cπ, we obtain as modeling basis
Π∆1 = Cπ∆x2i δ3u i δx3i (48) Σ∆1 i = C σ i uj∆x2j δ3ui δx3 j + ∆x2i δ3p δx3 i − ν 2∆x 2 j δ4ui δx4 j (49) For a 2D flow, this results in
Π∆1 = Cπ ∆x2δ 3u δx3 + ∆y 2δ3v δy3 (50) Σ∆1 u = C σ u ∆x2 uδ 3u δx3 + δ3p δx3 − ν 2 δ4u δx4 +∆y2 vδ 3u δy3 − ν 2 δ4u δy4 (51) Σ∆1 v = Cvσ ∆x2 uδ 3v δx3 − ν 2 δ4v δx4 +∆y2 vδ 3v δy3 + δ3p δy3 − ν 2 δ4v δy4 (52) This model is closely related to the exact expression of the discretization error. The major drawback is the requirement of a broader 5-point stencil to evaluate the 3rd and
4th derivatives. This can be very unpleasant near walls, where excentric 6-point stencils
(also second order accurate) have to be constructed to maintain the overall accuracy. It would be more convenient if the higher order derivatives could be reduced to maximum 2nd order, the highest appearing order in the physical Navier-Stokes equations.
The general idea behind the reduction of the high order derivatives for modeling purposes, is based on the observation that in Fourier space, higher order derivatives and certain products of lower order derivatives have identical modified wavenumbers. In physical space, such an approximation preserves the same character although it is not identical. We propose therefore the approximation
In order to obtain the same modified wavenumber in (53), the condition n = p · q + r should be satisfied. If n = 3, and if we choose p = 1, q = 1 and r = 2, a Smagorinsky-like approximation appears. The convective terms in (47) can be transformed, using (53) into
uj δ3ui δx3j ≈ uj ui δui δxj δ2ui δx2j ≈ δui δxj δ2ui δx2j (54) For the truncation error of the continuity equation (46), the pressure terms and the viscous terms in the truncation error of the Navier-Stokes equations (47), the foregoing hypothesis (53) does not lead directly to a suitable Smagorinksy-like model. Therefore, we decide, in a first attempt to model the truncation terms in a Smagorinsky-like fashion, to leave the truncation model for the continuity equation in its form (48), and to neglect the pressure terms and viscous terms in (47). This finally results in a Smagorinsky-like model for the momentum equations Σ∆1 i = C σ∗ i ∆x2j δui δxj δ2u i δx2 j (55) where an additional approximation is made by giving all terms in the truncation error model the same coefficient. For a 2D flow this results in
Π∆1 = Cπ ∆x2δ 3u δx3 + ∆y 2δ3v δy3 (56) Σ∆1 u = C σ∗ u ∆x2 δu δx δ2u δx2 + ∆y 2 δu δy δ2u δy2 (57) Σ∆1 v = C σ∗ v ∆x2 δv δx δ2v δx2 + ∆y 2 δv δy δ2v δy2 (58) 5 RESULTS
In order to evaluate the numerical qualities of the procedure, in the absence of ambi-guities caused by turbulence modeling aspects, we choose to test it on a laminar flow. We first analyse the method on a 1D convection-diffusion equation. Then we study the more complex flow of a 2D driven cavity at Reynolds number Re = 400
5.1 1D convection-diffusion equation
The 1D convection-diffusion equation is a very interesting canonical problem, which represents boundary layer behavior. Within the domain Ω = [0, L] ⊂ IR, this equation supplemented with a set of Dirichlet boundary conditions, is given by
Code Σ∆1 i κt n cd.a Cσ∆x2 3 cd.b Cσ∆x2 κ 2 cd.c 1 6∆x 2 ∂3 u ∂x3 − κ 12∆x 2 ∂4 u ∂x4 Cσ∆x 2 cd.d Cσ∆x2 c δu δx 2
Table 1: Summary of the different simulations for an 2nd
order central scheme
with c the convection speed and ν the viscosity. If we define the P´eclet number P e = Lκ , with κ = νc , the analytical solution reads
u(x) = u0+ (uL− u0)
1 − exκ
1 − eP e; (60)
Discretization on a uniform grid gives δu δx = κ δ2u δx2 + Σ ∆1 (61) u(x1) = u0 u(xN) = uL We propose for Σ∆1
the basic form Σ∆1
= κt
δnu
δxn (62)
An overview of the different models for an overall 2nd order central discretization are given
in Table 1. (Tests with a 1st order upwind and a 3rd order Kappa-Van Leer scheme for
the convective terms were also performed. These results are not included here.) Note that for Cσ = 1
12 an exact 4
th order accuracy is reached for model cd.a. Because of the
relation ∂u ∂x = κ ∂2u ∂x2 = κ 2∂3u ∂x3 = κ 3∂4u ∂x4 = . . . (63)
0 1 2 3 4 5 6 7 8 9 10 10−9 10−8 10−7 10−6 10−5 10−4 10−3 x (a) 0 1 2 3 4 5 6 7 8 9 10 10−10 10−9 10−8 10−7 10−6 10−5 10−4 x (b)
Figure 2: (a) absolute error εufor N=51; (b) absolute error εufor N=101. (N), 2ndorder; (H), 4thorder;
(△), model cd.a, f = 1
5;(▽), model cd.a, f = 1; (◦), model cd.b; (⋄), model cd.c; (), model cd.d; (+),
classical Richardson extrapolation; (×), differential Richardson extrapolation.
This artifact is only present in this simple linear problem. However, it allows us to in-vestigate the effects of broader stencils and their evaluation at the wall. Inconsistencies appear when evaluating the stencils at different grids (Ω∆1
and Ω∆2
) in near wall nodes, because of the use of excentric stencil formulations respecting the internal order of ac-curacy. As a consequence, an incorrect value for Cσ is generated in the first coarse grid
node near the wall and the wall node itself if n ≥ 3. Therefore these values cannot be used for interpolation to the fine grid. In order to avoid these ambiguities we extrapolate L, M and C in the compromised points. A piecewise cubic Hermite extrapolation is used. We like to emphasize that we intentionally propose the inconsistent model cd.c, and the Smagorinsky-like model cd.d. Simulations are performed on a uniform mesh with 50 cells and 100 cells. The adopted parameters are c = 0.1, ν = 0.1, L = 10 so that P e = 10. We define the rms error as
εu =
q
(u − uanalytic)2 (67)
The results in figure 2 are within the expectations for this simple canonical problem. The dynamic procedure with model cd.a and f = 1
5 obtains higher accuracy then the 4
thorder
solution and thus differential Richardson extrapolation (f = 0), close to 6th order. The
classical Richardson extrapolation, combining 2 resolved solutions on 2 different grids, seems to be less performant than 4th order accuracy. The analytical reduced model cd.b
reaches only 4th order, which is not surprising as the used stencils are smaller. Further it
Figure 3: Geometry of the square lid-driven cavity
is observed. Finally if f = 1, poor results are obtained with model cd.a, as was expected form the spectral analysis, although stability was surprisingly not affected.
5.2 2D lid-driven cavity
The driven cavity is also a benchmark test case. An internal recirculating flow is generated by a uniform moving wall in a 2D square closed domain (figure 3). The flow is representative for more complex situations with vortices, and is a challenging test case. No-slip Dirichlet boundary conditions for the velocity are imposed at all walls. The flow at the lid moves with the lid, leading to corner-singularities as result of the discontinuous boundary conditions. The pressure field is extrapolated at the wall using the 4th order
accurate Neumann condition ∂n∂3p3 with n the wall-normal direction, and the mean value is
kept at zero. The Reynolds number is Re = 400, for wich the flow is laminar. A pseudo-compressible code is used with a 3rd order Runge Kutta method for stepping in pseudo
time. Spatial discretization is 2nd order central. No artificial stabilisation was used in the
continuity equation. Therefore, the fields contain a minimal spurious pressure mode, that however does not affect the velocity results. The equations solved in the physical domain Ω = [0, L] × [0, L] ⊂ IR2, are: 1 β2 ∂p ∂t + ∂ui ∂xi = 0 (68) ∂ui ∂t + uj ∂ui ∂xj = −∂p ∂xi + ν∂ 2u i ∂x2 j (69) with β the artificial speed of sound. The Dirichlet boundary conditions are ui = 0 on Γ1
en ui = ulid on Γ2. For a uniform grid with spacing ∆x and ∆y the equations read:
The exact truncation errors (13-14) are given in previous paragraphs. We use a 6th order
central scheme on a 180 × 180 mesh, to generate the reference solution. The error for a variable φ is defined as
εφ= φref erence− φresolved (72)
The different simulations are given in Table 2.
Code Mass Momentum
ldc.2o60×60 2nd order central
ldc.4o60×60 4th order central
ldc.6o60×60 6th order central
ldc.1a60×60 exact truncation
ldc.1b60×60 exact truncation with f = 1
ldc.260×60 exact truncation no correction
ldc.360×60 no correction Smagorinsky
ldc.460×60 exact truncation Smagorinsky Table 2: Summary of the different simulations
Using an exact truncation model in combination with the dynamic procedure, a global least-square averaging was chosen for Cπ en Cσ. This is justified because of the uniform
grid and because the most important basic mode, is expected to be a the theoretical constant. For the Smagorinsky-like model however, we use a local least-square averaging over a 3 × 3 neighbourhood. Moreover, a clipping to Cσ > −0.2 for excessive negative
values of Cσ was necessary. For both models a blending with f = 1
5 as in formula (44)
was used, unless stated otherwise.
The results for the test sections x = L2 and y = L2, which are given in figures 4, correspond to the general expectations. Again it can be noticed that the classical Richardson extrap-olation (RE) does not reach 4th order accuracy, unlike the differential form (DRE), which
obtains almost 4th order. As expected, the dynamic procedure, with the exact truncation
model (ldc.1a60×60) obtains approximately 6th order accuracy. Surprisingly, in contrast
to the 1D test case, results of ldc.1b60×60 with f = 1, do not significantly deteriorate the
results. From simulations ldc.260×60, ldc.360×60 en ldc.460×60, we can conclude first that
correction of the continuity equation is crucial for overall accuracy in the driven cavity (mainly due to the large corner-pressure gradients of the driven lid). Secondly, we see that the Smagorinsky model for the momentum equations leads to noticeable improvement of accuracy, however it is less performant than the exact model.
6 CONCLUSIONS
Navier-0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.025 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 0.025 εu 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −2.5 −2 −1.5 −1 −0.5 0 0.5 1 1.5 2 2.5 x 10−3 εv 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −5 −4 −3 −2 −1 0 1 2 3 4 5 x 10−3 εp 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −6 −5 −4 −3 −2 −1 0 1 2 3 4 x 10−3 εu 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.03 −0.025 −0.02 −0.015 −0.01 −0.005 0 0.005 0.01 0.015 0.02 ε v 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −5 −4 −3 −2 −1 0 1 2 x 10−3 εp
Figure 4: Error levels section x = L2: (△), ldc.6o60×60; (), ldc.4o60×60; (◦), ldc.2o60×60; (▽),
ldc.RE60×60; (⊲), ldc.DRE60×60, (×), ldc.1a60×60;(+), ldc.1b60×60; (–.), ldc.260×60; (—), ldc.360×60; (
Stokes equations. The dynamic procedure applied to a 2nd order central scheme,
obtains a 6th order accurate solution.
2. The differential Richardson extrapolation is a special formulation of this technique. Moreover this point of view leads to an alternative implementation of the differential Richardson extrapolation.
3. The dynamic procedure, in combination with the exact truncation model for the mass and momentum equations leads to excellent results. Numerically 6th order
accuracy could be reached. The differential Richardson extrapolation reaches 4th
order accuracy.
4. The dynamic procedure, in combination with a Smagorinsky-like model for the truncation errors in the Navier-Stokes equations leads to an accuracy improvement within certain limitations.
5. The dynamic procedure for accuracy improvement can be used on any differential equation, including scalar transport equations.
ACKNOWLEDGMENTS
Research funded by a Ph.D grant of the Institute for the Promotion of Innovation through Science and Technology in Flanders (IWT-Vlaanderen).
REFERENCES
[1] S. Ghosal. An analysis of numerical errors in large-eddy simulations of turbulence. J. Comp. Phys., 125:187–206, January 1996.
[2] A. G. Kravchenko and P. Moin. On the effect of numerical errors in large-eddy simu-lations of turbulent flows. J. Comp. Phys., 131:310–322, September 1997.
[3] F. K. Chow and P. Moin. A further study of numerical errors in large-eddy simulations. J. Comp. Phys., 184:366–380, September 2003.
[4] O. Debliquy, B. Knaepen, D. Carati, and A. A. Wray. Sampling versus filtering in large-eddy simulations. In P. Moin and N. N. Mansour, editors, Proceedings of the Summer Program 2004, pages 133–144, Stanford, December 2004. Stanford University and Nasa Ames Research Center.
[5] B. Knaepen, O. Debliquy, and D. Carati. Large-eddy simulation without filter. J. Comp. Phys., 205:98–107, 2005.
[7] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot. A dynamic subgrid-scale eddy viscosity. Phys. Fluids A, 3(7):1760–1765, July 1991.
