- Email: [email protected]
A local mesh-refinement technique for incompressible flows
Computers & Fluids Vol. 14, No. I, pp. 69-81, 1986
0045-7930/86 $3.00 + .IX) © 1986 Pergamon Press Ltd.
Printed in Great Bdtain.
A LOCAL MESH-REFINEMENT TECHNIQUE FOR INCOMPRESSIBLE FLOWS L. FUCHS Delmrtment of Gasdyuamics, The Royal Institute of Technology, Stockholm, Sweden
(Received 26 November 1984; in revisedform 8 February 1985) Abstract--As part of a Multi-Grid scheme for the solution of the Navier-Stokes equations in primitive variables, we introduce a local mesh refinement procedure. New cartesian sub-grids are introduced into regions where the estimated truncation errors are too large. Through the Multi-Grid processing, informations is transferred among the grids in a stable and efficient manner. A simple pointer system allows the storage of the dependent variables, without increasing in the required computer memory. Two computed examples of incompressible flow problems are discussed.
INTRODUCTION
One of the complicating factors in the solution of the Navier-Stokes equations is due to the existence of different (length) scales in the flow field. The dependent variables are made dimensionless by some (globally valid) characteristic scale. This scale may not be represen= tative in all parts of the flow field and locally one may have much smaller or larger characteristic scales. The accuracy of finite-differences (or any other discrete approximation) depends not only on the formal order of approximation but also on the ratio of the mesh size (h) and the local scale (L). For accuracy, one has to resolve the local scales so that the variations in the dependent variables are described on as many mesh intervals as possible (h/L = o(1)). Generally, it is impractical to use uniformly fine meshes because of computer memory restrictions and long computational times. Furthermore, such fine meshes are not needed for accuracy, in larger parts of the flow field. A frequently used method is to construct a non-uniform (body fitted) grid so that the mesh spacing is small in regions with expectedly large gradients (e.g. [1, 2]). The main disadvantage of this method is that one has to know, a priori, the details of the solution when the mesh is generated. Furthermore, the transformation of the coordinates requires additional storage (of at least the location of the node points) and additional computational operations (on the additional terms in the governing equations). Alternative approaches to global transformation of coordinates are grid patching, mesh embedding and adaptive grid generation using either orthogonal or non-orthogonal grid systems. The method of grid patching uses the idea of dividing the domain into sub-domains which accommodate independently generated meshes. These grids are joined together along the boundaries of the sub-domains. This method has been used for the computation of some transonic flows [5, 6]. The approach requires careful treatment of the information transfer across the boundaries of the sub-domains. Mesh embedding can be considered as a special case of patching with constraints imposed on the grids on both sides of the 'internal' (sub-domain) boundaries. Embedded meshes have, usually, the same basic structure as the rest of the mesh in which they are being embedded. This method has also been applied to transonic problems with success [7, 8]. A more generalized mesh embedding/patching has been proposed recently by Benek, Steger and Dougherty [9]. Both in mesh patching and mesh embedding techniques the only interaction among the different sub-domains (sub-grids) is done by updating (transferring) boundary data. Such procedures reduce in general, the efficiency and in some cases also the stability properties of the numerical scheme. Adaptive schemes are the most desirable for obtaining top efficiency in terms of accuracy. On the other hand such optimal grids requires by themselves considerable amount of corn69 CAF 14:1-F
70
L. FUCHS
putational effort. For I-D problems (or even for time-dependent 1-D cases) several methods have been developed [10-12]. An adaptive grid generation procedure for multidimensional problems has been developed and described by Dwyer, Kee and Sanders [ 13]. Their adaptive scheme is based upon the placement of grid points in proportion to the gradients in the dependent variables. Here, we do not attempt to produce 'optimal' grids. Instead, we use a system of cartesian grids, arranged in a grid hierarchy. New sub-grids are introduced, at places where the estimated truncation error is too large. In this scheme, the new local grids form new levels in the hierarchy of grids. This approach is close in spirit to the approach of Berger and Oliger [14]. However, by realising that adaptive grid generation procedures and the solution algorithms are strongly connected, we construct the grid-system dynamically during a MultiGrid (MG) solution procedure. In fact, several elements of the current scheme have been proposed by Brandt [ 15]. Recently, the combination of mesh embedding with a MG algorithm have been reported by Usab and Murman [16] and by Brown [17]. In both papers the emphasis is on the gains in computational time compared to global mesh refinement computations, and no attempt is made to control the location or the size of the locally refined grids. In the present work we study the use of locally refined meshes as an integral part of a MG scheme for the solution of incompressible flow problems. As a further step toward a fully adaptive grid generation, we study the behavior of the truncation errors for different levels of refinement. The values of the appropriate dependent variables which must be specified on the boundaries of the sub-grids are computed by interpolations from the coarser grids. These interpolation formulae must, however, be such that the mass fluxes are conserved. In non-MG schemes the interaction among the patched/embedded grids takes place only by the updating procedure of the (internal) boundary conditions (see e.g. [5, 9, 14]). In full MG relaxations, information is passed from the (global) coarser grids to the locally refined grids by the interpolated corrections and the 'internal boundary conditions'. Information, which gives rise to improved global accuracy, is passed back from the fine local grid to the coarser grids by adding a truncation error related term to the equations. The scheme can identify regions with large truncation errors and by using local rectangular meshes, these errors can be reduced to a desirable (uniform) level. The administration of the local grids is very simple and requires minimal 'book-keeping'. In contrast to non-uniform, 'optimal' grids, one has to store only six integer values for each sub-grid. In the following we describe some details of the numerical scheme in the MG context. The method has been applied to the computation of the flow in the square driven cavity and a simplified separator model. By using very fine local meshes near the (singular) comers of the driven cavity, it is found that the truncation error has a singular behavior (O(h-2)) which is in agreement with the Stokes flow in such regions. In both cases, however, the level of the truncation errors can be controlled very efficiently by using locally refined grids. GOVERNING EQUATIONS AND THEIR APPROXIMATION The incompressible Navier-Stokes equations, in terms of the velocity vector u and the pressure p are given by V2u - Reu. Vu - Vp = 0 V.u
= 0.
(l.a) (l.b)
The Reynolds number Re is based on a characteristic velocity (which is also used to make the velocity vector dimensionless) a global characteristic length and the kinematic viscosity coefficient. To solve the system (1) one has to specify boundary condition on the components of the velocity vector. In our computations here, we have specified the velocity vector itself. In such cases one has to make sure that the total mass flux into the domain is zero. The system of differential equations is discretized on a rectangular staggered grid as shown in Fig. 1. The governing equations are satisfied inside or on the boundaries of each
Local mesh-refinement technique for incompressible flows
71
Equation
Uarlable w
x
U
x A
A
U
x
P
x
×
×
x
x
x
x-momentum y-momentum contlnulty
Fig. i. The location of the dependent variables on the staggered grid.
computational cell. Equations (1) are approximated by the following finite differences [18, 19]: V~Ui, j -- R e u . ~Ui, j -- Gpi.j = 0
(2.a)
D . u~,j = 0
(2.b)
where V2, G and D denote the central difference approximation to the Laplacian, the gradient and the divergence operators in (1), respectively. The non-linear terms of u. Vu are approximated by upwind differencing: e.g. uu:, is either 6uij(uij-
(2.c)
ui-~,j)/h
or
~ui, j(3ui, j - 4ui-~d + ui-26j)/2h
(2.d)
where ~ = sign(u/j). The formal accuracy of system (2) is determined by the finite difference approximation to the convective terms [(2.c) or (2.d)]. When the boundary conditions are applied as shown in Fig. l, the error on the velocity component parallel to the boundary is O(h), whereas the formal truncation error of the corresponding difference equations (2.a), at cells near the boundary, is O(h-I). To improve accuracy we introduce near the boundaries modified (first order accurate) schemes for the first and the second derivatives. Using the notation of Fig. 2, the normal second derivative is given by: 4
V ~ = (2Vb + II3 -- 3II2) 3h 2 .
(3)
The other normal derivatives at cells adjacent to the boundaries are approximated by corresponding formulae. The system of finite difference equations has been solved by using the MG procedure described in previous reports [ 18-21 ]. In the following we shall discuss only those aspects of the scheme which are related to the local mesh refinement technique.
ut
½
½
u4
Fig. 2. The finite-differences near the boundary (eqn (3)).
72
L. FucHs
Fig. 3. The streamline pattern in the lid-drivencavity,using local mesh refinements.R e
= 100.
LOCAL MESH REFINEMENTS AND MG-RELAXATIONS The M G relaxation procedure includes two steps of inter-grid information exchange. Fine grid approximations are used to correct (in a defect correction like manner) the coarse grid equations for their larger truncation errors. Coarse grid approximations correct the fine grid solution for global effects, such as boundary condition effects and slowly varying components of the solution. The former step gives automatically an estimate of the local truncation errors. This estimate cannot, however, be used directly as an adaptive grid refinement criterion. The truncation error estimates have to be normalized (by some constant) in such a way that when these errors are less than some given number, the corresponding solution
Fig. 4(a). T h e truncation error levels o n the u n i f o r m grid (l/h = 48) using (2.c). Re = 100. [] = +_1; • = + 1 0 ; × = __.100.
Local mesh-refinement technique for incompressible flows
73
Fig. 4(b). The truncation error levels on a refined sub-grid (I/h = 96) using (2.c). R e = 100. [] = +_I; •
= +10;
x = _+100.
attains a prescribed level of accuracy. Such an adaptive normalization shall be studied in the future. Here, we assume that the residuals are normalized such that the truncation errors should be less than unity everywhere, except possibly at some isolated computational cells. After identifying the subdomains where the truncation errors are too large one can 'cover' these sub-regions by locally refined grids. These new sub-grids are also rectangular with a mesh spacing equal to half of the mesh size in the coarser grid. All newly introduced sub-grids belong to the same level of mesh fineness. If additional levels of refinements are required, the new sub-grids must be contained completely in the corresponding sub-grid of •the coarser level.
Fig. 4(c). The truncation error levels on a refined sub-grid ( I / h = 192) using (2.c). R e = I00. [] = ±1; • = -+10; × = ±100.
74
L. FLICHS
Fig. 4(d). The truncation error levels on two refined sub-grids ( l / h = 384) using (2.c). R e = I00. [] = ±1; • = +10; X = _+100.
By the continuity equation (l.b) the total mass fluxes on each sub-grid must vanish identically. According to eqn (2.b) the mass fluxes are computed by the mid-point rule. The consistent interpolation formula (for the normal fluxes) from coarse to fine grids is by linear interpolations. The tangential velocity component, on the other hand, can be computed by other interpolation formulae without spoiling the total mass conservation. The effects of some different interpolation formulae for the velocity components shall be reported elsewhere. In the current computations we have used linear interpolations even for the tangential velocity component. It should be noted that no boundary conditions has to be, or can be,
II
Fig. 4(e). The truncation error levels on the uniform grid (1/h = 128) using (2.c). R e = 100. [] = _+1; •
= ±10; x = ±100.
Local mesh-refinement technique for incompressible flows
75
Fig. 5(a). The truncation error levels on the uniform grid (l/h = 48) using (2.d). Re = 100. [] = _+1; •
=
_+10; ×
=
+_100.
specified on the pressure. The form of transferring the pressure from the coarse to the fine grids does not effect the converged solution. ORGANISATION
OF
THE
GRID-SYSTEM
The grid hierarchy is defined in the following way: Each subgrid is identified by two integers; the level number (which indicates the mesh fineness) and the subgrid number (which indicates the chronological ordering of the subgrids belonging to the same level). The same subgrid can be identified even by a single integer, namely the grid-index. The
Fig. 5(b). The truncation error levels on a refined sub-grid ( l / h = 96) using (2.d). Re = 100. [] ~ ± 1; • = +-I0; x = :1:100.
76
L. FUCHS
Fig. 5(c). The t r u n c a t i o n error levels on a refined sub-grid ( l / h = 192) using (2.d). R e = 100. [] = - 1 ; • = - 1 0 ; X = +100.
later integer corresponds uniquely to the former pair of integers. For each grid one has to store the following data: (i) The coordinates of the origin of the local mesh and the indices of the corresponding mesh point on the next coarse subgrid. (ii) The subgrid-index of the next coarser grid. (iii) The number of computational cells in all space directions. The dependent variables (on all subgrids) are stored in vectors. Neighboring points in the
I
------.-----
•
Fig. 5(d). The t r u n c a t i o n error levels on the u n i f o r m grid (1/h = 128) using (2.d). R e = 100. [] = _+1; • = +10; x = +100.
Local mesh-refinementtechnique for incompressibleflows
77
Table 1. The maximal values of the truncation error (2.a) on different grids. l/h
eqn (2.c)
eqn (2.d)
48 64 96 128 192 384
690 870 2600 3300 8740 28300
610 840 2250 3260 7280
y-direction are also neighbors in the stored vector. A simple pointer system is used to identify each computational cell in each subgrid. For each subgrid-index a pointer indicates the first cell of each mesh line, x = constant. From this pointer vector, each cell location in the stored vector can be computed by a single addition. This system requires the storage of only 6 integers for each subgrid and an integer vector pointer which should be as long as the largest n u m b e r of cells in the x-direction in any of the subgrids. Since we use only cartesian coordinates the only information which has to be stored, besides the dependent variables, is the correction terms (for eqns. (2.a), but not for (2.b)) on each cell belonging to all levels but the finest one. The total m e m o r y requirement is m u c h smaller than in the corresponding case using non-uniform grids (even when these are optimally distributed). Transformation of coordinates may be necessary on some of the subgrids when the geometry cannot be approximated with an adequate accuracy on cartesian grids.
Fig. 6. The streamline pattern in the separator model, using local mesh refinements. Re = 50.
78
L. F o c r l s
ii
•
!t (a)
(b)
(c)
(d)
Fig. 7(a). T h e t r u n c a t i o n e r r o r levels o n the u n i f o r m g r i d (I/h = 24). Re = 50. [ ] = ___1; • = _+10; × = _+100. (b). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d (l/h = 48). Re = 50. [ ] = -+1; • = + 1 0 ; x = ± 1 0 0 . (c). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d (l/h = 96). Re = 50. [ ] = _+1; • = __+10; × = _+100. (d). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d (l/h = 96) u s i n g (2.d). Re = 50. [ ] = _ 1 ; • = _ 1 0 ; x = _+100.
COMPUTED
EXAMPLES
A. Lid-driven cavity The computation of this problem became a standard for testing and comparing numerical methods. The geometry of the computational domain is the simplest possible but the flow field is not smooth near the comers of the driving wall. First, we compute the flow field using the first order accurate finite-difference approximation ( 1. c ) . The streamline pattern on a uniform mesh with (48 X 48) intervals (but which further contains three finer levels, at Re --- 100), is shown in Fig. 3. The corresponding truncation error field is shown in Fig. 4(a). As the mesh is refined near the moving wall, regions with large truncation errors decrease in their size (Fig. 4(a)--4(d)). The last figure shows two subregions belonging to the same level. The total number of levels in this computation is 8 (including 9 subgrids). On the finest mesh h = 1/384. The truncation errors on a uniformly fine grid (h = 1/128) are shown in Fig. 4(e). Comparing this figure with Figs. 4(b) and 4(c), it is seen that the local meshes reduce the truncation errors as effectively as globally refined grids.
Local mesh-refinementtechnique for incompressibleflows
79
Fig. 8. The streamline pattern in the separator model, using local mesh refinements.R e = 100.
The results corresponding to the second order finite-differences (2.d) reveal an expected faster reduction in the truncation errors as the mesh size decreases locally (Figs. 5(a)-5(c)) or globally (Fig. 5(d)). As expected higher order schemes have the advantage that smaller subregions may be used to obtain the required level of truncation errors. It is interesting to note that for this problem the extreme values of truncation errors do not converge to zero as the mesh size is reduced. Table 1 displays the maximal value of the truncation errors for different cases. The truncation errors (at cells closest to the comers of the lid) increase at a rate close to O(h-2). This type of behavior can be explained by Batchelor's theory [22] for Stokes flow in a c o m e r between a moving and a stationary wall. By that theory, the pressure gradients at a distance r from the comer, behave as r 2. The truncation error terms, using central finite-differences are then O ( h 2 / r 4 ) . Thus, for the cells near the comers for which r ~ h the truncation errors are O(h-2). The Stokes approximation is valid only sufficiently close to the comer i.e. where r . R e ~ I. In our case this condition is not satisfied. Nevertheless, the Stokes approximation explains, qualitatively, the comer singularity of the truncation errors. This truncation error singularity decays as f 4 (where r is the distance from the comer). Thus, for any fixed point the level of the truncation errors decrease as the mesh is refined. On the other hand, the level of the truncation errors at the computational cell placed at the singular comers of the cavity, increase with mesh refinement, since the centers of the cells are moved closer to the comers. However, this singular behavior of the solution (and that of the truncation errors) has only a local effect. Detailed studies of computed cavity flows [23] show that the numerical solutions converge pointwise as the mesh is refined, in all subdomains excluding the comers.
80
L. FUCHS
B. Separator model In this case, fluid enters the domain (almost) parallel to the x-axis and leaves the region parallel to the y-axis. Particles with different masses are separated by the centrifugal forces. The computed streamline pattern (Re = 50) is shown in Fig. 6. Figs. 7(a)-7(c) show the truncation error levels as new levels are introduced. The corresponding cases for Re = 100 are shown in Figs. 8 and 9. For this problem, local mesh refinements lead to a uniform reduction of the truncation errors. The truncation error for the second order scheme with local mesh refinements (corresponding to Fig. 7(c) is shown in Fig. 7(d). The superiority of the second order scheme over the first order, is evident. For Re = 100, we note an increase in the truncation error level near the outflow boundary. This is a result o f imposing free-stream conditions at a distance too close to the entrance. When the length of the channel is increased the level of the 'truncation error' decreases. This case is an example of other useful information which can be obtained in the adaptive scheme: Large truncation errors near the boundaries may imply either singularity in the solution or errors in the specified boundary conditions. These types of "truncation errors" cannot be reduced by mesh refinements.
{'\
(a)
(b)
(c)
Fig. 9(a). T h e t r u n c a t i o n e r r o r levels o n the u n i f o r m g r i d ( l / h = 24). R e = 100. [ ] = +_1; • = _+!0; × = _+100. (b). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d ( l / h = 48). R e = 100. [ ] = __.l; • = _+10; × = _+100. (c). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d ( l / h = 96). R e = 100. [ ] = _+1; • = + 1 0 ; × = _ 1 0 0 .
Local mesh-refinement technique for incompressible flows
81
CONCLUDING REMARKS T h e m e t h o d o f d y n a m i c a l l y defining local subgrids d u r i n g the process o f M G - r e l a x a t i o n s was studied. New, locally-refined grids c a n be i n t r o d u c e d in a r a t h e r s i m p l e m a n n e r a n d t h e convergence rate o f t h e M G s c h e m e is as g o o d as with global grids. T h e g r i d - p o i n t e r system w h i c h has been o u t l i n e d above, requires very little a d d i t i o n a l c o m p u t e r m e m o r y . T h e local subgrids c a n reduce the t r u n c a t i o n errors in well identified regions. F u r t h e r m o r e , b y s t u d y i n g the t r u n c a t i o n errors (which are c o m p u t e d a n y w a y d u r i n g the M G - c y c l i n g ) o n e m a y even get i n f o r m a t i o n a b o u t inaccuracies in a p p l y i n g the b o u n d a r y c o n d i t i o n s . T h e i n c o m p r e s s i b l e N a v i e r - S t o k e s equations, in p r i m i t i v e variables form, have been solved b y the s c h e m e for t w o cases. F o r the lid-driven cavity p r o b l e m , local grids were i n t r o d u c e d near the driving wall. T h e t r u n c a t i o n errors were reduced except near the singular corners, where they increased as finer meshes were introduced. T h e behavior o f the t r u n c a t i o n errors n e a r the corners is in g o o d a g r e e m e n t with t h a t p r e d i c a t e d by Stokes equations. In the case o f the s e p a r a t o r m o d e l the s c h e m e c o u l d reduce substantially, the n u m b e r o f c o m p u t a t i o n a l n o d e p o i n t s while the d i s t r i b u t i o n o f the t r u n c a t i o n errors b e c a m e m o r e u n i f o r m . REFERENCES I. J. F. Thompson and Z. U. A. Warsi, Boundary fitted coordinate systems for numerical solution of partial differential equations. J. Comp. Phys. 47, l (1982). 2. A. S. Benjamin and V. E. Denny, On the convergence of numerical solutions for 2-D flows in a cavity at large Re. J. Comp. Phys. 33, 340 (1979). 3. C. F. Shieh, "Three-Dimensioual Grid Generation Using Elliptic Equations With Direct Grid Distribution Control." AIAA-83-448 (1983). 4. C. W. Mastin and J. F. Thompson, "Adaptive Grids Generated by Elliptic Systems." AIAA-83-451 (1983). 5. M. M. Rai, "A Conservative Treatment of Zonal Boundaries for Euler Equation Calculations." AIAA-84-164 (1984). 6. M. M. Rai and D. A. Anderson, Grid evolution in time asymptotic problems. J. Comp. Phys. 43, 327 (1981). 7. D. S. Thompson, A Mesh Embedded Approachfor Prediction of Transonic wing~body~storeFlow Fields. Numerical Boundary Condition Procedures Symposium, NASA CP-2201, 403 (1981). 8. S. G. Hedman, "An Application of Embedded Grid Technique to the Calculation of Transonic Flow Past Wings." FFA, Technical Note AU-1600, Stockholm (1980). 9. J. A. Benek, J. L. Steger and F. C. Dougherty, "A Flexible Grid Embedding Technique With Application to the Euler Equations." AIAA-83-1944 (1983). 10. V. Pereyra and E. G. SCweU,Mesh selection for discrete solution of boundary problems in ordinary differential equations. Numer. Math. 23, 261 (1975). 11. T. H. Chong, A variable mesh finite-difference method for solving a class of parabolic differential equations in one space variable. SIAM J. Numer. Anal. 15, 835 (1978). 12. S. F. Davis and J. E. Flaherty, An adaptive finite-difference method for initial-boundary value problems for partial differential equations. SIAMJ. Sci. Stat. Comput. 3, 6 (1982). 13. H. A. Dwyer, R. J. Kee and B. R. Sanders, Adaptive grid method for problems in fluid mechanics and heat transfer. AI.4A J. 18, 1205 (1980). 14. M. J. Berger and J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations. J. Comp. Phys. 53, 484 (1984). 15. A. Brandt, Multigrid Techniques: 1984 Guide. Computational Fluid Dynamics Lecture Series at the vonKarman Institute, March, 1984. 16. W. J. Usab and E. M. Murman, "Embedded Mesh Solutions of the Euler Equations Using a Multiple-Grid Method." AIAA-83-1946, (1983). 17. J. J. Brown, An embedded potential flow analysis. AIAA J. 22, 174 (1983). 18. T. Thunell and L. Fuchs, Numerical solution of the Navier-Stokes equations by multi-grid techniques, in Numerical Methods in Laminar and Turbulent Flow-H (Edited by C. Taylor and B. A. Schrefler). Pineridge Press, 141 (1981). 19. L. Fuchs, Multi-grid schemes for incompressible flows. Proc. of the GAMM-workshop on 'Et~cient Solvers for Elliptic Systems' (Edited by W. Hackbusch). Notes on Numerical Fluid Mechanics, 10, 38 (1984). 20. L. Fuchs, "Multi-Grid Solution of the Navier-Stokes Equations on Non-Uniform Grids. Multigrid Methods." NASA-CP 2202, 83 (1981). 21. L. Fuchs and H-S. Zhao, Solution of three-dimensional viscous incompressible flows by a multi-grid method. Int. J. Numerical Methods in Fluids. 4, 539 (1984). 22. G. K. Batchelor, An Introduction to Fluid Dynamics. Cambridge University Press (1967). 23. L. Fuchs, Higher Order Finite-Difference Schemes for Driven Cavity Flows. to appear (1985).
0045-7930/86 $3.00 + .IX) © 1986 Pergamon Press Ltd.
Printed in Great Bdtain.
A LOCAL MESH-REFINEMENT TECHNIQUE FOR INCOMPRESSIBLE FLOWS L. FUCHS Delmrtment of Gasdyuamics, The Royal Institute of Technology, Stockholm, Sweden
(Received 26 November 1984; in revisedform 8 February 1985) Abstract--As part of a Multi-Grid scheme for the solution of the Navier-Stokes equations in primitive variables, we introduce a local mesh refinement procedure. New cartesian sub-grids are introduced into regions where the estimated truncation errors are too large. Through the Multi-Grid processing, informations is transferred among the grids in a stable and efficient manner. A simple pointer system allows the storage of the dependent variables, without increasing in the required computer memory. Two computed examples of incompressible flow problems are discussed.
INTRODUCTION
One of the complicating factors in the solution of the Navier-Stokes equations is due to the existence of different (length) scales in the flow field. The dependent variables are made dimensionless by some (globally valid) characteristic scale. This scale may not be represen= tative in all parts of the flow field and locally one may have much smaller or larger characteristic scales. The accuracy of finite-differences (or any other discrete approximation) depends not only on the formal order of approximation but also on the ratio of the mesh size (h) and the local scale (L). For accuracy, one has to resolve the local scales so that the variations in the dependent variables are described on as many mesh intervals as possible (h/L = o(1)). Generally, it is impractical to use uniformly fine meshes because of computer memory restrictions and long computational times. Furthermore, such fine meshes are not needed for accuracy, in larger parts of the flow field. A frequently used method is to construct a non-uniform (body fitted) grid so that the mesh spacing is small in regions with expectedly large gradients (e.g. [1, 2]). The main disadvantage of this method is that one has to know, a priori, the details of the solution when the mesh is generated. Furthermore, the transformation of the coordinates requires additional storage (of at least the location of the node points) and additional computational operations (on the additional terms in the governing equations). Alternative approaches to global transformation of coordinates are grid patching, mesh embedding and adaptive grid generation using either orthogonal or non-orthogonal grid systems. The method of grid patching uses the idea of dividing the domain into sub-domains which accommodate independently generated meshes. These grids are joined together along the boundaries of the sub-domains. This method has been used for the computation of some transonic flows [5, 6]. The approach requires careful treatment of the information transfer across the boundaries of the sub-domains. Mesh embedding can be considered as a special case of patching with constraints imposed on the grids on both sides of the 'internal' (sub-domain) boundaries. Embedded meshes have, usually, the same basic structure as the rest of the mesh in which they are being embedded. This method has also been applied to transonic problems with success [7, 8]. A more generalized mesh embedding/patching has been proposed recently by Benek, Steger and Dougherty [9]. Both in mesh patching and mesh embedding techniques the only interaction among the different sub-domains (sub-grids) is done by updating (transferring) boundary data. Such procedures reduce in general, the efficiency and in some cases also the stability properties of the numerical scheme. Adaptive schemes are the most desirable for obtaining top efficiency in terms of accuracy. On the other hand such optimal grids requires by themselves considerable amount of corn69 CAF 14:1-F
70
L. FUCHS
putational effort. For I-D problems (or even for time-dependent 1-D cases) several methods have been developed [10-12]. An adaptive grid generation procedure for multidimensional problems has been developed and described by Dwyer, Kee and Sanders [ 13]. Their adaptive scheme is based upon the placement of grid points in proportion to the gradients in the dependent variables. Here, we do not attempt to produce 'optimal' grids. Instead, we use a system of cartesian grids, arranged in a grid hierarchy. New sub-grids are introduced, at places where the estimated truncation error is too large. In this scheme, the new local grids form new levels in the hierarchy of grids. This approach is close in spirit to the approach of Berger and Oliger [14]. However, by realising that adaptive grid generation procedures and the solution algorithms are strongly connected, we construct the grid-system dynamically during a MultiGrid (MG) solution procedure. In fact, several elements of the current scheme have been proposed by Brandt [ 15]. Recently, the combination of mesh embedding with a MG algorithm have been reported by Usab and Murman [16] and by Brown [17]. In both papers the emphasis is on the gains in computational time compared to global mesh refinement computations, and no attempt is made to control the location or the size of the locally refined grids. In the present work we study the use of locally refined meshes as an integral part of a MG scheme for the solution of incompressible flow problems. As a further step toward a fully adaptive grid generation, we study the behavior of the truncation errors for different levels of refinement. The values of the appropriate dependent variables which must be specified on the boundaries of the sub-grids are computed by interpolations from the coarser grids. These interpolation formulae must, however, be such that the mass fluxes are conserved. In non-MG schemes the interaction among the patched/embedded grids takes place only by the updating procedure of the (internal) boundary conditions (see e.g. [5, 9, 14]). In full MG relaxations, information is passed from the (global) coarser grids to the locally refined grids by the interpolated corrections and the 'internal boundary conditions'. Information, which gives rise to improved global accuracy, is passed back from the fine local grid to the coarser grids by adding a truncation error related term to the equations. The scheme can identify regions with large truncation errors and by using local rectangular meshes, these errors can be reduced to a desirable (uniform) level. The administration of the local grids is very simple and requires minimal 'book-keeping'. In contrast to non-uniform, 'optimal' grids, one has to store only six integer values for each sub-grid. In the following we describe some details of the numerical scheme in the MG context. The method has been applied to the computation of the flow in the square driven cavity and a simplified separator model. By using very fine local meshes near the (singular) comers of the driven cavity, it is found that the truncation error has a singular behavior (O(h-2)) which is in agreement with the Stokes flow in such regions. In both cases, however, the level of the truncation errors can be controlled very efficiently by using locally refined grids. GOVERNING EQUATIONS AND THEIR APPROXIMATION The incompressible Navier-Stokes equations, in terms of the velocity vector u and the pressure p are given by V2u - Reu. Vu - Vp = 0 V.u
= 0.
(l.a) (l.b)
The Reynolds number Re is based on a characteristic velocity (which is also used to make the velocity vector dimensionless) a global characteristic length and the kinematic viscosity coefficient. To solve the system (1) one has to specify boundary condition on the components of the velocity vector. In our computations here, we have specified the velocity vector itself. In such cases one has to make sure that the total mass flux into the domain is zero. The system of differential equations is discretized on a rectangular staggered grid as shown in Fig. 1. The governing equations are satisfied inside or on the boundaries of each
Local mesh-refinement technique for incompressible flows
71
Equation
Uarlable w
x
U
x A
A
U
x
P
x
×
×
x
x
x
x-momentum y-momentum contlnulty
Fig. i. The location of the dependent variables on the staggered grid.
computational cell. Equations (1) are approximated by the following finite differences [18, 19]: V~Ui, j -- R e u . ~Ui, j -- Gpi.j = 0
(2.a)
D . u~,j = 0
(2.b)
where V2, G and D denote the central difference approximation to the Laplacian, the gradient and the divergence operators in (1), respectively. The non-linear terms of u. Vu are approximated by upwind differencing: e.g. uu:, is either 6uij(uij-
(2.c)
ui-~,j)/h
or
~ui, j(3ui, j - 4ui-~d + ui-26j)/2h
(2.d)
where ~ = sign(u/j). The formal accuracy of system (2) is determined by the finite difference approximation to the convective terms [(2.c) or (2.d)]. When the boundary conditions are applied as shown in Fig. l, the error on the velocity component parallel to the boundary is O(h), whereas the formal truncation error of the corresponding difference equations (2.a), at cells near the boundary, is O(h-I). To improve accuracy we introduce near the boundaries modified (first order accurate) schemes for the first and the second derivatives. Using the notation of Fig. 2, the normal second derivative is given by: 4
V ~ = (2Vb + II3 -- 3II2) 3h 2 .
(3)
The other normal derivatives at cells adjacent to the boundaries are approximated by corresponding formulae. The system of finite difference equations has been solved by using the MG procedure described in previous reports [ 18-21 ]. In the following we shall discuss only those aspects of the scheme which are related to the local mesh refinement technique.
ut
½
½
u4
Fig. 2. The finite-differences near the boundary (eqn (3)).
72
L. FucHs
Fig. 3. The streamline pattern in the lid-drivencavity,using local mesh refinements.R e
= 100.
LOCAL MESH REFINEMENTS AND MG-RELAXATIONS The M G relaxation procedure includes two steps of inter-grid information exchange. Fine grid approximations are used to correct (in a defect correction like manner) the coarse grid equations for their larger truncation errors. Coarse grid approximations correct the fine grid solution for global effects, such as boundary condition effects and slowly varying components of the solution. The former step gives automatically an estimate of the local truncation errors. This estimate cannot, however, be used directly as an adaptive grid refinement criterion. The truncation error estimates have to be normalized (by some constant) in such a way that when these errors are less than some given number, the corresponding solution
Fig. 4(a). T h e truncation error levels o n the u n i f o r m grid (l/h = 48) using (2.c). Re = 100. [] = +_1; • = + 1 0 ; × = __.100.
Local mesh-refinement technique for incompressible flows
73
Fig. 4(b). The truncation error levels on a refined sub-grid (I/h = 96) using (2.c). R e = 100. [] = +_I; •
= +10;
x = _+100.
attains a prescribed level of accuracy. Such an adaptive normalization shall be studied in the future. Here, we assume that the residuals are normalized such that the truncation errors should be less than unity everywhere, except possibly at some isolated computational cells. After identifying the subdomains where the truncation errors are too large one can 'cover' these sub-regions by locally refined grids. These new sub-grids are also rectangular with a mesh spacing equal to half of the mesh size in the coarser grid. All newly introduced sub-grids belong to the same level of mesh fineness. If additional levels of refinements are required, the new sub-grids must be contained completely in the corresponding sub-grid of •the coarser level.
Fig. 4(c). The truncation error levels on a refined sub-grid ( I / h = 192) using (2.c). R e = I00. [] = ±1; • = -+10; × = ±100.
74
L. FLICHS
Fig. 4(d). The truncation error levels on two refined sub-grids ( l / h = 384) using (2.c). R e = I00. [] = ±1; • = +10; X = _+100.
By the continuity equation (l.b) the total mass fluxes on each sub-grid must vanish identically. According to eqn (2.b) the mass fluxes are computed by the mid-point rule. The consistent interpolation formula (for the normal fluxes) from coarse to fine grids is by linear interpolations. The tangential velocity component, on the other hand, can be computed by other interpolation formulae without spoiling the total mass conservation. The effects of some different interpolation formulae for the velocity components shall be reported elsewhere. In the current computations we have used linear interpolations even for the tangential velocity component. It should be noted that no boundary conditions has to be, or can be,
II
Fig. 4(e). The truncation error levels on the uniform grid (1/h = 128) using (2.c). R e = 100. [] = _+1; •
= ±10; x = ±100.
Local mesh-refinement technique for incompressible flows
75
Fig. 5(a). The truncation error levels on the uniform grid (l/h = 48) using (2.d). Re = 100. [] = _+1; •
=
_+10; ×
=
+_100.
specified on the pressure. The form of transferring the pressure from the coarse to the fine grids does not effect the converged solution. ORGANISATION
OF
THE
GRID-SYSTEM
The grid hierarchy is defined in the following way: Each subgrid is identified by two integers; the level number (which indicates the mesh fineness) and the subgrid number (which indicates the chronological ordering of the subgrids belonging to the same level). The same subgrid can be identified even by a single integer, namely the grid-index. The
Fig. 5(b). The truncation error levels on a refined sub-grid ( l / h = 96) using (2.d). Re = 100. [] ~ ± 1; • = +-I0; x = :1:100.
76
L. FUCHS
Fig. 5(c). The t r u n c a t i o n error levels on a refined sub-grid ( l / h = 192) using (2.d). R e = 100. [] = - 1 ; • = - 1 0 ; X = +100.
later integer corresponds uniquely to the former pair of integers. For each grid one has to store the following data: (i) The coordinates of the origin of the local mesh and the indices of the corresponding mesh point on the next coarse subgrid. (ii) The subgrid-index of the next coarser grid. (iii) The number of computational cells in all space directions. The dependent variables (on all subgrids) are stored in vectors. Neighboring points in the
I
------.-----
•
Fig. 5(d). The t r u n c a t i o n error levels on the u n i f o r m grid (1/h = 128) using (2.d). R e = 100. [] = _+1; • = +10; x = +100.
Local mesh-refinementtechnique for incompressibleflows
77
Table 1. The maximal values of the truncation error (2.a) on different grids. l/h
eqn (2.c)
eqn (2.d)
48 64 96 128 192 384
690 870 2600 3300 8740 28300
610 840 2250 3260 7280
y-direction are also neighbors in the stored vector. A simple pointer system is used to identify each computational cell in each subgrid. For each subgrid-index a pointer indicates the first cell of each mesh line, x = constant. From this pointer vector, each cell location in the stored vector can be computed by a single addition. This system requires the storage of only 6 integers for each subgrid and an integer vector pointer which should be as long as the largest n u m b e r of cells in the x-direction in any of the subgrids. Since we use only cartesian coordinates the only information which has to be stored, besides the dependent variables, is the correction terms (for eqns. (2.a), but not for (2.b)) on each cell belonging to all levels but the finest one. The total m e m o r y requirement is m u c h smaller than in the corresponding case using non-uniform grids (even when these are optimally distributed). Transformation of coordinates may be necessary on some of the subgrids when the geometry cannot be approximated with an adequate accuracy on cartesian grids.
Fig. 6. The streamline pattern in the separator model, using local mesh refinements. Re = 50.
78
L. F o c r l s
ii
•
!t (a)
(b)
(c)
(d)
Fig. 7(a). T h e t r u n c a t i o n e r r o r levels o n the u n i f o r m g r i d (I/h = 24). Re = 50. [ ] = ___1; • = _+10; × = _+100. (b). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d (l/h = 48). Re = 50. [ ] = -+1; • = + 1 0 ; x = ± 1 0 0 . (c). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d (l/h = 96). Re = 50. [ ] = _+1; • = __+10; × = _+100. (d). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d (l/h = 96) u s i n g (2.d). Re = 50. [ ] = _ 1 ; • = _ 1 0 ; x = _+100.
COMPUTED
EXAMPLES
A. Lid-driven cavity The computation of this problem became a standard for testing and comparing numerical methods. The geometry of the computational domain is the simplest possible but the flow field is not smooth near the comers of the driving wall. First, we compute the flow field using the first order accurate finite-difference approximation ( 1. c ) . The streamline pattern on a uniform mesh with (48 X 48) intervals (but which further contains three finer levels, at Re --- 100), is shown in Fig. 3. The corresponding truncation error field is shown in Fig. 4(a). As the mesh is refined near the moving wall, regions with large truncation errors decrease in their size (Fig. 4(a)--4(d)). The last figure shows two subregions belonging to the same level. The total number of levels in this computation is 8 (including 9 subgrids). On the finest mesh h = 1/384. The truncation errors on a uniformly fine grid (h = 1/128) are shown in Fig. 4(e). Comparing this figure with Figs. 4(b) and 4(c), it is seen that the local meshes reduce the truncation errors as effectively as globally refined grids.
Local mesh-refinementtechnique for incompressibleflows
79
Fig. 8. The streamline pattern in the separator model, using local mesh refinements.R e = 100.
The results corresponding to the second order finite-differences (2.d) reveal an expected faster reduction in the truncation errors as the mesh size decreases locally (Figs. 5(a)-5(c)) or globally (Fig. 5(d)). As expected higher order schemes have the advantage that smaller subregions may be used to obtain the required level of truncation errors. It is interesting to note that for this problem the extreme values of truncation errors do not converge to zero as the mesh size is reduced. Table 1 displays the maximal value of the truncation errors for different cases. The truncation errors (at cells closest to the comers of the lid) increase at a rate close to O(h-2). This type of behavior can be explained by Batchelor's theory [22] for Stokes flow in a c o m e r between a moving and a stationary wall. By that theory, the pressure gradients at a distance r from the comer, behave as r 2. The truncation error terms, using central finite-differences are then O ( h 2 / r 4 ) . Thus, for the cells near the comers for which r ~ h the truncation errors are O(h-2). The Stokes approximation is valid only sufficiently close to the comer i.e. where r . R e ~ I. In our case this condition is not satisfied. Nevertheless, the Stokes approximation explains, qualitatively, the comer singularity of the truncation errors. This truncation error singularity decays as f 4 (where r is the distance from the comer). Thus, for any fixed point the level of the truncation errors decrease as the mesh is refined. On the other hand, the level of the truncation errors at the computational cell placed at the singular comers of the cavity, increase with mesh refinement, since the centers of the cells are moved closer to the comers. However, this singular behavior of the solution (and that of the truncation errors) has only a local effect. Detailed studies of computed cavity flows [23] show that the numerical solutions converge pointwise as the mesh is refined, in all subdomains excluding the comers.
80
L. FUCHS
B. Separator model In this case, fluid enters the domain (almost) parallel to the x-axis and leaves the region parallel to the y-axis. Particles with different masses are separated by the centrifugal forces. The computed streamline pattern (Re = 50) is shown in Fig. 6. Figs. 7(a)-7(c) show the truncation error levels as new levels are introduced. The corresponding cases for Re = 100 are shown in Figs. 8 and 9. For this problem, local mesh refinements lead to a uniform reduction of the truncation errors. The truncation error for the second order scheme with local mesh refinements (corresponding to Fig. 7(c) is shown in Fig. 7(d). The superiority of the second order scheme over the first order, is evident. For Re = 100, we note an increase in the truncation error level near the outflow boundary. This is a result o f imposing free-stream conditions at a distance too close to the entrance. When the length of the channel is increased the level of the 'truncation error' decreases. This case is an example of other useful information which can be obtained in the adaptive scheme: Large truncation errors near the boundaries may imply either singularity in the solution or errors in the specified boundary conditions. These types of "truncation errors" cannot be reduced by mesh refinements.
{'\
(a)
(b)
(c)
Fig. 9(a). T h e t r u n c a t i o n e r r o r levels o n the u n i f o r m g r i d ( l / h = 24). R e = 100. [ ] = +_1; • = _+!0; × = _+100. (b). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d ( l / h = 48). R e = 100. [ ] = __.l; • = _+10; × = _+100. (c). T h e t r u n c a t i o n e r r o r levels o n a refined s u b - g r i d ( l / h = 96). R e = 100. [ ] = _+1; • = + 1 0 ; × = _ 1 0 0 .
Local mesh-refinement technique for incompressible flows
81
CONCLUDING REMARKS T h e m e t h o d o f d y n a m i c a l l y defining local subgrids d u r i n g the process o f M G - r e l a x a t i o n s was studied. New, locally-refined grids c a n be i n t r o d u c e d in a r a t h e r s i m p l e m a n n e r a n d t h e convergence rate o f t h e M G s c h e m e is as g o o d as with global grids. T h e g r i d - p o i n t e r system w h i c h has been o u t l i n e d above, requires very little a d d i t i o n a l c o m p u t e r m e m o r y . T h e local subgrids c a n reduce the t r u n c a t i o n errors in well identified regions. F u r t h e r m o r e , b y s t u d y i n g the t r u n c a t i o n errors (which are c o m p u t e d a n y w a y d u r i n g the M G - c y c l i n g ) o n e m a y even get i n f o r m a t i o n a b o u t inaccuracies in a p p l y i n g the b o u n d a r y c o n d i t i o n s . T h e i n c o m p r e s s i b l e N a v i e r - S t o k e s equations, in p r i m i t i v e variables form, have been solved b y the s c h e m e for t w o cases. F o r the lid-driven cavity p r o b l e m , local grids were i n t r o d u c e d near the driving wall. T h e t r u n c a t i o n errors were reduced except near the singular corners, where they increased as finer meshes were introduced. T h e behavior o f the t r u n c a t i o n errors n e a r the corners is in g o o d a g r e e m e n t with t h a t p r e d i c a t e d by Stokes equations. In the case o f the s e p a r a t o r m o d e l the s c h e m e c o u l d reduce substantially, the n u m b e r o f c o m p u t a t i o n a l n o d e p o i n t s while the d i s t r i b u t i o n o f the t r u n c a t i o n errors b e c a m e m o r e u n i f o r m . REFERENCES I. J. F. Thompson and Z. U. A. Warsi, Boundary fitted coordinate systems for numerical solution of partial differential equations. J. Comp. Phys. 47, l (1982). 2. A. S. Benjamin and V. E. Denny, On the convergence of numerical solutions for 2-D flows in a cavity at large Re. J. Comp. Phys. 33, 340 (1979). 3. C. F. Shieh, "Three-Dimensioual Grid Generation Using Elliptic Equations With Direct Grid Distribution Control." AIAA-83-448 (1983). 4. C. W. Mastin and J. F. Thompson, "Adaptive Grids Generated by Elliptic Systems." AIAA-83-451 (1983). 5. M. M. Rai, "A Conservative Treatment of Zonal Boundaries for Euler Equation Calculations." AIAA-84-164 (1984). 6. M. M. Rai and D. A. Anderson, Grid evolution in time asymptotic problems. J. Comp. Phys. 43, 327 (1981). 7. D. S. Thompson, A Mesh Embedded Approachfor Prediction of Transonic wing~body~storeFlow Fields. Numerical Boundary Condition Procedures Symposium, NASA CP-2201, 403 (1981). 8. S. G. Hedman, "An Application of Embedded Grid Technique to the Calculation of Transonic Flow Past Wings." FFA, Technical Note AU-1600, Stockholm (1980). 9. J. A. Benek, J. L. Steger and F. C. Dougherty, "A Flexible Grid Embedding Technique With Application to the Euler Equations." AIAA-83-1944 (1983). 10. V. Pereyra and E. G. SCweU,Mesh selection for discrete solution of boundary problems in ordinary differential equations. Numer. Math. 23, 261 (1975). 11. T. H. Chong, A variable mesh finite-difference method for solving a class of parabolic differential equations in one space variable. SIAM J. Numer. Anal. 15, 835 (1978). 12. S. F. Davis and J. E. Flaherty, An adaptive finite-difference method for initial-boundary value problems for partial differential equations. SIAMJ. Sci. Stat. Comput. 3, 6 (1982). 13. H. A. Dwyer, R. J. Kee and B. R. Sanders, Adaptive grid method for problems in fluid mechanics and heat transfer. AI.4A J. 18, 1205 (1980). 14. M. J. Berger and J. Oliger, Adaptive mesh refinement for hyperbolic partial differential equations. J. Comp. Phys. 53, 484 (1984). 15. A. Brandt, Multigrid Techniques: 1984 Guide. Computational Fluid Dynamics Lecture Series at the vonKarman Institute, March, 1984. 16. W. J. Usab and E. M. Murman, "Embedded Mesh Solutions of the Euler Equations Using a Multiple-Grid Method." AIAA-83-1946, (1983). 17. J. J. Brown, An embedded potential flow analysis. AIAA J. 22, 174 (1983). 18. T. Thunell and L. Fuchs, Numerical solution of the Navier-Stokes equations by multi-grid techniques, in Numerical Methods in Laminar and Turbulent Flow-H (Edited by C. Taylor and B. A. Schrefler). Pineridge Press, 141 (1981). 19. L. Fuchs, Multi-grid schemes for incompressible flows. Proc. of the GAMM-workshop on 'Et~cient Solvers for Elliptic Systems' (Edited by W. Hackbusch). Notes on Numerical Fluid Mechanics, 10, 38 (1984). 20. L. Fuchs, "Multi-Grid Solution of the Navier-Stokes Equations on Non-Uniform Grids. Multigrid Methods." NASA-CP 2202, 83 (1981). 21. L. Fuchs and H-S. Zhao, Solution of three-dimensional viscous incompressible flows by a multi-grid method. Int. J. Numerical Methods in Fluids. 4, 539 (1984). 22. G. K. Batchelor, An Introduction to Fluid Dynamics. Cambridge University Press (1967). 23. L. Fuchs, Higher Order Finite-Difference Schemes for Driven Cavity Flows. to appear (1985).