9512.net

甜梦文库

甜梦文库

当前位置：首页 >> >> # A least square extrapolation method for the a posteriori error estimate of cfd and heat tra

A Least Square Extrapolation Method for the A Posteriori Error Estimate of the Incompressible Navier Stokes Problem

M. Garbey , and W. Shyy Department of Computer Science University of Houston Houston, TX, 77204, USA http://www.cs.uh.edu Technical Report Number UH-CS-04-03 November 3, 2004 Keywords: Partial Differential Equations, Least Square Method, Richardson extrapolation, a posteriori error estimate. Abstract A posteriori error estimators are fundamental tools for providing con?dence in the numerical computation of PDEs. To date, the main theories of a posteriori estimators have been developed largely in the ?nite element framework, for either linear elliptic operators or non-linear PDEs in the absence of disparate length scales. On the other hand, there is a strong interest in using grid re?nement combined with Richardson extrapolation to produce CFD solutions with improved accuracy and, therefore, a posteriori error estimates. But in practice, the effective order of a numerical method often depends on space location and is not uniform, rendering the Richardson extrapolation method unreliable. We have recently introduced [Garbey 13th international conference on domain decomposition and Garbey & Shyy JCP 2003] a new method which estimates the order of convergence of a computation as the solution of a least square minimization problem on the residual. This method, called least square extrapolation, introduces a framework facilitating multi-level extrapolation, improves accuracy and provides a posteriori error estimate. This method can accommodate different grid arrangements. The goal of this paper is to investigate the power and limits of this method via incompressible Navier Stokes ?ow computations.

Department of Computer Science, University of Houston, Houston, TX 77204, USA Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, FL 32611, USA

1

A Least Square Extrapolation Method for the A Posteriori Error Estimate of the Incompressible Navier Stokes Problem

M. Garbey , and W. Shyy

Abstract A posteriori error estimators are fundamental tools for providing con?dence in the numerical computation of PDEs. To date, the main theories of a posteriori estimators have been developed largely in the ?nite element framework, for either linear elliptic operators or non-linear PDEs in the absence of disparate length scales. On the other hand, there is a strong interest in using grid re?nement combined with Richardson extrapolation to produce CFD solutions with improved accuracy and, therefore, a posteriori error estimates. But in practice, the effective order of a numerical method often depends on space location and is not uniform, rendering the Richardson extrapolation method unreliable. We have recently introduced [Garbey 13 th international conference on domain decomposition and Garbey & Shyy JCP 2003] a new method which estimates the order of convergence of a computation as the solution of a least square minimization problem on the residual. This method, called least square extrapolation, introduces a framework facilitating multi-level extrapolation, improves accuracy and provides a posteriori error estimate. This method can accommodate different grid arrangements. The goal of this paper is to investigate the power and limits of this method via incompressible Navier Stokes ?ow computations. Index Terms Partial Differential Equations, Least Square Method, Richardson extrapolation, a posteriori error estimate.

I. I NTRODUCTION

AND

M OTIVATION

Richardson extrapolation (RE) is a simple, elegant and general mathematical idea that works for numerical quadrature with the Romberg method or ODE integrations that have smooth enough solution with the Bulirsch-Stoer method. Its use in Computational Fluid Dynamics (CFD) [3], [4], [10], [12], [13], [16], [17], [18], [21] is limited by the fact that meshes might not be ?ne enough to satisfy accurately the a priori convergence estimates that are only asymptotic in nature. Furthermore the order of convergence of a CFD code is often space dependent and eventually parameters, such as the Reynolds number, dependent. To cope with these limitations of RE, we have introduced recently [8], [9] the so-called Least Square Extrapolation method (LSE) that is based on the idea of ?nding automatically the order of a method as the solution of a least square minimization problem on the residual. Our LSE method is based on the post-processing of data produced by existing PDE codes. The method has been described in detailed in [9]. From a practical point of view, we have used a two dimensional turning point problem [14] exhibiting a sharp transition layer as well as a ?nite difference approximation of the cavity ?ow problem in ω ? ψ formulation [15] to show that our method is more reliable than RE while the implementation is still fairly easy and the numerical procedure inexpensive. Our objective is to use any PDE or CFD solvers, independent of their inner working algorithm and procedures, provided that they can offer the information including the residual of the numerical approximation, stability estimates, and varying grid resolutions and numerical solutions, to accomplish the following goals: (i) a posteriori estimates of PDEs that are more reliable and robust than straightforward Richardson extrapolation-based methods with low cost in additional CPU time, (ii) a solution with improved accuracy, (iii) arithmetic ef?ciency of the PDE multilevel solution procedure by providing a good starting point for iterative solvers [6], and (iv) a dynamic solution veri?cation software.

Department of Computer Science, University of Houston, Houston, TX 77204, USA Department of Mechanical and Aerospace Engineering, University of Florida, Gainesville, FL 32611, USA

2

From the applied mathematics point of view, a posteriori estimates have been around for many years [23] [1]. Most work has been done in the framework of ?nite element analysis on linear elliptic problems in order to drive adaptive mesh re?nement. More recently a general framework for ?nite element a posteriori error control that can be applied to linear and non-linear elliptic problem has been introduced by Patera et al [19]. A posteriori Finite-Element free constant output bounds can be constructed for the incompressible Navier Stokes equation [11]. We propose to use least square extrapolation to produce a posteriori estimate using grid solutions that can be produced by any discretization. This approach might be combined to existing a posteriori estimate when they are available, but is still applicable as a better alternative to straightforward RE when none such stability estimate is available. The extrapolation procedure is simple to implement and can be incorporated into any computer codes without requiring detailed knowledge of the source code. Its arithmetic cost should be modest compare to a direct computation of the ?ne grid solution. Finally the procedure should overall enhance the accuracy and trust of a CFD application in the context of solution veri?cation. In this paper, we pursue the research initiated in [8], [9] to investigate the power and limit of the LSE for the incompressible Navier Stokes equation written in the u ? v ? p formulation. This test case has a number of interesting features compared to the test case in the ω ? ψ formulation already studied in [9]. As a matter of fact, the ω ? ψ formulation is essentially a fourth order problem with one unknown, i.e the stream function ψ, while we will see that the LSE method applied to the u ? v ? p formulation has to deal with coupled equations, i.e, the momentum equation, and search for an extrapolation solution constrained by the divergence free condition. Further, for problem with multiple scale, the relation ship between the residual and the numerical error can be fairly complex. In other words minimizing the residual via LSE does not guarantee that the error is minimized, unless the space of search is close to the exact solution [22]. We will present in this paper a postprocessing procedure on our coarse grid solution that makes LSE numerically ef?cient. The plan of this paper is as follows. In Section 2, we ?rst summarize basic properties of RE and summarize the general idea of the LSE method for PDEs. In Section 3 we present in detail the algorithm for the incompressible Navier Stokes equation. In Section 4, we discuss the numerical results for the cavity ?ow problem. Section 5 is our conclusion and refer to ongoing research. II. BASIC P ROPERTIES

OF

R ICHARDSON E XTRAPOLATION

AND

L EAST S QUARE E XTRAPOLATION

Let E be a normed linear space, || || its norm, v ∈ E , p > 0, and h ∈ (0, h 0 ). ui ∈ E, i = 1..3 have the following asymptotic expansion, h ui = v + C( i?1 )p + δ, (1) 2 with C constant independent of h, and ||δ|| = o(h p ). For known p, RE formula, 2p ui+1 ? ui i vr = , i = 1, 2 (2) 2p ? 1

i provides improved convergence: ||v ? vr || = o(hp ). An a posteriori error estimate on ui is then i ||ui ? vr ||.

(3)

In CFD practice, one applies RE to grid functions rather than to continuous functions. Let E i be a family of normed linear space, associated with a mesh Mh/2i?1 . We suppose a set of equations,

U i = v + Ci ( h 2i?1 )p + δ i ,

(4)

with Ci = (1 + i )C, and i = o(1). δi is a model for the h independent numerical perturbation induced by consistency errors and/or arithmetic error. The Richardson extrapolate

Vr2 = 2p U 3 ? U 2 , 2p ? 1

(5)

3

de?ned on grid points of M2 has then for error in E2 ,

v ? Vr2 = 2p 1 ((δ2 ? 2p δ3 ) + C ( ?1

p

2

?

3 )(

h p ) ). 2

(6)

2 +1 The numerical perturbation is ampli?ed by a factor 2p ?1 . This RE gives then an a posteriori error estimate on U i that is simply ||Vr2 ? U i ||, i = 1..3.

RE can then be used also to approximate a ?ner grid solution, for example

U 4 ≈ Vr2 + C( h p ) , 23

(7)

where C is obtained from the identity

h p ) . (8) 22 For applications in CFD calculation, the asymptotic order of convergence is in general not known or not closely satis?ed on the computational grid. One may use the estimate: U 3 = Vr2 + C( p ? log 2 ||u1 ? u2 || . ||u2 ? u3 ||

(9)

An entirely similar analysis can be applied to non-embedded re?ned grid solution U i in a normed linear space Ei , associated with a mesh Mhi , provided that one projects all grid functions to a ?ne grid M 0 with an interpolation procedure. However this interpolation should introduce an additional error term integrated in the δ i term of (4) kept much less than the expected convergence accuracy h p . i In practice, all pointwise RE extrapolation formula, particularly (4), are sensitive to numerical perturbation. RE is a common tool for solution quality assessment in CFD. In our experience [9] [21], we have observed that RE can improve the order of accuracy, but not consistently. If the quality of the solution is poor then RE may provide worse approximations. These conclusions are reached based on extensive solution veri?cation with two different Navier Stokes approximations for the steady state, 2-D laminar incompressible lid-driven square cavity ?ow with the Reynolds number (Re) in the range of 20 to 1000. Squared regular meshes using the ω ? ψ formulation and Finite Difference (FD) [15] or the Finite Volume (FV) version of the u ? v ? p formulation with centered cells [20] have been tested. Further experiments with turbulent ?ows on a back step has demonstrated the critical issue of multiscales [21]. In a recent stream of work of Eca et al, see [4], [5] and its references, one can use a least square model of the error provided that enough grid solutions are computed. To be more speci?c most variant of RE suppose that the error is represented by an asymptotic expansion as follows:

e = v ? ui = Σk=1..p ak βk (h) + o(hp ), i

(10)

in a normed space (Ei , || ||). In RE procedure, one neglects the o(hp ) residual and use the following identity pointwise: i

e = v ? ui = Σk=1..p ak βk (h),

(11)

We observe that this pointwise equality not only neglect the higher order term o(h p ) but also disregards the nature i of the asymptotic expansion that is depending on the norm associated to E i . We have then two standard situations: ? If the basis function β are given, then one need p + 1 grid solutions v j to derive the unknown coef?cients ak , , k = 1..p and the (approximated) true solution v. γ ? If the set of basis functions βk (h) is a one parameter family of functions, for example h k , with arbitrary exponent γ, one needs 2 p + 1 grid solution to solve the error model, i.e compute v, a k , k = 1..p, and γk , k = 1..n, pointwise. To cope with the fact that RE is very sensitive to noisy data, L.Eca et Al retrieve the error model with a least square ?t instead of enforcing equalities. This procedure is indeed less sensitive to noisy data, but requires many more grid solutions than with the standard RE procedure.

4

We have developed a completely different technique to optimize RE. Our criterium to select the best extrapolate solution is to minimize an objective function such as the l 2 norm of the residual for the discrete solution on a very ?ne grid. This ?ne grid must be chosen to resolve the ?ne scale of the problem. We use no more than two or three coarse grids solution in our procedure. We do not try to compute directly an approximation of the exact solution either. We rather try to extract the best information from these two or three coarse grid solutions by reintroducing in the construction of the extrapolation formula, the discretization of the PDE, instead of assuming any kind of asymptotic model for the error. Let us review brie?y the LSE method for the numerical approximation of scalar function ?rst. 1 2 Let E = L2 (0, 1), u ∈ E. Let vh and vh be two approximations of u in E:

1 2 vh , vh → u in E as h → 0.

A consistent linear extrapolation formula writes

1 2 αvh + (1 ? α)vh = u.

In RE the α function is a constant. In the LSE method we formulate the following problem for the unknown function α that is in general a non-constant function 1 2 Pα : Find α ∈ Λ(0, 1) ? L∞ such that (α vh + (1 ? α) vh ? u) is minimum in L2 (0, 1). Typically we choose for the space Λ(0, 1) a set of polynomial trigonometric functions, but this is not necessary. We have shown

i Theorem [9]: if u, vh , ∈ C 1 (0, 1), i = 1, 2 , if an 0(M ?1 ) × 0(hp ) approximation of u. 1 1 2 vh ?vh 2 1 1 2 ∈ L∞ (0, 1) and vh ? vh = 0(hp ) then αvh + (1 ? α)vh is

1 2 2 Special care must be done if vh ? vh << u ? vh , in some set of non-zero measure. These outliers should not affect globally the least square extrapolation and we impose α to be a bounded function independent of h. A potentially more robust approximation procedure consists of using three levels of grid solution as follows: 1 2 3 Pα,β : Find α, β ∈ Λ(0, 1) such that (α vh + β vh + (1 ? α ? β) vh ? u) is minimum in L2 (0, 1). j As a matter of fact, all vh , j = 1..3, may coincide at the same grid points only if there is no grid convergence locally. In such a situation, one cannot expect improved local accuracy from any extrapolation technique. The robustness of the LSE method comes from the fact that the extrapolated solution does not deteriorate the accuracy of the coarse grid solution while it may not be the case for RE, especially when one uses (9). In practice, we work with grid functions solution of discretized PDE problem. The idea is now to use the PDE in the RE process to ?nd an improved solution on a given ?ne grid M 0 . Let us denote formally the linear PDE

L[u] = f, with u ∈ (Ea , || ||a ) and f ∈ (Eb , || ||b ),

and its numerical approximation,

Lh [U ] = fh , with U ∈ (Eh , || ||a ) and fh ∈ (Eh , || ||b ), a b

parameterized by a mesh step h. We suppose that we have a priori a stability estimate for these norms

||U ||a ≤ C hs (||fh ||b ),

(12)

with s real not necessarily positive. Let Gi , i = 1..3, be three embedded grids that do not necessarily match, and their corresponding grid solutions ? Ui . Let M 0 be a regular grid that is ?ner than the grids Gi . Let Ui be the coarse grid solutions interpolated on the 0. ?ne grid M

5

The main idea of the LSE method is to look for a consistent extrapolation formula based on the interpolated ? ? coarse grid solutions Ui that minimizes the residual, resulting from Ui on a grid M 0 that is ?ne enough to capture a good approximation of the continuous solution. Let us restrict for simplicity to a two-point boundary value problems in (0, 1). Our least square extrapolation is now de?ned as follows: ? ? Pα : Find α ∈ Λ(0, 1) ? L∞ such that (Lh [αU 1 + (1 ? α)U 2 ] ? fh ) is minimum in L2 (M 0 ). The three-level version is analogous to the two-level one. To focus on the practical use of this method, we should make the following observations. It is essential that the interpolation operator gives a smooth interpolant depending on the order of the differential operator and the regularity of the solution of the differential problem. For conservation laws, one may require that the interpolation operator satis?es the same conservation properties. For reacting ?ow problems, one may require that the interpolant preserves the positivity of species. For elliptic ? problems, it is convenient to postprocess the interpolated functions U i , by few steps of the relaxation scheme

V k+1 ? V k ? = Lh [V k ] ? fh , V 0 = U i , (13) δt with appropriate arti?cial time step δt. This will readily smooth out the interpolant. Let ej , j = 1..m be a set of basis function of Λ(0, 1). The solution process of P α and/or P(α,β) can be decomposed into three consecutive steps. ? First, interpolation of the coarse grid solution from G i , i = 1..3 to M 0 . ? ? ? ? Second, evaluation of the residual L h [ej (U i ? U i+1 )], j = 1..m, and Lh [U 3 ] on the ?ne grid M 0 . ? Third, the solution of the linear least square problem that has m unknowns. In practice, we keep m low by using a spectral representation of the unknown weight functions α and eventually β. The arithmetic complexity of the overall procedure is then still of order Card(M 0 ), i.e., it is linear. The application to nonlinear PDE problem is done via a Newton-like loop [9]. The algorithm might be coded in a stand alone program independent of the main code application. We are going now to describe the algorithm for the incompressible Navier Stokes set of equations.

III. A PPLICATION

TO THE CAVITY FLOW PROBLEM

Let us consider the velocity-pressure formulation of the square cavity problem in two space dimensions. The steady problem writes in ? = (0, 1)2 ,

1 ?u ?u ?p ?u + u +v + = 0, (x, y) ∈ ? (14) Re ?x ?y ?x ?v ?v ?p 1 +v + = 0, (x, y) ∈ ?, (15) N2 [u, v, p] = ? ?v + u Re ?x ?y ?y submitted to the constraint ?u ?v Div(u, v) = + = 0, (x, y) ∈ ?. (16) ?x ?y In this system of equations Re is the Reynolds number. Furthermore this set of equations is supplemented with the no-slip boundary conditions on the walls of the cavity. The ?ow speed is zero on all walls except on the sliding wall u(x, 1) = g(x), x ∈ (0, 1). (17) N1 [u, v, p] = ?

In applying the LSE method with u ? v ? p formulation, we deal with three new dif?culties that were not present in the Navier Stokes calculation of [9] ? We have a system of coupled non linear PDEs. The cavity ?ow problem with ? Ψ formulation is really a fourth order non linear elliptic problem on Ψ only. ? The LSE on the velocity ?eld should satisfy the divergence free constraint. ? Thanks to the discontinuity of the boundary condition on the velocity ?eld at the corners, there is no valid pointwise model of the error that follows a standard Taylor expansion. The grid functions (ui , vi , pi ) on Gi are computed with a standard FD code using a projection method and staggered grids [15].

6

Let us consider a set of three-grid solutions (? i , vi , pi )i=1..3 projected onto the ?ne grid M 0 via a high order u ? ? 0 0 smooth interpolation procedure. Let us denote N1 , N2 , Div 0 the corresponding discretized operator. For ?nite differences that we will consider from now on, M 0 is a staggered grid system, and the discretized operator are given by central second order ?nite differences. The projected ?ow ?eld (?i , vi ) does not satisfy a priori the divergence-free condition u ?

Div 0 (?i , vi ) = 0, (x, y) ∈ ?. u ?

(18)

In the unlikely case where (18) is satis?ed, the extrapolated value of the ?ow ?eld

α(?i , vi ) + (1 ? α)(?i , vi ), i = j, u ? u ?

will not be divergence-free anyway for the Div 0 operator, unless α is a constant. We de?ne then the following mapping (u, v) → Ψ → (U, V ), where ?Ψ ?Ψ U= , V =? , (x, y) ∈ ?, and ?y ?x ?u ?v ?Ψ = ? , (x, y) ∈ ?, Ψ = 0 on ??. ?y ?x ? ? We de?ne also its discrete analogue (?, v ) → ( U , V ), in M 0 . u ? Thanks to this mapping, we can retrieve from the grid function (u i , vi ) in Gi a divergence free approximation ? ? (Ui , Vi ) in M 0 . The least square extrapolation problem with two levels writes then Pα1 ,α2 : Find α1 and α2 ∈ Λ(?) ? L∞ (?) such that

? ? N 0 [α1 Ψ1 + (1 ? α1 )Ψ2 , α2 p1 + (1 ? α2 )?2 ] ? p

0 0 is minimum in L2 (M 0 , M 0 ), with N 0 [Ψ, p] = (N1 [U, V, p], N2 [U, V, p]).

Since this problem is nonlinear, we use a Newton loop to construct a sequence of weight functions (α n , αn ) 1 2 that may converge to the solution. The iterative procedure starts from the ?nest coarse grid solution at our disposal. Convergence is not guaranteed and may depend on how close the initial guess is to the true M 0 grid solution. If (U 0 , V 0 , p0 ) represents the current solution, the next iterate is found by applying the least square extrapolation procedure to the linear operator

L0 (U 0 , V 0 )[Ψ, p] = (L0 (U 0 , V 0 )[U, V, p], L0 (U 0 , V 0 )[U, V, p]), 1 2

(19)

with

L1 (U 0 , V 0 )[u, v, p] = ? ?U 0 ?U 0 ?p ?u ?u ?U 0 ?U 0 1 ?u + U 0 +u +V0 +v + ? U0 ?V0 , Re ?x ?x ?y ?y ?x ?x ?y L2 (U 0 , V 0 )[u, v, p] = 1 ?V 0 ?V 0 ?p ?v ?v ?V 0 ?V 0 ?v + U 0 +u +V0 +v + ? U0 ?V0 , Re ?x ?x ?y ?y ?y ?x ?y A similar algorithm is derived for the three-level case. The space of unknown weight function is chosen as in [9] to be the set of trigonometric polynomial functions ? α = Σi=1..m,j=1..mαi,j ei ej ,

with e0 = 1, e1 = cos(πx) and ei = sin((i ? 2)πx), for i = 3..m. This set of trigonometric functions allows us to approximate at second order in L 2 norm any smooth nonperiodic functions of C 1 [(0, 1)2 ], [7]. The main advantage of this choice of approximation space for the weight function is that it will allow us to easily interpret our numerical result in the frequency space. However, for Navier Stokes computation with large Reynolds number, we are currently investigating the use of wavelets since the convergence order might be closely related to the multiscale properties of the solution. We are going now to present our numerical experiments with the LSE method.

7

IV. R ESULTS

AND

D ISCUSSION

To illustrate the numerical result, we restrict ourselves to the test case of the square cavity with a constant sliding wall velocity that is g(x) = ?1, and a Reynolds number Re = 400. This test case is representative of the results obtained with our method. In particular, the ?rst component u of the speed is singular at the corner, as well as the pressure. From this numerical experiment and many others we can draw the following conclusions: The three-level extrapolation method is more robust and more accurate than the two-level extrapolation method. Figure 1 illustrates the cancelation phenomenon with two grid solutions 51 × 51 and 61 × 61. We plot in this picture the local minimum per vertical and/or horizontal lines of the difference between two coarse grid solutions projected on the ?ne grid. These minima are such that any a priori bounded weight coef?cient α will have no in?uence on the extrapolated solution. The two-level extrapolated solution cannot therefore improve the accuracy of the solution. Further we cannot decide if the grid solutions are fully converged at these points or, on the contrary, if the numerical methods lack convergence locally, unless we use a very grid solution to compare. We have veri?ed then that the RE formula based on three-level extrapolated solution with the computed convergence order as in (9) exhibits large local errors on the convergence order. Following the approach of [3], one may select adaptively the points where RE has better chances to be valid. We are going to present the numerical experiments for LSE using three grids only. In the following we will present a fair comparison between RE and LSE, i.e we will compare RE and LSE to predict the same ?ne grid solution, instead of using RE to predict the continuous true solution. In each graph (2) to (6) the horizontal axis indicates the number of grid points N in each space direction for the ?ne grid M 0 used to evaluate the residual. The vertical axis gives in log 10 scale the relative error in L2 norm’s. Labels of curves are as follows: ’o’ for the grid solution G2 solution, ’v’ for the the ?nest coarse grid solution G 3 solution, for RE, for LSE. We recall that M 0 is the ?ner grid on which the extrapolation is conducted. Let us take M 0 to be a grid with a slightly smaller space step than G3 . RE as well as LSE predicts very accurately the solution on M 0 . In ?gure 2, the three coarse grid solutions are 51 × 51 for G1, 61 × 61 for G2, and 71 × 71 for G3. We see that as M 0 gets ?ner the LSE deteriorates, while the RE assuming second order improved. In general LSE seems to be reliable to bounds from below the true error on G 3 by comparing the prediction 0 done by LSE and the coarse grid solution G3 interpolated on M 0 . This seems to be a promising tool for on M routine solution veri?cation. RE gives similar performance for this speci?c benchmark, and gives much better result than LSE for ?ner grid prediction. We observe a numerical locking of the LSE method to predict solution when the grid M 0 gets signi?cantly ?ner than G3 . In other words increasing the number of Fourier modes m to approximate the weight function provides little improvement on the minimum of the residual. Further, the optimal weight function does not correspond necessarily to the minimum of the residual. This phenomenon has been clearly demonstrated in [22], where LSE was restricted to the search for constant α values with two grid levels only. In this speci?c case, the best weight coef?cient that minimizes the residual in the least square sense, can give poorer result than the RE method assuming second order of convergence. Figure 3 illustrates the poor performance of the LSE method with three grid levels and m = 4 Fourier modes in each space direction to approximate the weight function. In this ?gure, we have compared the solution with the 121 × 121 grid solution declared as the true solution, while the LSE uses a ?ne grid M 0 with growing size from 81 × 81 to 121 × 121. The three coarse grids solution are still 51 × 51 for G1, 61 × 61 for G2, and 71 × 71 for G3. We propose the following explanation of this phenomenon: we observe that the coarse grid solution interpolated on a ?ne grid M 0 has high spurious wave number terms brought by the interpolation procedure. We know that the high wave number components of the coarse grid solution relative to the coarse grid itself are inaccurate. The interpolation procedure combined with the weight function expansion worsen the phenomenon. These high wave numbers components are ampli?ed in the computation of the residual for a nonlinear stiff problem as the cavity ?ow with large Reynolds number. The LSE method minimizes therefore the L 2 norm of a residual polluted by high wave number components. To get the minimum of this residual does not guarantee therefore that the error on the low wave number component of the solution is minimum. This phenomenon is not visible when LSE is used to predict the solution on a near by G3 ?ner grid, because the gap in frequency between G 3 and M 0 is small.

8

To validate our heuristic analysis, we postprocess each coarse grid solution with the NS code on the ?ne grid. We use explicit time stepping with dt constraint by the CFL condition as well as the explicit treatment of the diffusion term. Ten time steps does not allow the NS to converge on M 0 by all means, but relax ef?ciently the high frequency components of the projected coarse grids due to the interpolation. We further apply a least square low mode approximation of the computation of each residual computed in the LSE method. LSE is therefore now computing the weight functions that minimize the low mode approximation of the residual. In other words high frequency components are completely ?ltered out, and the effect of the singularity at the corner in the computation of the residual is somehow weakly weighted. The same test case as the one illustrated in Figure 3 has been used. Figure 4 shows that keeping a 8 Fourier modes in each space direction for the residual approximation, improves signi?cantly the result. This result is fairly insensitive to the number of NS iterates on the ?ne grid, once the spurious oscillations introduced by the interpolation of the coarse grids are damped out on M 0 . We checked for example that to take 100 NS iterates instead of 10 improve the accuracy of LSE marginally only. The same result can be reproduced with higher accuracy for ?ner meshes as in Figure 5. In this last case the declared true solution is for the grid 181 × 181 and the three coarse grid solutions are 101 × 101, 111 × 111 and 121 × 121 grids. It demonstrates some practical convergence of the LSE method with very good error estimate on the solution. Further LSE gives the best performance to predict the grid solution on M 0 when the coarse grid solution are projected on the same M 0 . However, how good should be the coarse grid solution is still an open issue. Finally, it should be noticed that as the coarse grid solutions get ?ner, the LSE accuracy is always signi?cantly better than the RE prediction. This is shown in Figure 6 where LSE is computed with simultaneous increasing resolution of the coarse grid solutions (N ? 20)2 for G1, (N ? 10)2 for G2 and N 2 for G3, for N = 70 up to N = 110. Finally it can be observed on this test case that the ?ow speed is discontinuous at the two corners of the sliding wall. This singular behavior of the pressure and the velocity components at these corner points leads to locally low order accuracy of the numerical solution [21]. This impacts the ef?ciency of the RE indeed. We obtain then better a posteriori estimate with LSE than with RE, thanks to the postprocessing of the coarse grid solution. It is not yet clear if a wavelet representation of the weight function α and β, will better approximate sharp variation of the convergence order at the corner and give a signi?cantly better result for the singular case obtain with g(x) = ?1. We have also observed that spline interpolation of the coarse grid solution on M 0 smears out the singularity at the corner and might be also one of the barriers for recovering accurate solution in the L 2 norm with very coarse grids. This is currently the subject of further investigations. V. C ONCLUSIONS We have studied a new extrapolation method for PDEs that is more robust and accurate than RE applied to numerical solutions with inexact or varying convergence order. Our method provides a better tool to establish a posteriori error estimate than RE when the convergence order of a CFD code is space dependent. However there are still many open questions concerning mainly how ?ne should be the coarse grid solution to provide accurate a posteriori error estimate. We are currently investigating the use of wavelet approximation instead of regular trigonometric polynomial to track the multiscale properties of the solution re?ected in sharp variation of the convergence order. Further let us mention that there are many variants of the least square method [2]. We have so far considered the most straightforward method with an unreliable estimator. We may therefore need to ?nd an optimal weight to the Least square, and eventually use better methods such as the generalized least square or nonlinear least square method. Acknowledgment:We would like to thank one of the referees for pointing out the fact that one can use RE to estimate a ?ne grid solution. The work of M.Garbey was sponsored by Sandia Nat. Lab. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000. Part of the work of W.Shyy was sponsored by NASA.

9

R EFERENCES

[1] M.Ainsworth and J.T.Oden, A Posteriori Error Estimation in Finite Element Analysis, Wiley Press, 2000 [2] A.Bjorck, Numerical Method for Least Squares Problems, SIAM, 1996. [3] G.De Vahl Davis, Natural Convection of Air in a Square Cavity: a Bench Mark Numerical Solution Int. J. for Num. Methods in Fluids, Vol. 3, pp 249-264, 1983. [4] L.Eca and M.Hoekstra, An Evaluation of Veri?cation Procedures for CFD Applications, 24th Symposium on Naval Hydrodynamics, Fukuoka, Japan, 8-13 July 2002. [5] L.Eca and M.Hoekstra, Veri?cation Procedures for CFD, IST Report D72-14, 2002. [6] A.Ecer, M.Garbey and M.Hervin, On the Design of Robust and Ef?cient Algorithms that Combine Schwarz Method and Multilevel Grids 12th Int. Conf. Parallel CFD2000, , C.B.Jenssen et al Edt. Elsevier, pp 165-172, 2000. [7] D. Gottlieb and C.W.Shu, On the Gibbs Phenomenon and its Resolution, SIAM Review, Vol 39, No4,pp 644-668, 1997. [8] M. Garbey, Some Remarks on Multilevel Method, Extrapolation and Code Veri?cation 13th Int. Conf. on Domain Decomposition DD13, Domain Decomposition Methods in Science and Engineering, CIMNE, Barcelona, N.Debit et Al edt, pp379-386, 2002. [9] M.Garbey and W.Shyy, A Least Square Extrapolation Method for improving solution accuracy of PDE computations, J. of Comput. Physic, 186, pp1-23, 2003. [10] A.G.Hutton and M.V.Casey, Quality and Trust in Industrial CFD - A European Initiative, 39th AIAA Aerospace Sciences Meeting, 8-11 January, 2001/ Reno, NV, AIAA Paper 2001-0656. [11] L.Machiels, J.Peraire and A.T.Patera, A Posteriori Finite Element Output Bounds for the Incompressible Navier Stokes Equations: Application to a Natural Convection Problem, JCP, Vol 172, No 2, 2001. [12] W.L. Oberkampf, F.G. Blottner and D.Aeshliman, Methodology for Computational Fluid Dynamics Code Veri?cation and Validation, 26th AIAA Fluide Dynamic Conference, June 19-22, 1995/ San Diego, CA, AIAA Paper 95-2226. [13] W.L.Oberkampf and T. G.Trucano, Veri?vation and Validation in Computational Fluid Dynamics, Sandia Report 2002-0529, March 2002. [14] R.J.O’Malley, Singular Perturbation Methods for ODE, Applied Math. Sciences, Vol 89, Springer Verlag 1991. [15] R.Peyret and T.T.Taylor, Computational Methods for Fluid Flow, Springer series in Computational Physics, 1985. [16] P.J.Roache, Veri?cation and Validation in Computational Science and Engineering, Hermosa Publishers, Albuquerque, New Mexico, 1998. [17] C.J.Roy, M.A.McWherter-Payne and W.L.Oberkampf, Veri?cation and Validation for Laminar Hypersonic Flow?elds, AIAA2000-2550, Fluids 2000 Conference, Denver, CO, 2000. [18] C.J.Roy and M.M.Hopkins, Discretization Error Estimates Using Exact Solution to Nearby Problems AIAA2003-0629, Fluids 2003 Conference, Denver, CO, 2003. [19] J.Sarrate, J.Peraire and A.Patera, A Posteriori Finite Element Error Bounds for Nonlinear Outputs of the Helmholtz Equation, to appear in Int. J. Numer. Meth. in Fluids. [20] W.Shyy, S.S.Thakur, H.Ouyang, J.Liu and E.Blosch, Computational Techniques for Complex Transport Phenomena, Cambridge University Press 1997. [21] W. Shyy, M. Garbey, A. Appukuttan and J. Wu, Evaluation of Richardson Extrapolation in Computational Fluid Dynamics Numerical Heat Transfer, Part B: Fundamentals 41 (2): 139-164, 2002. [22] R.Vaidyanathan, W.Shyy, M.Garbey and R.Haftka, CFD Code Veri?cation Using Least Square Extrapolation Method, AIAA Reno 2004. [23] R.Verfurth, A Review of A Posteriori Estimation and Adaptive Mesh Re?nement Techniques, Wiley-Teubner, 1996.

10

120 100 80 60 40 20

From pressure on G1 and G2 100 80 60 40 20 20 40 60 80 100 120

From stream function on G1 and G2

20

40

60

80

100

The difference of pressure between grid solutions G1 and G2 cancels along these curves

Same phenomenon for the difference in stream functions

Fig. 1. places.

? ? Location of some local minima of p2 ? p1 on the left and ψ2 ? ψ1 on the right where cancellation with the two level’s LSE take ? ?

?1 ?1.1 ?1.2 ?1.3 ?1.4 ?1.5 ?1.6 ?1.7 ?1.8 ?1.9 ?2 80 o for G2 v for G3 square for LSE diamond for RE 85 90 95 100 105 Pressure 110 115 120 125

?1.4

?1.4

?1.6

?1.6

?1.8

?1.8

?2

?2

?2.2

?2.2

?2.4

?2.6 80

o for G2 v for G3 square for LSE diamond for RE 85 90 95 100 U 105 110 115 120 125

?2.4

?2.6 80

o for G2 v for G3 square for LSE diamond for RE 85 90 95 100 V 105 110 115 120 125

Fig. 2. Performance of LSE (without postprocessing of the residual) to predict the solution on the grid 121 × 121. The coarse grid solution G1 , G2 , G3 , are respectively 51 × 51, 61 × 61, 71 × 71.

11

?1

?1.2

?1.4

?1.6

?1.8

?2 o for G2 v for G3 square for LSE diamond for RE 85 90 95 100 105 Pressure 110 115 120 125

?2.2 80

?1.4

?1.4

?1.6

?1.6

?1.8

?1.8

?2

?2

?2.2

?2.2

?2.4

?2.4

?2.6 o for G2 v for G3 square for LSE diamond for RE 85 90 95 100 U 105 110 115 120 125

?2.6 o for G2 v for G3 square for LSE diamond for RE 85 90 95 100 V 105 110 115 120 125

?2.8 80

?2.8 80

Fig. 3.

Performance of LSE to predict the M 0 grid solution with the same notation and computation as in Figure 2.

12

?1

?1.2

?1.4

?1.6

?1.8

?2

?2.2 80

o for G2 v for G3 square for LSE diamond for RE 85 90 95 100 105 Pressure 110 115 120 125

?1.4 ?1.6 ?1.8 ?2 ?2.2 ?2.4 ?2.6 ?2.8 ?3 ?3.2 80 o for G2 v for G3 square for LSE diamond for RE 85 90 95 100 U 105 110 115 120 125

?1.4 ?1.6 ?1.8 ?2 ?2.2 ?2.4 ?2.6 ?2.8 ?3 ?3.2 80 o for G2 v for G3 square for LSE diamond for RE 85 90 95 100 V 105 110 115 120 125

Fig. 4. Performance of LSE to predict the M 0 grid solution with the same notation and computation as in Figure 2 but with postprocessing of the coarse grid solution and the residual to ?lter out the high waves components.

13

?1.4

?1.6

?1.8

?2

?2.2

?2.4 o for G2 v for G3 square for LSE diamond for RE 130 ?2 140 150 Pressure 160 ?2 170 180

?2.6

?2.5

?2.5

?3

?3

?3.5

?3.5

?4 o for G2 v for G3 square for LSE diamond for RE 130 140 150 U 160 170 180

?4 o for G2 v for G3 square for LSE diamond for RE 130 140 150 V 160 170 180

?4.5

?4.5

Fig. 5. Performance of LSE, with postprocessing of the residual, to predict the solution on the grid 181 × 181. The coarse grid solution G1 , G2 , G3 , are respectively 101 × 101, 111 × 111, 121 × 121.

14

?1

?1.5

?2

?2.5

?3

?3.5

?4 65

o for G2 v for G3 square for LSE diamond for RE 70 75 80 85 90 Pressure 95 100 105 110 115

?1.5

?1.5

?2

?2

?2.5

?2.5

?3

?3

?3.5

?3.5

?4

?4

?4.5

?4.5

?5

?5.5 65

o for G2 v for G3 square for LSE diamond for RE 70 75 80 85 90 U 95 100 105 110 115

?5

?5.5 65

o for G2 v for G3 square for LSE diamond for RE 70 75 80 85 90 V 95 100 105 110 115

Fig. 6. Comparison of LSE versus RE with varying accuracy for the coarse grid solution. G1 , G2 , G3 , are respectively (N ? 20)2 , (N ? 10)2 , N 2 grids. Horizontal axis gives the number of grid points N in each space direction for G3.

赞助商链接

- The Method of Least Squares
- MethodLeastSquares
- Steganalysis of two least significant bits embedding based on least square method
- Power system load forecasting using partial least square method
- CONSTRAINED OPTIMIZATION OVER DISCRETE SETS VIA SPSA WITH APPLICATION TO NON-SEPARABLE RESO
- A posteriori error estimate for the mixed finite element method
- Steganalysis of two least significant bits embedding based on least square method
- A CFD-based test method for control of indoor environment and space ventilation_虚拟传感器,PID控制
- CFD simulation of coupled heat and mass transfer through porous foods during vacuum cooling process
- CFD Flow and Heat Transfer Simulation for Empty and Packed Fixed Bed Reactor in Catalytic C
- CFD modeling of flow and heat transfer in a thermosyphon
- CFD MODELLING OF MIXING AND HEAT TRANSFERe
- Analysis of fluid flow and heat transfer in a rib grit roughened surface solar air heater using CFD
- CFD analysis of heat collection in a glazed gallery
- Power system load forecasting using partial least square method

更多相关文章：
更多相关标签：