- Email: [email protected]
Adaptive Cartesian grid with mesh-less zones for compressible flow calculations
Computers & Fluids 54 (2012) 10–17
Contents lists available at SciVerse ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Adaptive Cartesian grid with mesh-less zones for compressible flow calculations q A. Jahangirian ⇑, M.Y. Hashemi Department of Aerospace Engineering, Amirkabir University of Technology, 424 Hafez Avenue, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 1 May 2010 Received in revised form 13 August 2011 Accepted 16 August 2011 Available online 6 September 2011 Keywords: Adaptive Cartesian grid Mesh-less method Compressible flow
a b s t r a c t A hybrid grid method that combines the computational efficiency of adaptive Cartesian grids with the flexibility of mesh-less scheme is developed for calculation of compressible flows. The mesh-less zone is created around the geometry by producing layers of nodes along normal direction vectors. Last layer is used as virtual geometry for automatic generation of unstructured Cartesian grid around mesh-less zone. An efficient explicit central difference scheme with artificial dissipation terms and convergence acceleration techniques is developed for both Cartesian grid and mesh-less zones. The method is used for computation of inviscid and viscous flows around single and multi-element airfoils at transonic flow conditions. Results indicate good agreements with those of available unstructured finite-volume flow solver (Uns2D) and other reference numerical data. The new method is shown to reduce the computational cost up to 70% compared with the fully unstructured flow solver. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction The main problem of computational fluid dynamics (CFD) for configurations with complex geometry is mesh generation. The structured grid methods are not efficient for such cases [1]. The main advantage of unstructured grid methods is the facility of grid generation for complex configurations [2]. However, the computational costs and memory requirements are generally higher than their structured grid counterparts. Alternatively, the Cartesian grid methods can be used that offer ease of grid generation, lower computational storage requirements, and significantly less operational count per cell compared to body fitted curvilinear grids [3]. The convergence properties of the solver are also potentially better because there are no problems related to skewness or distortion of cells and implementation of high order discretization schemes is simple. Embedded and adaptive grids, can also be used to provide better resolution of geometry and flow features in Cartesian grid [4]. However, the main challenge in using this method deals with arbitrary boundaries. As the grids are not body aligned, Cartesian cells near the body can extend through surfaces of solid components. Hence, accurate means of representations for surface boundary conditions are essential for the success of Cartesian schemes. Different methods have been proposed in the literature to resolve the boundary conditions problem in Cartesian grids. Broadly the methods proposed involve either cut cells for a finite volume treatment [5] or ghost points for a finite difference construction
[6]. A cut-cell method ensures conservation and uses the same flux computation algorithm as for the interior cells. However, the task of computing the volume and fluxes for all of the irregularly shaped cut cells and hanging nodes entail a considerable increase in complexity. More crucially, the method may lead to the creation of very small (irregular) cells at the boundary. Additionally, in the viscous flows the favorable normal lines to the surface and regular cells are not possible to generate which pose problems of numerical stability and computational cost. Recently, mesh-less methods have been developed to conveniently and accurately treats boundary conditions on curved boundaries based on a least-square approach [7]. However, the mesh-less method needs further development in automating the procedure for any geometry particularly in viscous applications. The main objective of the present work is to combine the advantages of recent works regarding the Cartesian grid method [4] and the mesh-less strategy [8] for computation of compressible flows. The method includes generating a layer of nodes around the body while Cartesian grid method is used elsewhere. Since the majority of computational domain is solved using the Cartesian grid method, the resulting procedure is very efficient in terms of both computational cost and storage requirements. In addition, since the mesh-less method is used to obtain the solution for the solid boundary points the developed method can easily be used for complex geometries without need to generate the body-fitted mesh. 2. Mesh-less zone point creation
q
This paper was originally submitted to the ICFD 2010 special issue which was published in Volume 46, Issue 1 of Computers and Fluids. ⇑ Corresponding author. E-mail address: [email protected] (A. Jahangirian). 0045-7930/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2011.08.010
In the first step, the surface must be broken up into elements (edges in 2D and faces in 3D) for a full description of a complex geometry. To generate the efficient distribution of points in the
11
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
Fig. 1. Mesh-less zone with its boundary for SKF1.1 two element airfoil.
mesh-less zone around the geometry the unit normal vector for all vertices on the geometry are calculated. At the next step, the layers of nodes are produced along the normal lines until a user specified numbers of layers (Fig. 1). For efficient layer generation the normal vectors can be smoothed where needed and/or truncated where overlapping the layers occurs. The neighbors of any node in the mesh-less zone need to be identified for the discretization of equations. Thus, a node-based data structure is formed which defines the neighbors of each node. In the present work, since the nodes inside the mesh-less zones are generated in a regular manner, the neighbors of each node can be easily defined. In the area that mesh-less zones overlap with each other, neighbors can be obtained by a simple local search. 3. Adaptive Cartesian grid generation Generation of the unstructured Cartesian grid is carried out following the work of Jahangirian and Shoraka [4] which is explained here briefly as follows: (1) Define an outer boundary box for the flow domain. (2) Define surface geometries (here the outer boundary of the mesh-less zone). (3) Subdivide the box and subsequent boxes if they contain or intersect the geometry until the required level of embedded grids is reached. (4) Detect and remove cells inside the mesh-less zone boundary. To complement the mesh-less zone for solution of the governing flow equations a series of halo points are added outside the mesh-less zone. After numerical experiments the first two layers of Cartesian cells are selected and their centers are considered for mesh-less zone computations. The close view of the generated final grid (including halo points) around a wing-flap configuration is shown in Fig. 2. One of the advantages of the Cartesian-based grid generation methods is the ease in which the refinement and coarsening is performed by the tree data structure. The use of binary tree, i.e. allows refined cells to be added to the domain by simply creating new sub-branches logically below the refined cell. Cell coarsening will become more important especially when a transient solution is needed because the position of various phenomena may change during transient solution. In the present work, the solution-based grid adaptation is carried out using the following algorithm: (1) Start the adaptation process when the selected equation residual is lower than a specified limit. (2) Calculate the indicator parameter for all Cartesian cells. (3) Mark cells for refining and coarsening. (4) Perform coarsening processes on the Cartesian grid. (5) Perform refining processes on the Cartesian grid. (6) Check the balancing constraint.
Fig. 2. Close view of cove region of SKF1.1 wing–flap section.
(7) Interpolate flow properties on the new Cartesian cells. (8) Return to the flow solution algorithm. The error indicator parameter based on variable W, for each cell k is calculated as:
n1 o ek ¼ max l r jDWj
ð1Þ
where DW is the W difference between the left and right cells of edge j. l, is the distance between the centers of the left and right cells of edge j and r determines the weight applies to the W variable difference that was defined in order to prevent over-adaptation near discontinuities like shock waves. For all cases presented in this paper, r is taken to be 2. Finally, the maximum value of the above parameter for all edges of cell k is chosen as the error indicator for that cell. Since the refinement and coarsening of the grid is performed simultaneously in each adaptation cycle, defining the refinement and coarsening threshold has a great impact on the efficiency of the method. In the present work, a statistical approach is used and the standard deviation of error about zero is computed as:
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn 2 k¼1 ek r¼ n
ð2Þ
where n is the total number of Cartesian cells. A cell is then marked for refinement if ek > (pr) and is marked for coarsening if ek < (r/q). In the present work the values of p and q are chosen as 1.5 and 1.75, respectively. Usually 3–6 adaptation cycles are used depending upon the flow features that are present in the domain. In the present work, pressure and Mach number variables are used for shock wave and shear layer detection.
4. Numerical method Two-dimensional compressible viscous flow equations consisting of the mass, momentum, and energy conservation laws can be written in the differential form as:
@W @FI @GI M1 þ þ ¼ @t @x @y Re1 where
"
@FV @GV þ @x @y
# ð3Þ
12
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
3
2
3
2
3
2
q qu qv 7 6 qu 7 6 qu2 þ p 7 6 qv u 7 7 7 6 6 6 I I W¼6 7; F ¼ 6 7; G ¼ 6 2 7; 5 4 qv 5 4 quv 4 qv þ p 5 qE quH qv H 2
3
0
7 7 7; 5
6s 6 xx FV ¼ 6 4 sxy usxx þ v sxy qx
2
0
6 sxy 6 GV ¼ 6 4 syy
usxy þ v syy qy
3
ð4Þ
7 7 7 5
Here P, q, u, v, E and H denote the pressure, density, Cartesian velocity components in x and y directions, total energy and total enthalpy, respectively. The superscripts I and V correspond to the inviscid and viscous terms of the equations. For a perfect gas E ¼ P 1 2 2 P qðc1Þ þ 2 ðu þ v Þ and H ¼ E þ q, where c is the ratio of specific heats. sxx, syy and sxy are the viscous stress tensor components and qx and qy are the heat flux components. 4.1. Finite volume spatial discretization for Cartesian grid The flow equations are solved using a cell-centered finite volume scheme in Cartesian grid zone [9]. The finite volume approximation of the cell k becomes:
Xk
@Wk þ Rk ðWÞ ¼ 0 @t
ð5Þ
where Xk is the area of the cell k and Rk(W) is obtained by evaluating the convective and diffusive flux integrals with a central difference scheme in Eq. (3). In order to prevent oscillations in the neighborhood of shock waves and to provide background dissipation to suppress odd–even modes that are present due to central differencing a blend of first and third order dissipative terms (Dk(W)) is added [9]. The set of ordinary differential equations to be solved can then be written as:
@Wk Xk þ Rk ðWÞ Dk ðWÞ ¼ 0 @t m X 2 2 Dk ðWÞ ¼ aki ekið2Þ ðWp Wk Þ aki eð4Þ ki ðD Wp D Wk Þ
ð6Þ
i¼1
where index i denotes the edge delimiting cells k and p, and m is the number of edges for each cell that in the present work can be between 4 to 8. D2W is obtained summing the differences of conserved variables across the edges forming each cell.
D2 Wk ¼
m X i¼1
ðWp Wk Þ; aki ¼ jui Dyi v i Dxi j þ ci
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dx2i þ Dy2i
Pk P p Pk þ Pp
ð7Þ
e ¼ k2 ðmki Þ; ekið4Þ ¼ maxð0; k4 eð2Þ ki Þ; mki ¼ ð2Þ ki
where ui and vi are the velocity components in the x and y directions and ci is the speed of sound on edge i. k2 and k4 are two empirically chosen constants, which typically have values in the range 1 1 1 < k2 < 1 and 256 < k4 < 32 . 2 4.2. Spatial discretization for mesh-less zone Like finite difference methods, the mesh-less algorithm is applied directly to the differential form of the governing equations. Similar to all mesh-less methods it makes use of a least-square formulation but differs in the way it introduces artificial dissipation, which is essential and necessary for the governing hyperbolic equations. The set of cloud points for a given point i is shown in Fig. 3. Using Taylor’s formula about i to any of its cloud points lead to none square coefficient matrix. Therefore, the least-square method is used to solve this set of equations. The spatial deriva-
Fig. 3. A typical edge (ij) in a least squares cloud.
tives of function u at point i in x and y directions can then be obtained using the least-square method [10] as:
m m X X @u @u ¼ aij ðuj ui Þ; ¼ bij ðuj ui Þ @x i @y i j¼1 j¼1
ð8Þ
where aij and bij are least squares coefficients vectors and m is the number of cloud points for point i. Applying the least-square approximations given by Eq. (8) to each component of flux functions in Eq. (3), a semi-discrete form of flow equations is obtained at point i:
" # m m m m X X @Wi X M1 X I V I V aij DFij þ bij DGij ¼ aij DFij þ bij DGij þ @t Re1 j¼1 j¼1 j¼1 j¼1
ð9Þ
where DFij and DGij are defined as:
DFij ¼ Fj Fi ; DGij ¼ Gj Gi
ð10Þ
A flux H = aF + bG can then be defined in the direction of the leastsquare coefficient vector for an edge ij similar to a directional flux through a face area on an unstructured mesh.
" # m m @Wi X M1 X DHIij ¼ DHVij ; DHij ¼ ðHj Hi Þ þ @t Re1 j¼1 j¼1
ð11Þ
It is well known that such discretization lead to instabilities and must be augmented by stabilizing terms. In this work similar dissipation method as for the Cartesian grid zone is used for mesh-less zone. The dissipation terms consisting a blend of the second and fourth differences of conserved variables (W) that is added to Eq. (11) in order to prevent the oscillations particularly in the neighborhood of shock waves and sharp gradients.
" #! m m X @Wi M1 X DHIij DHVij Di ðWÞ ¼ 0 þ @t Re1 j¼1 j¼1
ð12Þ
Following the above approach the equivalent dissipation terms in mesh-less method are defined as:
Di ðWÞ ¼ ½rðeð2Þ kÞrW r2 ðeð4Þ kÞD2 Wi m X rðeð2Þ kÞrW i ¼ 2 fij ðeð2Þ kÞjþ1 ðWj Wi Þ 2
j¼1
h
i
r2 ðeð4Þ kÞD2 W ¼ 2 i
m X fij ðeð4Þ kÞjþ1 ðD2 Wj D2 Wi Þ 2
j¼1
m m X X D2 Wk ¼ ðWj Wk Þ ¼ ðWj Þ mWk ; fij ¼ j¼1 !
ð13Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a2ij þ bij !
j Lij j
j¼1 !
where Lij is the line vector i to j and j Lij j is the distance between points i and j. Here (2) and (4) are local adaptive coefficients that
13
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
rUjþ12
Uj Ui ~ ¼ rUjþ1 rUjþ1 ~ sij sij 2 2 jlij j
ð17Þ
where rUjþ1 is the average of gradient at mid-point j þ 12 or 2 ! rUi þrUj L and ~ sij is the unit vector between i and j or !ij . 2 j Lij j
4.3. Data exchange between zones For coupling the calculations between mesh-less zone and Cartesian grid zone, the flow variables (W) are computed in the joint area as following:
Fig. 4. Schematic of defined influence domain for zone boundary interpolation.
use the pressure as a sensor to detect sharp gradients. They are formulated as:
ð2Þ ¼ k2 rij ð4Þ ¼ maxð0; k4 ð2Þ Þ
ð14Þ
where rij is a pressure sensor for shocks at any edge ij that is evaluated as:
rij
p p j i ¼ pj þ pi
ð15Þ
and kij is:
! kij ¼ u jþ1 Dxij þ v jþ1 Dyij þ c jþ1 j Lij j 2
2
2
ð16Þ
ijþ1 and v ijþ1 are the velocity components and cijþ1 is the where u
2 2 2 speed of sound at the mid-point j þ 12 of edge ij. Viscous terms appeared in the right-hand side of Eq. (9) are cal
culated at the mid-point station j þ 12 . These terms contain gradients of velocity components and temperature. If the average of these gradients at midpoint are used to calculate the viscous terms, the viscous flux discretization lead to even/odd decoupling. Therefore, it is necessary to modify the gradients at mid-points to remove solution instability. This can be carried out for any variable U as following [8].
(a) In Cartesian grid zone, the variables along boundary faces are interpolated using neighboring mesh-less nodes and halo points that place inside a circle with center k (Fig. 4). The interpolation is carried out using a weighting average on edge as follows:
Pm j¼1 ðxkj Wj Þ Wk ¼ Pm j¼1 ðxkj Þ
ð18Þ
(b) In the mesh-less zone, the halo points (Cartesian cell centers of halo layers) together with mesh-less nodes that situated in the cloud of point L are considered for interpolation of flow variables (Fig. 4). 4.4. Time integration and boundary conditions The integration in time of Eqs. (5) and (12) are obtained using an explicit 4-stage Runge–Kutta scheme where, for economy, the dissipation function (D) is only calculated at the first and third stages. To accelerate convergence, local time stepping and residual smoothing are employed in the present work [9]. For inviscid flow at a solid boundary, no mass or other convective fluxes can penetrate the solid body and no-slip boundary condition is imposed in viscous flow case. Ghost points are used to balance the local clouds of the surface boundary points and improve the accuracy. Solution values at ghost points may be set consistent with physical boundary conditions to obtain a solution. 5. Results To demonstrate the capabilities of the present hybrid method the inviscid and viscous flow computations are carried out around the single and multi-element configurations at transonic flow
Fig. 5. (a) Initial hybrid grid and (b) hybrid adapted grid after four levels of adaptation for NACA 0012, M1 = 0.8, a = 1.25°.
14
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
Fig. 6. Pressure contours over NACA 0012 at M1 = 0.8, a = 1.25° (a) before adaptation and (b) after adaptation.
conditions. A triangular unstructured grid flow solver (Uns2D) has also been used for comparison of results [9]. The computations are carried out on a Pentium PC Dual core with 2.00 GHz speed. The values of k2 and k4 dissipation parameters are set to 0.65 and 0.015 for inviscid test cases.
For the first case a transonic flow with Mach number 0.8 and angle of incidence of 1.25° is considered. Fig. 5 shows the initial and final hybrid grid after four levels of adaptation. The initial grid consisting 1393 points in the mesh-less zone and 2139 cells in the Cartesian grid area. The final grid includes 1393 points in the
Fig. 7. Smooth variation of contours between mesh-less and Cartesian grid zones.
Fig. 8. (a) Surface pressure coefficient distribution and (b) convergence history for NACA 0012 at M1 = 0.8, a = 1.25°.
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17 Table 1 Lift and drag coefficients for NACA 0012 at M1 = 0.8, a = 1.25°.
Present method Uns2D (finite-volume) AGARD [11]
Cl
Cd
0.34904 0.34653 0.359 ± 0.01
0.02303 0.02236 0.0229 ± 0.0008
mesh-less zone and 5802 cells in the Cartesian grid area. The triangular unstructured grid for Uns2D flow solver has 12,678 cells. Fig. 6 shows the pressure contours before and after adaptation. As illustrated the influence of flow adaptation is crucial especially in accurate capturing of shock waves and other flow phenomena that are present in the domain. Smooth variation of contours particularly, between mesh-less and Cartesian grid zones, that are shown in Fig. 7, confirms the validity and accuracy of the proposed method in this area. Fig. 8a shows the surface pressure coefficient distributions obtained from hybrid method computations compared with the results of Uns2D. Fig. 8b illustrates the convergence history for hybrid method in comparison with the Uns2D flow solver. Considerable reduction in computational time (around 60%) has been achieved. The calculated lift and drag coefficients are also shown in Table 1 demonstrating good agreements with the other reliable results. For the second case an inviscid transonic flow with Mach number 0.65 and angle of incidence 2.06° is considered over SKF1.1 multi-element airfoil. Fig. 9a shows the final grid after four levels of adaptation that has 2405 points in the mesh-less zone and 10,784 cells in the Cartesian grid area. Fig. 9b shows the corresponding pressure contours after adaptation. Once again accuracy and smoothness of the results can be observed for a more complex case. Fig. 10 shows the triangular unstructured grid over SKF1.1 wing section that includes 16,361 triangular cells. Fig. 11 illustrates the surface pressure coefficient distributions obtained from present method compared with the results of Uns2D. Desirable agreements have been achieved for surface pressure distributions and shock capturing. The last test case is a transonic laminar flow at Mach number 0.8, angle of attack 10° and Reynolds number 500 over NACA 0012 airfoil. This is the case number A2 of GAMM workshop on numerical simulation of compressible Navier–Stokes flow [12]. In this case the values of 0.25 and 0.015 are used for k2 and k4 dissipation parameters in both zones. Fig. 12a shows the final grid after four levels of adaptation. The hybrid adapted grid contains 3841 nodes in the mesh-less zone and 8529 Cartesian cells. The triangular unstructured grid that is used for Uns2D flow solver has 26,030 cells. Fig. 12b shows the streamlines over the airfoil including a separation zone in the upper surface of the airfoil. The pressure and skin friction coefficients
Fig. 10. Triangular grid over SKF1.1 for Uns2D flow solver.
Fig. 11. Surface pressure coefficient distribution.
Fig. 9. (a) Hybrid adapted grid and (b) pressure contours over SKF1.1 at M1 = 0.65, a = 2.06°.
15
16
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
Fig. 12. (a) Hybrid adapted grid and (b) streamlines of separated flow over NACA 0012 at M1 = 0.8, a = 10° and Re1 = 500.
Fig. 13. (a) Surface pressure coefficient distribution and (b) skin friction coefficient distribution for M1 = 0.8, a = 10° and Re1 = 500.
Fig. 14. (a) Separation point at 37% of chord on the surface and (b) convergence history for NACA 0012 at M1 = 0.8, a = 10° and Re1 = 500.
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
distributions over the airfoil are compared with numerical data [13] in Fig. 13a and b, respectively. The calculated lift and drag coefficients are 0.475 and 0.261 respectively that are again within the range of results reported in Ref. [12]. In particular, the separation point for this flow is about 37% of the chord with maximum 1.0% error compared with the numerical results reported in Ref. [13]. A close view of boundary layer velocity profiles close to the separation point are shown in Fig. 14a. Convergence histories for the new hybrid and fully unstructured control volume (Uns2D) methods are presented in Fig. 14b. As illustrated the computational time is reduced by more than 70% when using the new method compared with the fully unstructured grid finite volume flow solver. 6. Conclusion A hybrid adaptive Cartesian grid with mesh-less zones was developed for calculation of compressible flows. The developed procedure was used for computation of inviscid and viscous flows to validate the accuracy and efficiency of the method. Results indicated good agreements with reliable numerical data. The new method was shown to reduce the computational cost up to 70% especially in viscous flow calculations compared with the fully unstructured flow solver. The future work is focused on extension of the method to the unsteady moving boundary applications. References [1] Thompson JF, Thomas FC, Mastin CW. Automated numerical generation of body fitted curvilinear co-ordinate system for field containing any number of arbitrary 2D bodies. J Comput Phys 1974;15:299–319.
17
[2] Hassan O, Probert EJ, Morgan K. Unstructured mesh procedures for the simulation of three-dimensional transient compressible inviscid flows with moving boundary components. Int J Numer Methods Fluids 1998;27:41–55. [3] Aftosmis MJ. Solution adaptive cartesian grid methods for aerodynamic flows with complex geometries. von Karman Institute for Fluid Dynamics, March 1997, Lecture Series 1997–02. [4] Jahangirian A, Shoraka Y. Adaptive unstructured grid generation for engineering computation of aerodynamic flows. Math Comput Simul 2008;78:627–44. [5] Coirier WJ, Powell KG. Solution-adaptive Cartesian cell approach for viscous and inviscid flows. AIAA J 1996;34(5):938–45. [6] Ding H, Shu C, Yeo KS, Xu D. Development of least-square-based twodimensional finite-difference schemes and their application to simulate natural convection in a cavity. Comput Fluids 2004;33:137–54. [7] Luo H, Baum JD, Lohner R. A hybrid building-block and gridless method for compressible flows. Int J Numer Methods Fluids 2009;59:459–74. [8] Hashemi Y, Jahangirian A. Implicit fully mesh-less method for compressible viscous flow calculations. J Comput Appl Math 2011;235:4687–700. [9] Jahangirian A, Hadidoolabi M. Unstructured moving grids for implicit calculation of unsteady compressible viscous flows. Int J Numer Methods Fluids 2005;47:1107–13. [10] Nguyen VP, Rabczuk T, Bordas S, Duflot MM. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79:763–813. [11] AGARD Fluid Dynamics Panel. Test cases for inviscid flow field methods. AGARD Advisory Report AR-211, AGARD, May 1985. [12] Bristeau M, Glowinski R, Periaux J, Viviand H. Numerical simulation of compressible Navier–Stokes flow. Friedr. Viwey Sohn, Braunschweig, Wiesbaden; 1978. [13] Grasso F, Jameson A, Martinelli L. A multistage multigrid method for the compressible Navier–Stokes equations. In: Bristeau M, Glowinski R, Periaux J, Viviand H, editors. A GAAM-workshop on numerical simulation of compressible Navier–Stokes flow, Nice, France; 1985.
Contents lists available at SciVerse ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Adaptive Cartesian grid with mesh-less zones for compressible flow calculations q A. Jahangirian ⇑, M.Y. Hashemi Department of Aerospace Engineering, Amirkabir University of Technology, 424 Hafez Avenue, Tehran, Iran
a r t i c l e
i n f o
Article history: Received 1 May 2010 Received in revised form 13 August 2011 Accepted 16 August 2011 Available online 6 September 2011 Keywords: Adaptive Cartesian grid Mesh-less method Compressible flow
a b s t r a c t A hybrid grid method that combines the computational efficiency of adaptive Cartesian grids with the flexibility of mesh-less scheme is developed for calculation of compressible flows. The mesh-less zone is created around the geometry by producing layers of nodes along normal direction vectors. Last layer is used as virtual geometry for automatic generation of unstructured Cartesian grid around mesh-less zone. An efficient explicit central difference scheme with artificial dissipation terms and convergence acceleration techniques is developed for both Cartesian grid and mesh-less zones. The method is used for computation of inviscid and viscous flows around single and multi-element airfoils at transonic flow conditions. Results indicate good agreements with those of available unstructured finite-volume flow solver (Uns2D) and other reference numerical data. The new method is shown to reduce the computational cost up to 70% compared with the fully unstructured flow solver. Ó 2011 Elsevier Ltd. All rights reserved.
1. Introduction The main problem of computational fluid dynamics (CFD) for configurations with complex geometry is mesh generation. The structured grid methods are not efficient for such cases [1]. The main advantage of unstructured grid methods is the facility of grid generation for complex configurations [2]. However, the computational costs and memory requirements are generally higher than their structured grid counterparts. Alternatively, the Cartesian grid methods can be used that offer ease of grid generation, lower computational storage requirements, and significantly less operational count per cell compared to body fitted curvilinear grids [3]. The convergence properties of the solver are also potentially better because there are no problems related to skewness or distortion of cells and implementation of high order discretization schemes is simple. Embedded and adaptive grids, can also be used to provide better resolution of geometry and flow features in Cartesian grid [4]. However, the main challenge in using this method deals with arbitrary boundaries. As the grids are not body aligned, Cartesian cells near the body can extend through surfaces of solid components. Hence, accurate means of representations for surface boundary conditions are essential for the success of Cartesian schemes. Different methods have been proposed in the literature to resolve the boundary conditions problem in Cartesian grids. Broadly the methods proposed involve either cut cells for a finite volume treatment [5] or ghost points for a finite difference construction
[6]. A cut-cell method ensures conservation and uses the same flux computation algorithm as for the interior cells. However, the task of computing the volume and fluxes for all of the irregularly shaped cut cells and hanging nodes entail a considerable increase in complexity. More crucially, the method may lead to the creation of very small (irregular) cells at the boundary. Additionally, in the viscous flows the favorable normal lines to the surface and regular cells are not possible to generate which pose problems of numerical stability and computational cost. Recently, mesh-less methods have been developed to conveniently and accurately treats boundary conditions on curved boundaries based on a least-square approach [7]. However, the mesh-less method needs further development in automating the procedure for any geometry particularly in viscous applications. The main objective of the present work is to combine the advantages of recent works regarding the Cartesian grid method [4] and the mesh-less strategy [8] for computation of compressible flows. The method includes generating a layer of nodes around the body while Cartesian grid method is used elsewhere. Since the majority of computational domain is solved using the Cartesian grid method, the resulting procedure is very efficient in terms of both computational cost and storage requirements. In addition, since the mesh-less method is used to obtain the solution for the solid boundary points the developed method can easily be used for complex geometries without need to generate the body-fitted mesh. 2. Mesh-less zone point creation
q
This paper was originally submitted to the ICFD 2010 special issue which was published in Volume 46, Issue 1 of Computers and Fluids. ⇑ Corresponding author. E-mail address: [email protected] (A. Jahangirian). 0045-7930/$ - see front matter Ó 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2011.08.010
In the first step, the surface must be broken up into elements (edges in 2D and faces in 3D) for a full description of a complex geometry. To generate the efficient distribution of points in the
11
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
Fig. 1. Mesh-less zone with its boundary for SKF1.1 two element airfoil.
mesh-less zone around the geometry the unit normal vector for all vertices on the geometry are calculated. At the next step, the layers of nodes are produced along the normal lines until a user specified numbers of layers (Fig. 1). For efficient layer generation the normal vectors can be smoothed where needed and/or truncated where overlapping the layers occurs. The neighbors of any node in the mesh-less zone need to be identified for the discretization of equations. Thus, a node-based data structure is formed which defines the neighbors of each node. In the present work, since the nodes inside the mesh-less zones are generated in a regular manner, the neighbors of each node can be easily defined. In the area that mesh-less zones overlap with each other, neighbors can be obtained by a simple local search. 3. Adaptive Cartesian grid generation Generation of the unstructured Cartesian grid is carried out following the work of Jahangirian and Shoraka [4] which is explained here briefly as follows: (1) Define an outer boundary box for the flow domain. (2) Define surface geometries (here the outer boundary of the mesh-less zone). (3) Subdivide the box and subsequent boxes if they contain or intersect the geometry until the required level of embedded grids is reached. (4) Detect and remove cells inside the mesh-less zone boundary. To complement the mesh-less zone for solution of the governing flow equations a series of halo points are added outside the mesh-less zone. After numerical experiments the first two layers of Cartesian cells are selected and their centers are considered for mesh-less zone computations. The close view of the generated final grid (including halo points) around a wing-flap configuration is shown in Fig. 2. One of the advantages of the Cartesian-based grid generation methods is the ease in which the refinement and coarsening is performed by the tree data structure. The use of binary tree, i.e. allows refined cells to be added to the domain by simply creating new sub-branches logically below the refined cell. Cell coarsening will become more important especially when a transient solution is needed because the position of various phenomena may change during transient solution. In the present work, the solution-based grid adaptation is carried out using the following algorithm: (1) Start the adaptation process when the selected equation residual is lower than a specified limit. (2) Calculate the indicator parameter for all Cartesian cells. (3) Mark cells for refining and coarsening. (4) Perform coarsening processes on the Cartesian grid. (5) Perform refining processes on the Cartesian grid. (6) Check the balancing constraint.
Fig. 2. Close view of cove region of SKF1.1 wing–flap section.
(7) Interpolate flow properties on the new Cartesian cells. (8) Return to the flow solution algorithm. The error indicator parameter based on variable W, for each cell k is calculated as:
n1 o ek ¼ max l r jDWj
ð1Þ
where DW is the W difference between the left and right cells of edge j. l, is the distance between the centers of the left and right cells of edge j and r determines the weight applies to the W variable difference that was defined in order to prevent over-adaptation near discontinuities like shock waves. For all cases presented in this paper, r is taken to be 2. Finally, the maximum value of the above parameter for all edges of cell k is chosen as the error indicator for that cell. Since the refinement and coarsening of the grid is performed simultaneously in each adaptation cycle, defining the refinement and coarsening threshold has a great impact on the efficiency of the method. In the present work, a statistical approach is used and the standard deviation of error about zero is computed as:
sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Pn 2 k¼1 ek r¼ n
ð2Þ
where n is the total number of Cartesian cells. A cell is then marked for refinement if ek > (pr) and is marked for coarsening if ek < (r/q). In the present work the values of p and q are chosen as 1.5 and 1.75, respectively. Usually 3–6 adaptation cycles are used depending upon the flow features that are present in the domain. In the present work, pressure and Mach number variables are used for shock wave and shear layer detection.
4. Numerical method Two-dimensional compressible viscous flow equations consisting of the mass, momentum, and energy conservation laws can be written in the differential form as:
@W @FI @GI M1 þ þ ¼ @t @x @y Re1 where
"
@FV @GV þ @x @y
# ð3Þ
12
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
3
2
3
2
3
2
q qu qv 7 6 qu 7 6 qu2 þ p 7 6 qv u 7 7 7 6 6 6 I I W¼6 7; F ¼ 6 7; G ¼ 6 2 7; 5 4 qv 5 4 quv 4 qv þ p 5 qE quH qv H 2
3
0
7 7 7; 5
6s 6 xx FV ¼ 6 4 sxy usxx þ v sxy qx
2
0
6 sxy 6 GV ¼ 6 4 syy
usxy þ v syy qy
3
ð4Þ
7 7 7 5
Here P, q, u, v, E and H denote the pressure, density, Cartesian velocity components in x and y directions, total energy and total enthalpy, respectively. The superscripts I and V correspond to the inviscid and viscous terms of the equations. For a perfect gas E ¼ P 1 2 2 P qðc1Þ þ 2 ðu þ v Þ and H ¼ E þ q, where c is the ratio of specific heats. sxx, syy and sxy are the viscous stress tensor components and qx and qy are the heat flux components. 4.1. Finite volume spatial discretization for Cartesian grid The flow equations are solved using a cell-centered finite volume scheme in Cartesian grid zone [9]. The finite volume approximation of the cell k becomes:
Xk
@Wk þ Rk ðWÞ ¼ 0 @t
ð5Þ
where Xk is the area of the cell k and Rk(W) is obtained by evaluating the convective and diffusive flux integrals with a central difference scheme in Eq. (3). In order to prevent oscillations in the neighborhood of shock waves and to provide background dissipation to suppress odd–even modes that are present due to central differencing a blend of first and third order dissipative terms (Dk(W)) is added [9]. The set of ordinary differential equations to be solved can then be written as:
@Wk Xk þ Rk ðWÞ Dk ðWÞ ¼ 0 @t m X 2 2 Dk ðWÞ ¼ aki ekið2Þ ðWp Wk Þ aki eð4Þ ki ðD Wp D Wk Þ
ð6Þ
i¼1
where index i denotes the edge delimiting cells k and p, and m is the number of edges for each cell that in the present work can be between 4 to 8. D2W is obtained summing the differences of conserved variables across the edges forming each cell.
D2 Wk ¼
m X i¼1
ðWp Wk Þ; aki ¼ jui Dyi v i Dxi j þ ci
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Dx2i þ Dy2i
Pk P p Pk þ Pp
ð7Þ
e ¼ k2 ðmki Þ; ekið4Þ ¼ maxð0; k4 eð2Þ ki Þ; mki ¼ ð2Þ ki
where ui and vi are the velocity components in the x and y directions and ci is the speed of sound on edge i. k2 and k4 are two empirically chosen constants, which typically have values in the range 1 1 1 < k2 < 1 and 256 < k4 < 32 . 2 4.2. Spatial discretization for mesh-less zone Like finite difference methods, the mesh-less algorithm is applied directly to the differential form of the governing equations. Similar to all mesh-less methods it makes use of a least-square formulation but differs in the way it introduces artificial dissipation, which is essential and necessary for the governing hyperbolic equations. The set of cloud points for a given point i is shown in Fig. 3. Using Taylor’s formula about i to any of its cloud points lead to none square coefficient matrix. Therefore, the least-square method is used to solve this set of equations. The spatial deriva-
Fig. 3. A typical edge (ij) in a least squares cloud.
tives of function u at point i in x and y directions can then be obtained using the least-square method [10] as:
m m X X @u @u ¼ aij ðuj ui Þ; ¼ bij ðuj ui Þ @x i @y i j¼1 j¼1
ð8Þ
where aij and bij are least squares coefficients vectors and m is the number of cloud points for point i. Applying the least-square approximations given by Eq. (8) to each component of flux functions in Eq. (3), a semi-discrete form of flow equations is obtained at point i:
" # m m m m X X @Wi X M1 X I V I V aij DFij þ bij DGij ¼ aij DFij þ bij DGij þ @t Re1 j¼1 j¼1 j¼1 j¼1
ð9Þ
where DFij and DGij are defined as:
DFij ¼ Fj Fi ; DGij ¼ Gj Gi
ð10Þ
A flux H = aF + bG can then be defined in the direction of the leastsquare coefficient vector for an edge ij similar to a directional flux through a face area on an unstructured mesh.
" # m m @Wi X M1 X DHIij ¼ DHVij ; DHij ¼ ðHj Hi Þ þ @t Re1 j¼1 j¼1
ð11Þ
It is well known that such discretization lead to instabilities and must be augmented by stabilizing terms. In this work similar dissipation method as for the Cartesian grid zone is used for mesh-less zone. The dissipation terms consisting a blend of the second and fourth differences of conserved variables (W) that is added to Eq. (11) in order to prevent the oscillations particularly in the neighborhood of shock waves and sharp gradients.
" #! m m X @Wi M1 X DHIij DHVij Di ðWÞ ¼ 0 þ @t Re1 j¼1 j¼1
ð12Þ
Following the above approach the equivalent dissipation terms in mesh-less method are defined as:
Di ðWÞ ¼ ½rðeð2Þ kÞrW r2 ðeð4Þ kÞD2 Wi m X rðeð2Þ kÞrW i ¼ 2 fij ðeð2Þ kÞjþ1 ðWj Wi Þ 2
j¼1
h
i
r2 ðeð4Þ kÞD2 W ¼ 2 i
m X fij ðeð4Þ kÞjþ1 ðD2 Wj D2 Wi Þ 2
j¼1
m m X X D2 Wk ¼ ðWj Wk Þ ¼ ðWj Þ mWk ; fij ¼ j¼1 !
ð13Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 a2ij þ bij !
j Lij j
j¼1 !
where Lij is the line vector i to j and j Lij j is the distance between points i and j. Here (2) and (4) are local adaptive coefficients that
13
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
rUjþ12
Uj Ui ~ ¼ rUjþ1 rUjþ1 ~ sij sij 2 2 jlij j
ð17Þ
where rUjþ1 is the average of gradient at mid-point j þ 12 or 2 ! rUi þrUj L and ~ sij is the unit vector between i and j or !ij . 2 j Lij j
4.3. Data exchange between zones For coupling the calculations between mesh-less zone and Cartesian grid zone, the flow variables (W) are computed in the joint area as following:
Fig. 4. Schematic of defined influence domain for zone boundary interpolation.
use the pressure as a sensor to detect sharp gradients. They are formulated as:
ð2Þ ¼ k2 rij ð4Þ ¼ maxð0; k4 ð2Þ Þ
ð14Þ
where rij is a pressure sensor for shocks at any edge ij that is evaluated as:
rij
p p j i ¼ pj þ pi
ð15Þ
and kij is:
! kij ¼ u jþ1 Dxij þ v jþ1 Dyij þ c jþ1 j Lij j 2
2
2
ð16Þ
ijþ1 and v ijþ1 are the velocity components and cijþ1 is the where u
2 2 2 speed of sound at the mid-point j þ 12 of edge ij. Viscous terms appeared in the right-hand side of Eq. (9) are cal
culated at the mid-point station j þ 12 . These terms contain gradients of velocity components and temperature. If the average of these gradients at midpoint are used to calculate the viscous terms, the viscous flux discretization lead to even/odd decoupling. Therefore, it is necessary to modify the gradients at mid-points to remove solution instability. This can be carried out for any variable U as following [8].
(a) In Cartesian grid zone, the variables along boundary faces are interpolated using neighboring mesh-less nodes and halo points that place inside a circle with center k (Fig. 4). The interpolation is carried out using a weighting average on edge as follows:
Pm j¼1 ðxkj Wj Þ Wk ¼ Pm j¼1 ðxkj Þ
ð18Þ
(b) In the mesh-less zone, the halo points (Cartesian cell centers of halo layers) together with mesh-less nodes that situated in the cloud of point L are considered for interpolation of flow variables (Fig. 4). 4.4. Time integration and boundary conditions The integration in time of Eqs. (5) and (12) are obtained using an explicit 4-stage Runge–Kutta scheme where, for economy, the dissipation function (D) is only calculated at the first and third stages. To accelerate convergence, local time stepping and residual smoothing are employed in the present work [9]. For inviscid flow at a solid boundary, no mass or other convective fluxes can penetrate the solid body and no-slip boundary condition is imposed in viscous flow case. Ghost points are used to balance the local clouds of the surface boundary points and improve the accuracy. Solution values at ghost points may be set consistent with physical boundary conditions to obtain a solution. 5. Results To demonstrate the capabilities of the present hybrid method the inviscid and viscous flow computations are carried out around the single and multi-element configurations at transonic flow
Fig. 5. (a) Initial hybrid grid and (b) hybrid adapted grid after four levels of adaptation for NACA 0012, M1 = 0.8, a = 1.25°.
14
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
Fig. 6. Pressure contours over NACA 0012 at M1 = 0.8, a = 1.25° (a) before adaptation and (b) after adaptation.
conditions. A triangular unstructured grid flow solver (Uns2D) has also been used for comparison of results [9]. The computations are carried out on a Pentium PC Dual core with 2.00 GHz speed. The values of k2 and k4 dissipation parameters are set to 0.65 and 0.015 for inviscid test cases.
For the first case a transonic flow with Mach number 0.8 and angle of incidence of 1.25° is considered. Fig. 5 shows the initial and final hybrid grid after four levels of adaptation. The initial grid consisting 1393 points in the mesh-less zone and 2139 cells in the Cartesian grid area. The final grid includes 1393 points in the
Fig. 7. Smooth variation of contours between mesh-less and Cartesian grid zones.
Fig. 8. (a) Surface pressure coefficient distribution and (b) convergence history for NACA 0012 at M1 = 0.8, a = 1.25°.
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17 Table 1 Lift and drag coefficients for NACA 0012 at M1 = 0.8, a = 1.25°.
Present method Uns2D (finite-volume) AGARD [11]
Cl
Cd
0.34904 0.34653 0.359 ± 0.01
0.02303 0.02236 0.0229 ± 0.0008
mesh-less zone and 5802 cells in the Cartesian grid area. The triangular unstructured grid for Uns2D flow solver has 12,678 cells. Fig. 6 shows the pressure contours before and after adaptation. As illustrated the influence of flow adaptation is crucial especially in accurate capturing of shock waves and other flow phenomena that are present in the domain. Smooth variation of contours particularly, between mesh-less and Cartesian grid zones, that are shown in Fig. 7, confirms the validity and accuracy of the proposed method in this area. Fig. 8a shows the surface pressure coefficient distributions obtained from hybrid method computations compared with the results of Uns2D. Fig. 8b illustrates the convergence history for hybrid method in comparison with the Uns2D flow solver. Considerable reduction in computational time (around 60%) has been achieved. The calculated lift and drag coefficients are also shown in Table 1 demonstrating good agreements with the other reliable results. For the second case an inviscid transonic flow with Mach number 0.65 and angle of incidence 2.06° is considered over SKF1.1 multi-element airfoil. Fig. 9a shows the final grid after four levels of adaptation that has 2405 points in the mesh-less zone and 10,784 cells in the Cartesian grid area. Fig. 9b shows the corresponding pressure contours after adaptation. Once again accuracy and smoothness of the results can be observed for a more complex case. Fig. 10 shows the triangular unstructured grid over SKF1.1 wing section that includes 16,361 triangular cells. Fig. 11 illustrates the surface pressure coefficient distributions obtained from present method compared with the results of Uns2D. Desirable agreements have been achieved for surface pressure distributions and shock capturing. The last test case is a transonic laminar flow at Mach number 0.8, angle of attack 10° and Reynolds number 500 over NACA 0012 airfoil. This is the case number A2 of GAMM workshop on numerical simulation of compressible Navier–Stokes flow [12]. In this case the values of 0.25 and 0.015 are used for k2 and k4 dissipation parameters in both zones. Fig. 12a shows the final grid after four levels of adaptation. The hybrid adapted grid contains 3841 nodes in the mesh-less zone and 8529 Cartesian cells. The triangular unstructured grid that is used for Uns2D flow solver has 26,030 cells. Fig. 12b shows the streamlines over the airfoil including a separation zone in the upper surface of the airfoil. The pressure and skin friction coefficients
Fig. 10. Triangular grid over SKF1.1 for Uns2D flow solver.
Fig. 11. Surface pressure coefficient distribution.
Fig. 9. (a) Hybrid adapted grid and (b) pressure contours over SKF1.1 at M1 = 0.65, a = 2.06°.
15
16
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
Fig. 12. (a) Hybrid adapted grid and (b) streamlines of separated flow over NACA 0012 at M1 = 0.8, a = 10° and Re1 = 500.
Fig. 13. (a) Surface pressure coefficient distribution and (b) skin friction coefficient distribution for M1 = 0.8, a = 10° and Re1 = 500.
Fig. 14. (a) Separation point at 37% of chord on the surface and (b) convergence history for NACA 0012 at M1 = 0.8, a = 10° and Re1 = 500.
A. Jahangirian, M.Y. Hashemi / Computers & Fluids 54 (2012) 10–17
distributions over the airfoil are compared with numerical data [13] in Fig. 13a and b, respectively. The calculated lift and drag coefficients are 0.475 and 0.261 respectively that are again within the range of results reported in Ref. [12]. In particular, the separation point for this flow is about 37% of the chord with maximum 1.0% error compared with the numerical results reported in Ref. [13]. A close view of boundary layer velocity profiles close to the separation point are shown in Fig. 14a. Convergence histories for the new hybrid and fully unstructured control volume (Uns2D) methods are presented in Fig. 14b. As illustrated the computational time is reduced by more than 70% when using the new method compared with the fully unstructured grid finite volume flow solver. 6. Conclusion A hybrid adaptive Cartesian grid with mesh-less zones was developed for calculation of compressible flows. The developed procedure was used for computation of inviscid and viscous flows to validate the accuracy and efficiency of the method. Results indicated good agreements with reliable numerical data. The new method was shown to reduce the computational cost up to 70% especially in viscous flow calculations compared with the fully unstructured flow solver. The future work is focused on extension of the method to the unsteady moving boundary applications. References [1] Thompson JF, Thomas FC, Mastin CW. Automated numerical generation of body fitted curvilinear co-ordinate system for field containing any number of arbitrary 2D bodies. J Comput Phys 1974;15:299–319.
17
[2] Hassan O, Probert EJ, Morgan K. Unstructured mesh procedures for the simulation of three-dimensional transient compressible inviscid flows with moving boundary components. Int J Numer Methods Fluids 1998;27:41–55. [3] Aftosmis MJ. Solution adaptive cartesian grid methods for aerodynamic flows with complex geometries. von Karman Institute for Fluid Dynamics, March 1997, Lecture Series 1997–02. [4] Jahangirian A, Shoraka Y. Adaptive unstructured grid generation for engineering computation of aerodynamic flows. Math Comput Simul 2008;78:627–44. [5] Coirier WJ, Powell KG. Solution-adaptive Cartesian cell approach for viscous and inviscid flows. AIAA J 1996;34(5):938–45. [6] Ding H, Shu C, Yeo KS, Xu D. Development of least-square-based twodimensional finite-difference schemes and their application to simulate natural convection in a cavity. Comput Fluids 2004;33:137–54. [7] Luo H, Baum JD, Lohner R. A hybrid building-block and gridless method for compressible flows. Int J Numer Methods Fluids 2009;59:459–74. [8] Hashemi Y, Jahangirian A. Implicit fully mesh-less method for compressible viscous flow calculations. J Comput Appl Math 2011;235:4687–700. [9] Jahangirian A, Hadidoolabi M. Unstructured moving grids for implicit calculation of unsteady compressible viscous flows. Int J Numer Methods Fluids 2005;47:1107–13. [10] Nguyen VP, Rabczuk T, Bordas S, Duflot MM. Meshless methods: a review and computer implementation aspects. Math Comput Simul 2008;79:763–813. [11] AGARD Fluid Dynamics Panel. Test cases for inviscid flow field methods. AGARD Advisory Report AR-211, AGARD, May 1985. [12] Bristeau M, Glowinski R, Periaux J, Viviand H. Numerical simulation of compressible Navier–Stokes flow. Friedr. Viwey Sohn, Braunschweig, Wiesbaden; 1978. [13] Grasso F, Jameson A, Martinelli L. A multistage multigrid method for the compressible Navier–Stokes equations. In: Bristeau M, Glowinski R, Periaux J, Viviand H, editors. A GAAM-workshop on numerical simulation of compressible Navier–Stokes flow, Nice, France; 1985.