- Email: [email protected]
Overset adaptive-grid method with applications to compressible flows
PII:
Computers & Fluids Vol. 27, Nos 5±6, pp. 599±610, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain S0045-7930(97)00059-5 0045-7930/98 $19.00 + 0.00
OVERSET ADAPTIVE-GRID METHOD WITH APPLICATIONS TO COMPRESSIBLE FLOWS KENICHI MATSUNO{, MASASHI YAMAKAWA and NOBUYUKI SATOFUKA Department of Mechanical and System Engineering, Kyoto Institute of Technology Matsugasaki, Sakyo-ku, Kyoto 606 Japan (Received 6 March 1997) AbstractÐA new overset adaptive-grid method which combines a solution-adaptive grid method with an overset-grid method is proposed and applied to transonic ¯ows in this paper. The position on which the subgrid is located is adjusted to the region where the gradient of ¯ow properties is large. The solution-adaptive grid method using an elliptic equation is applied to the overset subgrid. Thus the present overset adaptive-grid method generates the grid adapted twice to the ¯ow solution. Two new ideas are also proposed for the basic adaptive-grid method of the elliptic type in this paper. One of the new ideas is to use a modi®ed distribution of ¯ow property instead of the actual one in order to estimate weighting functions which control grid spacing. The other is a locally-variable relaxation parameter technique which locally changes the relaxation parameter in the Jacobi iteration process according to the value of the coecient of the ®rst order derivatives in the elliptic equation. The present overset adaptive-grid method is applied to transonic ¯ows, and is shown to be very promising, especially for sharp shock capturing. # 1998 Elsevier Science Ltd. All rights reserved
1. INTRODUCTION
Developments of computers and numerical methods have made it possible to simulate practical ¯ow problems with reasonable accuracy. However, since the computer memory and processor speed are still limited we cannot use enough grid points for complex ¯ow problems. It is also well known that the accuracy and resolution of the numerical solutions are strongly dependent on the grid generated in the ¯ow®eld. Thus it is very important in Computational Fluid Dynamics to generate a grid of high quality. In generating a grid of high quality, there exist two important issues. One is how to eciently ®t the grid to the body boundaries of complex geometry. The overset (Chimera) grid method [1] is the one of the strategies of that problem. The overset-grid method ®rstly generates the main grid which covers the whole of the ¯ow®eld including the bodies, and secondly overlays the body-®tted subgrid on the main grid. This method has been successfully applied to various ¯ow problems in complex geometry [2]. The other important issue is how to concentrate the grid in the high-gradient regions of the ¯ow®eld. A solution-adaptive method is the one which make it possible to concentrate the grid automatically in such regions. There are two approaches of the grid adaption: moving grid and grid re®nement. The moving grid method is simpler in the data structure and easier to implement in the existing ¯ow solvers. There are many solution-adaptive methods of the moving-grid approach [3]. Among them, the adaptive methods using a system of elliptic equations are relatively popular since the elliptic equation generates smooth grid distribution due to the inherent maximum principle. The solution-adaptive grid method of elliptic type uses the source term to control the grid distribution so that the grid density is large at the large-gradient region. One of the drawbacks of the adaptive grid method of elliptic type is the grid skewness. The solution adaptive grid method generates possibly highly skewed grid, for example, in the case when the shock wave stands in the diagonal direction of the mesh lines. This often happens in the single grid structure. The overset-grid method can be used to achieve a ®ne grid locally. That is, the subgrid is overlaid on the region where the gradient of the ¯ow property is large on the main grid. This idea uses the over{Author for correspondence. 599
600
Kenichi Matsuno et al.
set-grid method in a solution-adaptive manner. A similar idea was used by Fujii [4] for tracking blast waves. Kao et al. [5] proposed the idea to use of the overset Chimera grid in large gradient regions to capture the salient features accurately during computations. These authors successfully applied the overset-grid method to ¯uid dynamics problems. However, they merely overlaid the subgrid in large gradient regions. They did not apply any moving-grid techniques to the subgrids and main grid. The purpose of this paper is to propose a new `overset adaptive-grid method'. The basic idea of this method is to apply the solution-adaptive grid method to the subgrid which is already overlaid on the high-gradient region. In the procedure, the elliptic adaptive-grid method proposed by Matsuno and Dwyer [6] is successfully applied to the overset subgrid. Two new ideas are also proposed to develop the present adaptive-grid method. One of the new ideas is to use a modi®ed distribution of ¯ow properties instead of the actual distribution in order to estimate weighting (control) functions which control grid spacing. This approach makes it possible to adapt the grid points in the high-gradient regions more smoothly and eciently. The other is a locally-variable relaxation parameter technique. The adaptive grid generation method using elliptic equations can possibly diverge when the value of control function becomes large in the iteration process. In that case, an under-relaxation is eective to prevent divergence. Thus, the locally-variable relaxation parameter method is proposed in this paper. This method is changing relaxation parameter locally from under-relaxation to over-relaxation due to the value of the control functions. These two new methods are successfully introduced in the adaptive grid generation procedures. Overset adaptive-grid solutions of transonic ¯ows are presented in this paper. Presently, the focus of adaptive-griding is on the regions of the shock waves. However, the present procedure is easily extended to more general cases. 2. ADAPTIVE-GRID METHOD
The present adaptive-grid method is based on the method proposed by Matsuno and Dwyer [6]. Here we describe the method brie¯y using one-dimensional space and after that we extend the method to two dimensions. Now we consider the one-dimensional space, whose coordinate is denoted by s. It is necessary for the adaptive-grid method to concentrate the grid in the region where the derivative of ¯ow physics is large. This property is described by the equation wDsi const:
1
where w is the weight function and Dsi=si + 1ÿsi is the grid spacing at grid index i. In this paper the following weight function, which consists of a ®rst derivative of the ¯ow property f, is adopted q 2 w 1 b fs 2 : Here, f is any ¯ow property such as pressure p, density r, Mach number M, and so on. The coecient b is a user-speci®ed free parameter, which is related to grid spacing due to the eect of ¯ow change. Equation (1) may be considered as a ®nite-dierence equation of the following dierential equation when we regard the term Dsi as Dsi=Dsi/Dx with Dx = 1. wsx const:
3
To avoid the estimation of the constant in the right-hand side of Equation (3), we dierentiate Equation (3) with respect to x. Then we get the following one-dimensional equation: wx 4 sx 0: sxx w Equation (4) is the dierential form for Equation (1). Thompson [7] developed a grid generation method using elliptic equations. This elliptic grid generation method has been applied to various types of geometry. The elliptic equation pro-
Overset adaptive-grid method
601
posed by Thomas and Middleco [8] is written in the one-dimensional form as: xss xs 2 P x:
5
Here P is a so-called `control function', which controls grid point distribution. This equation can be transformed to x coordinate by interchanging the role of dependent and independent variables. This yields the following elliptic equation for one-dimensional grid generation sxx Psx 0:
6
Comparing Equation (4) with Equation (6), the control function P relates with weight function w as follows: P wx =w:
7
Thus we get the elliptic equation for one-dimensional adaptive-grid generation sxx Psx 0, P wx =w
8
For a ¯ow-solution method, it is unfavorable from the viewpoint of solution-accuracy that the grid-spacing ratio of three consecutive points changes suddenly. Thus, it is natural that the user speci®es an allowable limit of grid-spacing ratio of the adaptive grid. Let ri be the grid-spacing ratio at index i ri Dsi =Dsiÿ1 , Dsi si1 ÿ si :
9
Now, let K r1 be the limit of the allowable grid-spacing ratio, then we get the following inequality 1=KRri RK:
10
Once we specify the allowable limit of grid-spacing ratio K, the parameter b in Equation (2) is no more free. The value of b must satisfy the inequality at grid point i 2 3 4 Kÿ1 K1 11 bR4 ÿ 2 ÿ 5: fs ÿ 4 Kÿ1 fs 2 x
K1
i
We have only one value of b. Therefore b is determined by taking the minimum of the largest bi, that is 12 b min bi : i 1, 2, . . . , imax where
2
3 4 Kÿ1 bi 4 ÿ 2 K1 ÿ 5: fs ÿ 4 Kÿ1 fs 2 x
K1
13
i
The extension of the method described above to two-dimensional is done by applying the one-dimensional formulation straightforwardly to the two arclength directions de®ned by the generalized coordinates x and Z. Let s and s' be the arclength along the Z-constant and x-constant grid lines, respectively. The basic criterion for the grid adaption is wsx const:, s2x x 2x y2x w 0 sZ 0 const:, s2Z 0 x 2Z y2Z : The two-dimensional forms of the Poisson equation [8] are ÿ xxx xyy P x, Z x2x x2y
14
602
Kenichi Matsuno et al.
ÿ Zxx Zyy Q x, Z Z2x Z2y
15
where P and Q are control functions de®ned by weight functions w and w' estimated respectively along the arclengths s(x) and s'(Z). P wx =w Q wZ 0 =w 0 :
16
Equation (15) can be inverted to give a form which can be used for numerical solution to obtain the grid ÿ ÿ 17 a rxx Prx ÿ 2brxZ g rZZ QrZ 0 where r = (x, y) and a x 2Z y2Z b x x x Z yx yZ g x 2x y2x :
18
The weight functions w and w' are de®ned the same manner as described in one-dimensional formulation. The ®rst derivative formulations along the two arclength directions are q ÿ 2 w 1 b fs w0
q ÿ 2 1 b 0 fs 0 :
19
The parameters b and b' are determined by specifying the maximum limits K and K' of the gridspacing ratio in two directions, respectively. When the control functions P and Q include the geometrical information beforehand, or if we desire to retain the eect of the solution-independent grid point distribution (usually an initial grid) determined by Pg and so on, then the eect of solution adaption may be included by de®ning control function as combination of the original-grid part Pg and adaptive-grid part Pa. The present control function, for example P, is composed of an original-grid part Pg and an adaptive-grid part Pa. P 1 ÿ yPg yPa , 0RyR1
20
where y is a parameter which controls the weight between the original and adaptive grids property. If the initial (original) grid is generated by another type of grid generation method, then we need to incorporate the information of the original grid distribution in the control functions Pg and Qg. To obtain the information of the initial grid distribution, we backsolve the control function Pi, j and Qi, j in Equation (17) as: Pgi, j
ÿ 1 ÿ yZ ÿ ax xx 2bx xZ ÿ gx ZZ ÿ x Z ÿ ayxx 2byxZ ÿ gyZZ Ja
Qgi, j
ÿ 1 ÿ yx ax xx ÿ 2bx xZ gx ZZ x x ÿ ayxx 2byxZ ÿ gyZZ Jg
21
where J = xxyZÿxZyx and derivatives of the right side are estimated at the grid point (i, j) by ®nite dierence approximations. To develop the above mentioned adaptive-grid method [6], two new techniques are proposed and introduced in this paper:
Overset adaptive-grid method
603
1. modi®ed distribution of the ¯ow properties instead of the actual distribution; 2. locally-variable relaxation parameter techniques for iterative solutions of elliptic equations.
3. MODIFIED DISTRIBUTION OF FLOW-PROPERTY
The use of the modi®ed distribution of the ¯ow properties is motivated by the consideration of making the variation of the control function smooth when the high-resolution Euler solver is applied to get the ¯ow solution. When we apply the high resolution scheme to solve inviscid (Euler) part of ¯uid dynamics equations, the shock wave is captured in two or three grid points. In this case, the control functions have large value only in a narrow region, and zero value in the other region. These control function cannot realize smooth concentration of the grid points. To avoid such situation, the present adaptive method uses a modi®ed distribution of the ¯ow variables. The modi®cation of the distribution of the ¯ow variables means in this paper to smear the obtained ¯ow solution. Suppose fi, j is an actual solution at the grid point (i, j), then the modi®ed solution f~i, j is obtained by averaging N 1 ÿfX fp , 0RfR1 f~i, j ffi, j N p1
22
where fp is a value of f at the pth grid point of N points surrounding (i, j), and f is the weight of the averaging. It is eective in getting the modi®ed smooth solution that the above averaging procedure is repeated several times. The control function Pa is estimated using the distribution of the modi®ed f~i, j .
4. LOCALLY-VARIABLE RELAXATION PARAMETER
The adaptive-grid generation Equation (17) can be solved by a relaxation method, usually the line Jacobi method, and written symbolically as Ldr ÿR r, rm1 rm o dr
23
where m is the index of iteration, and o is relaxation parameter of usually 1 Ro < 2. The under-relaxation parameter is eective for robustness in the present adaptive method when the elliptic Poisson equation is solved iteratively. It decreases the grid movement in each iteration process, especially in the case of extremely high-gradient regions such as shock wave presents. However, the under-relaxation parameter decreases the convergence rate. Hence we use the under-relaxation parameter locally and total convergence rate can be avoided decreasing. This is similar to the use of the local-time step techniques in the ¯ow solution method. In this paper, the relaxation parameter o is determined according to the local values of Pi, j and Qi, j as Ldr ÿR r, rm1 rm o i, j dr, o i, j func Pi, j , Qi, j , 0
24
25
5. OVERSET ADAPTIVE-GRID METHOD
The overset(Chimera) grid method [1] was originally developed to adjust the grid smoothly to the body surface in complex multi-body geometry by oversetting local subgrid around bodies on
604
Kenichi Matsuno et al.
the main grid covering whole of the ¯ow ®eld. The overset method can be used to achieve a ®ne grid locally. The subgrid is overlaid in the high-gradient region. This idea is to use the overset method in a solution-adaptive manner. In this paper, the following overset adaptive-grid method is proposed. First, we overlay the subgrid in the high-gradient regions. Next we apply the solution-adaptive method to the overset subgrid. Thus, the present method adapts the grid twice to the ¯ow solution. In other words, the overset-grid method is combined with the adaptive-grid method in order to get the solution adaptive grid of totally high quality and eciency. The adaptive-grid method on a single grid structure is often bothered with grid skewness. It is hard to generate orthogonal adaptive grid in multi-dimensions. However, the present overset adaptive-grid method can avoid such bad skewness by overlaying subgrid that adjusts to the ¯ow direction or gradient. The present research eort is motivated partly by consideration to avoid serious grid skewness in the adaptive griding. The present overset adaptive-grid generation procedure is as follows. The body-®tted coarse main grid is ®rstly generated by the elliptic grid generation method in the usual manner. The ¯ow is solved on the coarse main grid as a base ¯ow for overset griding and adaptive griding. The regions with large gradients of ¯ow properties, such as shock waves, are detected. The subgrids which cover the large-gradient regions are generated and put on the large-gradient regions in the main grid. In this step, each subgrid is not adapted to the ¯ow solution but merely put on the large-gradient region. Next, the solution-adaptive method using the elliptic equation described above is applied to the subgrids. Finally the ¯ow®eld is solved again on the main grid and subgrids simultaneously. The data transfer between the main grid and subgrid is done by bi-linear interpolation. To get the desirable adaptive properties, the solution-adaptive grid generation procedure is usually repeated or iterated twice or three times.
6. THE EULER EQUATION AND FLOW SOLVER
The two-dimensional compressible Euler equations written in conservation law forms are numerically solved by a cell-centered ®nite volume scheme. The scheme incorporates Roe's approximate Riemann solvers [9] for estimating ¯ux at cell interface. The time integration is performed by a diagonal ADI scheme or explicit forward Euler method. The space accuracy is increased by introducing MUSCL interpolation techniques with a minmod limiter. The present computations assume steady-state solutions. Thus, convergence acceleration techniques, for example, local time stepping, are introduced in the code.
7. APPLICATIONS
First of all, the two newly introduced techniques, the use of modi®ed ¯ow solution for the estimation of the weight function and the locally-variable relaxation parameter in the iteration, were tested using the same problem as of Ref. [6]. In this test, the adaptive solution was obtained without troubles even for the cases that the divergence were experienced in Ref. [6]. Next the test problem is the nonlinear Burgers' equation ÿ 26 ut u2 =2 x 0, 0RxR1 with the boundary condition, x = 0:u = 1 and x = 1:u = ÿ 1. The initial condition and exact steady-state solution with t 4 1 are 8 1 0RxR0:25 < 1 0RxR0:5 0:25RxR0:75 u x, 1 27 u x, 0 1 ÿ 2 xÿ0:25 0:5 ÿ1 0:5RxR1 : ÿ1 0:75RxR1
Overset adaptive-grid method
605
Fig. 1. Actual and modi®ed solution for Burgers'equation.
Note that the Roe's ¯ux-dierence splitting scheme of ®rst order accuracy captures discontinuity solution at steady state in one grid spacing in this case. The adaptive grid is generated for the steady-state solution. The modi®ed solution for estimating weight function is determined by applying the smoothing procedure, Equation (22), ten times with f = 0.5. Figure 1 shows the actual solution and its modi®ed solution. Figure 2(a) shows the adaptive grid distribution, where the weight function is estimated using actual solution. Figure 2(b) shows the adaptive grid distribution, where the modi®ed solution is used to estimate the weight function. The ®rst derivative is used to estimate the weighting function for adaptive griding. From these ®gures, it
Fig. 2. Adaptive grid for the Burgers' equation (21 points): (a) Pa by actual solution; and (b) Pa by modi®ed solution.
606
Kenichi Matsuno et al.
Fig. 3. Transonic ¯ow through channel with 108 thick circlar arc. Cp, Contours, Min¯ow=0.675. (a) Single grid solution (70 24), (b) overset grid solution (49 17 for the main grid, 20 17, 16 16, 16 16 for the subgrids), (c) overset adaptive-grid solution [the numbers of grid points are the same as (b)].
is shown that the use of modi®ed solution for estimation of weight function is very promising for adaptive griding. Next we show the results of the overset adaptive-grid method applied to transonic ¯ows. The present Euler solver is the cell-centered ®nite volume scheme and gives a solution on the cell-centered coordinate. Thus, the ®gures of ¯ow solutions are drawn using cell-centered coordinate. The ®rst example is a transonic channel ¯ow. The geometry is the channel with a 10% thick circular bump. The in¯ow Mach number is 0.675. Figure 3(a) shows the single non-adaptive grid and the solution. Figure 3(b) shows the solution-adaptive overset grid and the solution. Figure 3(c) shows the present overset-adaptive grid and the solution. It is seen from the ®gures that the resolution of the shock wave is excellent in the overset adaptive-grid solution. In Fig. 3(c), the adaptive griding is performed in one direction only. The next results are the Euler ¯ows around the BGK-1 airfoil. The freestream Mach number is 0.723, and the angle of attack is 1.598. Figure 4 shows the single grid solution. The number of grid points is 143 45. Figure 5 shows the overset grid and solution. While Fig. 6 shows the present overset adaptive-grid and its solution. The numbers of the grid points used for both the overset-grid (Fig. 5) and the overser adaptive-grid are 121 41 for the main grid, and 15 30 and 40 for the overset subgrids. The total numbers of the grid point are approximately the same for these three cases. As can be seen from the ®gures, the two subgrids are overlaid on the main grid. The adaptive griding was performed on one of the subgrids in Fig. 6. The adaptive gridings has been done using the ®rst-derivative of the pressure for the weighting functions. In this paper, the adaptive griding was applied to only one subgrid for the purpose of showing the grid quality at the shock wave. From these ®gures the present overset adaptive-grid method has been shown to give very promising results.
Overset adaptive-grid method
607
Fig. 4. Single grid (143 45) and ¯ow solution (Cp, contours, M1=0.723, a = 1.598).
8. CONCLUDING REMARKS
The overset adaptive-grid method, which combines the overset grid method with adaptive grid method, is presented in this paper. The present overset adaptive-grid method showed the promising results. The two new ideas, the use of the modi®ed distribution of the ¯ow properties and the locally-variable relaxation parameter technique, are proposed and introduced for the
608
Kenichi Matsuno et al.
Fig. 5. Overset grid (121 41, 40 25, 15 30) and ¯ow solution (Cp, contours, M1=0.723, a = 1.598).
adaptive grid method which uses elliptic equations. The numerical results show the high performance of the present method. Presently the optimum criterion for the changing locally-variable relaxation parameter is not established. The local-time-stepping-like criterion will be eective. The present locally-variable relaxation parameter method can be applied to any elliptic equations.
Overset adaptive-grid method
609
Fig. 6. Overset adaptive grid (121 41, 40 25, 15 30) and ¯ow solution (Cp, contours, M1=0.723, a = 1.598).
REFERENCES 1. Steger, J. L. and Benek, J. A., On the use of composite grid schemes in computational aerodynamics. Computer Methods in Applied Mechanics and Engineering, 1987, 64, 301±320. 2. Meakin, R. L., Moving body overset grid methods for complete aircraft tiltrotor simulation, AIAA Paper 93-3350, 1993. 3. Eiseman, P. R., Adaptive grid generation. Computer Methods in Applied Mechanics and Engineering, 1987, 64, 321± 376.
610
Kenichi Matsuno et al.
4. Fujii, K., A zonal method for practical ¯uid dynamics problems. Proceedings of the 4th International Symposium on Computational Fluid Dynamics, Davis, 1991, pp. 353±358. 5. Kao, K. H., Liou, M. S. and Chow, C. Y., Grid adaption using chimera composite overlapping meshes, AIAA Paper 93-3389, 1993. 6. Matsuno, K. and Dwyer, H. A., Adaptive methods for elliptic grid generation. Journal of Computational Physics, 1988, 77, 40±52. 7. Thompson, J. F., Warsi, Z. U. A. and Mastin, C. W., Numerical Grid Generation. Elsevier, New York, 1985. 8. Thomas, P. D. and Middleco, J. F., Direct control of the grid point distribution in meshes generated by elliptic equations. AIAA Journal, 1980, 18, 652±656. 9. Roe, P. L., Approximate Riemann solvers, parameter vectors and dierence scheme. Journal of Computational Physics, 1981, 43, 357±372.
Computers & Fluids Vol. 27, Nos 5±6, pp. 599±610, 1998 # 1998 Elsevier Science Ltd. All rights reserved Printed in Great Britain S0045-7930(97)00059-5 0045-7930/98 $19.00 + 0.00
OVERSET ADAPTIVE-GRID METHOD WITH APPLICATIONS TO COMPRESSIBLE FLOWS KENICHI MATSUNO{, MASASHI YAMAKAWA and NOBUYUKI SATOFUKA Department of Mechanical and System Engineering, Kyoto Institute of Technology Matsugasaki, Sakyo-ku, Kyoto 606 Japan (Received 6 March 1997) AbstractÐA new overset adaptive-grid method which combines a solution-adaptive grid method with an overset-grid method is proposed and applied to transonic ¯ows in this paper. The position on which the subgrid is located is adjusted to the region where the gradient of ¯ow properties is large. The solution-adaptive grid method using an elliptic equation is applied to the overset subgrid. Thus the present overset adaptive-grid method generates the grid adapted twice to the ¯ow solution. Two new ideas are also proposed for the basic adaptive-grid method of the elliptic type in this paper. One of the new ideas is to use a modi®ed distribution of ¯ow property instead of the actual one in order to estimate weighting functions which control grid spacing. The other is a locally-variable relaxation parameter technique which locally changes the relaxation parameter in the Jacobi iteration process according to the value of the coecient of the ®rst order derivatives in the elliptic equation. The present overset adaptive-grid method is applied to transonic ¯ows, and is shown to be very promising, especially for sharp shock capturing. # 1998 Elsevier Science Ltd. All rights reserved
1. INTRODUCTION
Developments of computers and numerical methods have made it possible to simulate practical ¯ow problems with reasonable accuracy. However, since the computer memory and processor speed are still limited we cannot use enough grid points for complex ¯ow problems. It is also well known that the accuracy and resolution of the numerical solutions are strongly dependent on the grid generated in the ¯ow®eld. Thus it is very important in Computational Fluid Dynamics to generate a grid of high quality. In generating a grid of high quality, there exist two important issues. One is how to eciently ®t the grid to the body boundaries of complex geometry. The overset (Chimera) grid method [1] is the one of the strategies of that problem. The overset-grid method ®rstly generates the main grid which covers the whole of the ¯ow®eld including the bodies, and secondly overlays the body-®tted subgrid on the main grid. This method has been successfully applied to various ¯ow problems in complex geometry [2]. The other important issue is how to concentrate the grid in the high-gradient regions of the ¯ow®eld. A solution-adaptive method is the one which make it possible to concentrate the grid automatically in such regions. There are two approaches of the grid adaption: moving grid and grid re®nement. The moving grid method is simpler in the data structure and easier to implement in the existing ¯ow solvers. There are many solution-adaptive methods of the moving-grid approach [3]. Among them, the adaptive methods using a system of elliptic equations are relatively popular since the elliptic equation generates smooth grid distribution due to the inherent maximum principle. The solution-adaptive grid method of elliptic type uses the source term to control the grid distribution so that the grid density is large at the large-gradient region. One of the drawbacks of the adaptive grid method of elliptic type is the grid skewness. The solution adaptive grid method generates possibly highly skewed grid, for example, in the case when the shock wave stands in the diagonal direction of the mesh lines. This often happens in the single grid structure. The overset-grid method can be used to achieve a ®ne grid locally. That is, the subgrid is overlaid on the region where the gradient of the ¯ow property is large on the main grid. This idea uses the over{Author for correspondence. 599
600
Kenichi Matsuno et al.
set-grid method in a solution-adaptive manner. A similar idea was used by Fujii [4] for tracking blast waves. Kao et al. [5] proposed the idea to use of the overset Chimera grid in large gradient regions to capture the salient features accurately during computations. These authors successfully applied the overset-grid method to ¯uid dynamics problems. However, they merely overlaid the subgrid in large gradient regions. They did not apply any moving-grid techniques to the subgrids and main grid. The purpose of this paper is to propose a new `overset adaptive-grid method'. The basic idea of this method is to apply the solution-adaptive grid method to the subgrid which is already overlaid on the high-gradient region. In the procedure, the elliptic adaptive-grid method proposed by Matsuno and Dwyer [6] is successfully applied to the overset subgrid. Two new ideas are also proposed to develop the present adaptive-grid method. One of the new ideas is to use a modi®ed distribution of ¯ow properties instead of the actual distribution in order to estimate weighting (control) functions which control grid spacing. This approach makes it possible to adapt the grid points in the high-gradient regions more smoothly and eciently. The other is a locally-variable relaxation parameter technique. The adaptive grid generation method using elliptic equations can possibly diverge when the value of control function becomes large in the iteration process. In that case, an under-relaxation is eective to prevent divergence. Thus, the locally-variable relaxation parameter method is proposed in this paper. This method is changing relaxation parameter locally from under-relaxation to over-relaxation due to the value of the control functions. These two new methods are successfully introduced in the adaptive grid generation procedures. Overset adaptive-grid solutions of transonic ¯ows are presented in this paper. Presently, the focus of adaptive-griding is on the regions of the shock waves. However, the present procedure is easily extended to more general cases. 2. ADAPTIVE-GRID METHOD
The present adaptive-grid method is based on the method proposed by Matsuno and Dwyer [6]. Here we describe the method brie¯y using one-dimensional space and after that we extend the method to two dimensions. Now we consider the one-dimensional space, whose coordinate is denoted by s. It is necessary for the adaptive-grid method to concentrate the grid in the region where the derivative of ¯ow physics is large. This property is described by the equation wDsi const:
1
where w is the weight function and Dsi=si + 1ÿsi is the grid spacing at grid index i. In this paper the following weight function, which consists of a ®rst derivative of the ¯ow property f, is adopted q 2 w 1 b fs 2 : Here, f is any ¯ow property such as pressure p, density r, Mach number M, and so on. The coecient b is a user-speci®ed free parameter, which is related to grid spacing due to the eect of ¯ow change. Equation (1) may be considered as a ®nite-dierence equation of the following dierential equation when we regard the term Dsi as Dsi=Dsi/Dx with Dx = 1. wsx const:
3
To avoid the estimation of the constant in the right-hand side of Equation (3), we dierentiate Equation (3) with respect to x. Then we get the following one-dimensional equation: wx 4 sx 0: sxx w Equation (4) is the dierential form for Equation (1). Thompson [7] developed a grid generation method using elliptic equations. This elliptic grid generation method has been applied to various types of geometry. The elliptic equation pro-
Overset adaptive-grid method
601
posed by Thomas and Middleco [8] is written in the one-dimensional form as: xss xs 2 P x:
5
Here P is a so-called `control function', which controls grid point distribution. This equation can be transformed to x coordinate by interchanging the role of dependent and independent variables. This yields the following elliptic equation for one-dimensional grid generation sxx Psx 0:
6
Comparing Equation (4) with Equation (6), the control function P relates with weight function w as follows: P wx =w:
7
Thus we get the elliptic equation for one-dimensional adaptive-grid generation sxx Psx 0, P wx =w
8
For a ¯ow-solution method, it is unfavorable from the viewpoint of solution-accuracy that the grid-spacing ratio of three consecutive points changes suddenly. Thus, it is natural that the user speci®es an allowable limit of grid-spacing ratio of the adaptive grid. Let ri be the grid-spacing ratio at index i ri Dsi =Dsiÿ1 , Dsi si1 ÿ si :
9
Now, let K r1 be the limit of the allowable grid-spacing ratio, then we get the following inequality 1=KRri RK:
10
Once we specify the allowable limit of grid-spacing ratio K, the parameter b in Equation (2) is no more free. The value of b must satisfy the inequality at grid point i 2 3 4 Kÿ1 K1 11 bR4 ÿ 2 ÿ 5: fs ÿ 4 Kÿ1 fs 2 x
K1
i
We have only one value of b. Therefore b is determined by taking the minimum of the largest bi, that is 12 b min bi : i 1, 2, . . . , imax where
2
3 4 Kÿ1 bi 4 ÿ 2 K1 ÿ 5: fs ÿ 4 Kÿ1 fs 2 x
K1
13
i
The extension of the method described above to two-dimensional is done by applying the one-dimensional formulation straightforwardly to the two arclength directions de®ned by the generalized coordinates x and Z. Let s and s' be the arclength along the Z-constant and x-constant grid lines, respectively. The basic criterion for the grid adaption is wsx const:, s2x x 2x y2x w 0 sZ 0 const:, s2Z 0 x 2Z y2Z : The two-dimensional forms of the Poisson equation [8] are ÿ xxx xyy P x, Z x2x x2y
14
602
Kenichi Matsuno et al.
ÿ Zxx Zyy Q x, Z Z2x Z2y
15
where P and Q are control functions de®ned by weight functions w and w' estimated respectively along the arclengths s(x) and s'(Z). P wx =w Q wZ 0 =w 0 :
16
Equation (15) can be inverted to give a form which can be used for numerical solution to obtain the grid ÿ ÿ 17 a rxx Prx ÿ 2brxZ g rZZ QrZ 0 where r = (x, y) and a x 2Z y2Z b x x x Z yx yZ g x 2x y2x :
18
The weight functions w and w' are de®ned the same manner as described in one-dimensional formulation. The ®rst derivative formulations along the two arclength directions are q ÿ 2 w 1 b fs w0
q ÿ 2 1 b 0 fs 0 :
19
The parameters b and b' are determined by specifying the maximum limits K and K' of the gridspacing ratio in two directions, respectively. When the control functions P and Q include the geometrical information beforehand, or if we desire to retain the eect of the solution-independent grid point distribution (usually an initial grid) determined by Pg and so on, then the eect of solution adaption may be included by de®ning control function as combination of the original-grid part Pg and adaptive-grid part Pa. The present control function, for example P, is composed of an original-grid part Pg and an adaptive-grid part Pa. P 1 ÿ yPg yPa , 0RyR1
20
where y is a parameter which controls the weight between the original and adaptive grids property. If the initial (original) grid is generated by another type of grid generation method, then we need to incorporate the information of the original grid distribution in the control functions Pg and Qg. To obtain the information of the initial grid distribution, we backsolve the control function Pi, j and Qi, j in Equation (17) as: Pgi, j
ÿ 1 ÿ yZ ÿ ax xx 2bx xZ ÿ gx ZZ ÿ x Z ÿ ayxx 2byxZ ÿ gyZZ Ja
Qgi, j
ÿ 1 ÿ yx ax xx ÿ 2bx xZ gx ZZ x x ÿ ayxx 2byxZ ÿ gyZZ Jg
21
where J = xxyZÿxZyx and derivatives of the right side are estimated at the grid point (i, j) by ®nite dierence approximations. To develop the above mentioned adaptive-grid method [6], two new techniques are proposed and introduced in this paper:
Overset adaptive-grid method
603
1. modi®ed distribution of the ¯ow properties instead of the actual distribution; 2. locally-variable relaxation parameter techniques for iterative solutions of elliptic equations.
3. MODIFIED DISTRIBUTION OF FLOW-PROPERTY
The use of the modi®ed distribution of the ¯ow properties is motivated by the consideration of making the variation of the control function smooth when the high-resolution Euler solver is applied to get the ¯ow solution. When we apply the high resolution scheme to solve inviscid (Euler) part of ¯uid dynamics equations, the shock wave is captured in two or three grid points. In this case, the control functions have large value only in a narrow region, and zero value in the other region. These control function cannot realize smooth concentration of the grid points. To avoid such situation, the present adaptive method uses a modi®ed distribution of the ¯ow variables. The modi®cation of the distribution of the ¯ow variables means in this paper to smear the obtained ¯ow solution. Suppose fi, j is an actual solution at the grid point (i, j), then the modi®ed solution f~i, j is obtained by averaging N 1 ÿfX fp , 0RfR1 f~i, j ffi, j N p1
22
where fp is a value of f at the pth grid point of N points surrounding (i, j), and f is the weight of the averaging. It is eective in getting the modi®ed smooth solution that the above averaging procedure is repeated several times. The control function Pa is estimated using the distribution of the modi®ed f~i, j .
4. LOCALLY-VARIABLE RELAXATION PARAMETER
The adaptive-grid generation Equation (17) can be solved by a relaxation method, usually the line Jacobi method, and written symbolically as Ldr ÿR r, rm1 rm o dr
23
where m is the index of iteration, and o is relaxation parameter of usually 1 Ro < 2. The under-relaxation parameter is eective for robustness in the present adaptive method when the elliptic Poisson equation is solved iteratively. It decreases the grid movement in each iteration process, especially in the case of extremely high-gradient regions such as shock wave presents. However, the under-relaxation parameter decreases the convergence rate. Hence we use the under-relaxation parameter locally and total convergence rate can be avoided decreasing. This is similar to the use of the local-time step techniques in the ¯ow solution method. In this paper, the relaxation parameter o is determined according to the local values of Pi, j and Qi, j as Ldr ÿR r, rm1 rm o i, j dr, o i, j func Pi, j , Qi, j , 0
24
25
5. OVERSET ADAPTIVE-GRID METHOD
The overset(Chimera) grid method [1] was originally developed to adjust the grid smoothly to the body surface in complex multi-body geometry by oversetting local subgrid around bodies on
604
Kenichi Matsuno et al.
the main grid covering whole of the ¯ow ®eld. The overset method can be used to achieve a ®ne grid locally. The subgrid is overlaid in the high-gradient region. This idea is to use the overset method in a solution-adaptive manner. In this paper, the following overset adaptive-grid method is proposed. First, we overlay the subgrid in the high-gradient regions. Next we apply the solution-adaptive method to the overset subgrid. Thus, the present method adapts the grid twice to the ¯ow solution. In other words, the overset-grid method is combined with the adaptive-grid method in order to get the solution adaptive grid of totally high quality and eciency. The adaptive-grid method on a single grid structure is often bothered with grid skewness. It is hard to generate orthogonal adaptive grid in multi-dimensions. However, the present overset adaptive-grid method can avoid such bad skewness by overlaying subgrid that adjusts to the ¯ow direction or gradient. The present research eort is motivated partly by consideration to avoid serious grid skewness in the adaptive griding. The present overset adaptive-grid generation procedure is as follows. The body-®tted coarse main grid is ®rstly generated by the elliptic grid generation method in the usual manner. The ¯ow is solved on the coarse main grid as a base ¯ow for overset griding and adaptive griding. The regions with large gradients of ¯ow properties, such as shock waves, are detected. The subgrids which cover the large-gradient regions are generated and put on the large-gradient regions in the main grid. In this step, each subgrid is not adapted to the ¯ow solution but merely put on the large-gradient region. Next, the solution-adaptive method using the elliptic equation described above is applied to the subgrids. Finally the ¯ow®eld is solved again on the main grid and subgrids simultaneously. The data transfer between the main grid and subgrid is done by bi-linear interpolation. To get the desirable adaptive properties, the solution-adaptive grid generation procedure is usually repeated or iterated twice or three times.
6. THE EULER EQUATION AND FLOW SOLVER
The two-dimensional compressible Euler equations written in conservation law forms are numerically solved by a cell-centered ®nite volume scheme. The scheme incorporates Roe's approximate Riemann solvers [9] for estimating ¯ux at cell interface. The time integration is performed by a diagonal ADI scheme or explicit forward Euler method. The space accuracy is increased by introducing MUSCL interpolation techniques with a minmod limiter. The present computations assume steady-state solutions. Thus, convergence acceleration techniques, for example, local time stepping, are introduced in the code.
7. APPLICATIONS
First of all, the two newly introduced techniques, the use of modi®ed ¯ow solution for the estimation of the weight function and the locally-variable relaxation parameter in the iteration, were tested using the same problem as of Ref. [6]. In this test, the adaptive solution was obtained without troubles even for the cases that the divergence were experienced in Ref. [6]. Next the test problem is the nonlinear Burgers' equation ÿ 26 ut u2 =2 x 0, 0RxR1 with the boundary condition, x = 0:u = 1 and x = 1:u = ÿ 1. The initial condition and exact steady-state solution with t 4 1 are 8 1 0RxR0:25 < 1 0RxR0:5 0:25RxR0:75 u x, 1 27 u x, 0 1 ÿ 2 xÿ0:25 0:5 ÿ1 0:5RxR1 : ÿ1 0:75RxR1
Overset adaptive-grid method
605
Fig. 1. Actual and modi®ed solution for Burgers'equation.
Note that the Roe's ¯ux-dierence splitting scheme of ®rst order accuracy captures discontinuity solution at steady state in one grid spacing in this case. The adaptive grid is generated for the steady-state solution. The modi®ed solution for estimating weight function is determined by applying the smoothing procedure, Equation (22), ten times with f = 0.5. Figure 1 shows the actual solution and its modi®ed solution. Figure 2(a) shows the adaptive grid distribution, where the weight function is estimated using actual solution. Figure 2(b) shows the adaptive grid distribution, where the modi®ed solution is used to estimate the weight function. The ®rst derivative is used to estimate the weighting function for adaptive griding. From these ®gures, it
Fig. 2. Adaptive grid for the Burgers' equation (21 points): (a) Pa by actual solution; and (b) Pa by modi®ed solution.
606
Kenichi Matsuno et al.
Fig. 3. Transonic ¯ow through channel with 108 thick circlar arc. Cp, Contours, Min¯ow=0.675. (a) Single grid solution (70 24), (b) overset grid solution (49 17 for the main grid, 20 17, 16 16, 16 16 for the subgrids), (c) overset adaptive-grid solution [the numbers of grid points are the same as (b)].
is shown that the use of modi®ed solution for estimation of weight function is very promising for adaptive griding. Next we show the results of the overset adaptive-grid method applied to transonic ¯ows. The present Euler solver is the cell-centered ®nite volume scheme and gives a solution on the cell-centered coordinate. Thus, the ®gures of ¯ow solutions are drawn using cell-centered coordinate. The ®rst example is a transonic channel ¯ow. The geometry is the channel with a 10% thick circular bump. The in¯ow Mach number is 0.675. Figure 3(a) shows the single non-adaptive grid and the solution. Figure 3(b) shows the solution-adaptive overset grid and the solution. Figure 3(c) shows the present overset-adaptive grid and the solution. It is seen from the ®gures that the resolution of the shock wave is excellent in the overset adaptive-grid solution. In Fig. 3(c), the adaptive griding is performed in one direction only. The next results are the Euler ¯ows around the BGK-1 airfoil. The freestream Mach number is 0.723, and the angle of attack is 1.598. Figure 4 shows the single grid solution. The number of grid points is 143 45. Figure 5 shows the overset grid and solution. While Fig. 6 shows the present overset adaptive-grid and its solution. The numbers of the grid points used for both the overset-grid (Fig. 5) and the overser adaptive-grid are 121 41 for the main grid, and 15 30 and 40 for the overset subgrids. The total numbers of the grid point are approximately the same for these three cases. As can be seen from the ®gures, the two subgrids are overlaid on the main grid. The adaptive griding was performed on one of the subgrids in Fig. 6. The adaptive gridings has been done using the ®rst-derivative of the pressure for the weighting functions. In this paper, the adaptive griding was applied to only one subgrid for the purpose of showing the grid quality at the shock wave. From these ®gures the present overset adaptive-grid method has been shown to give very promising results.
Overset adaptive-grid method
607
Fig. 4. Single grid (143 45) and ¯ow solution (Cp, contours, M1=0.723, a = 1.598).
8. CONCLUDING REMARKS
The overset adaptive-grid method, which combines the overset grid method with adaptive grid method, is presented in this paper. The present overset adaptive-grid method showed the promising results. The two new ideas, the use of the modi®ed distribution of the ¯ow properties and the locally-variable relaxation parameter technique, are proposed and introduced for the
608
Kenichi Matsuno et al.
Fig. 5. Overset grid (121 41, 40 25, 15 30) and ¯ow solution (Cp, contours, M1=0.723, a = 1.598).
adaptive grid method which uses elliptic equations. The numerical results show the high performance of the present method. Presently the optimum criterion for the changing locally-variable relaxation parameter is not established. The local-time-stepping-like criterion will be eective. The present locally-variable relaxation parameter method can be applied to any elliptic equations.
Overset adaptive-grid method
609
Fig. 6. Overset adaptive grid (121 41, 40 25, 15 30) and ¯ow solution (Cp, contours, M1=0.723, a = 1.598).
REFERENCES 1. Steger, J. L. and Benek, J. A., On the use of composite grid schemes in computational aerodynamics. Computer Methods in Applied Mechanics and Engineering, 1987, 64, 301±320. 2. Meakin, R. L., Moving body overset grid methods for complete aircraft tiltrotor simulation, AIAA Paper 93-3350, 1993. 3. Eiseman, P. R., Adaptive grid generation. Computer Methods in Applied Mechanics and Engineering, 1987, 64, 321± 376.
610
Kenichi Matsuno et al.
4. Fujii, K., A zonal method for practical ¯uid dynamics problems. Proceedings of the 4th International Symposium on Computational Fluid Dynamics, Davis, 1991, pp. 353±358. 5. Kao, K. H., Liou, M. S. and Chow, C. Y., Grid adaption using chimera composite overlapping meshes, AIAA Paper 93-3389, 1993. 6. Matsuno, K. and Dwyer, H. A., Adaptive methods for elliptic grid generation. Journal of Computational Physics, 1988, 77, 40±52. 7. Thompson, J. F., Warsi, Z. U. A. and Mastin, C. W., Numerical Grid Generation. Elsevier, New York, 1985. 8. Thomas, P. D. and Middleco, J. F., Direct control of the grid point distribution in meshes generated by elliptic equations. AIAA Journal, 1980, 18, 652±656. 9. Roe, P. L., Approximate Riemann solvers, parameter vectors and dierence scheme. Journal of Computational Physics, 1981, 43, 357±372.