Role of Coordinates in Computational Fluid Dynamics
W. H. Hui1, J. J. Hu2and Keh-Ming Shyue3
1Division of Mechanics, Research Centre for Applied Sciences Academia Sinica, Taipei, Taiwan
2Department of Information Management, Shu-Te University, Kaohsiung County, Taiwan
3Department of Mathematics, National Taiwan University, Taipei, Taiwan Email: whhui@ust.hk
A
BSTRACTComputational fluid dynamics uses large scale numerical computation to solve problems of fluid flow. It turns out that the numerical solution for a given flow depends on the coordinates (grid) used to compute the flow. The commonly used Eulerian and Lagrangian coordinate systems both have advantages and drawbacks. In this paper we first discuss the role of coordinates in computational fluid dynamics regarding the questions of: (a) conservation form partial differential equations, (b) numerical resolution of contact discontinuities, (c) grid generation, and (d) grid orthogonality. We then introduce a unified coordinate system which combines the advantages of both Eulerian and Lagrangian system and beyond, while avoiding their drawbacks. Examples include a transonic flow past an airfoil and a two-fluids flow with shocks.
Keywords: Unified coordinates; Automatic Grid Generation; Multi-Fluids flow
1. Eulerian versus Lagrangian Coordinates
The starting point in CFD is the physical conservation law
) ( )
(tEd tF ndS
t
(1)
In Eulerian coordinates, which are fixed in space, Eq. (1) can be written as partial differential equations (PDE) in conservation form. Thus, for a
-law perfect gas we have0
y G x F t
Ee e e (2)
. )]
( , , , [
, )]
( , , , [
, ] , , , [
2 2
T e
T e
T e
p e v p v uv v G
p e u uv p u u F
e v u E
In (2), t is time, p pressure,density, u and v are the x- and y-component of fluid velocity, and
v p u
e 1
) 1 2(
1 2 2
.
Shock-capturing computation is therefore easy to apply in the Eulerian system, but it tends to smear contact discontinuities badly. Furthermore, for computing flow past a body, which is a central problem in CFD, it is necessary to generate a body-fitted grid prior to flow computation.
By contrast, in Lagrangian coordinates, contact discontinuities can be resolved sharply.
However, the coordinates are flow dependent, hence Eq. (1) cannot easily be written in conservation PDE form (except in the special case of one-dimensional flow). This complicates the computation, for instance, a stagger grid is needed, causing numerical diffusion. Moreover, the grid deforms with the fluid, causing inaccuracy and even breakdown of computation.
2. The “ Opt i mal ” Coordinate System
Can we have a coordinate system that combines the advantages of the Eulerian and Lagrangian, while avoiding their drawbacks? Such a system would be “optimal” in some sense(whetheror not a system is optimal depends on the criteria, which are necessarily subjective). Specifically, we want the coordinate system to possess the following properties:(a) Conservation PDE form exists, as in Eulerian;
(b) Contact discontinuities are sharply resolved, as in Lagrangian;
(c) Grid can be automatically generated to fit given body shapes;
(d) Grid is orthogonal;
(e) Grid is uniform;
(f)……
We shall show that the unified coordinate system introduced in the next section possesses these properties.
3. The Unified Coordinate System
For simplicity, we consider 2-D flow and introduce a unified coordinate system (λ,ξ,η)[1, 2] via the transformation from Eulerian system (t, x, y)(3) From (3), we get
0
Dt
DQ (4)
where
t x
Dt
D
Q
Q
So, the coordinates (
,
), and hence computational cells, move with velocity Q = (U, V). It includes the Eulerian as a special case when Q = 0 and Lagrangian when Q = q = (u, v).Under the transformation, (2) becomes
0
G F
E (5)
where
0 0
) , (
0 0
) , (
V U
uB vA p Ye
pA Yv
pB Yu
Y
G V
U vL uM p Xe
pL Xv
pM Xu
X
F
M L B A Je Jv Ju J
E
Here J = AM–BL, X = (u - U)M–(v - V)L and Y = (v – V)A – (u – U)B. The first four equations in (4) are the physical conservation laws. But since they involve the coefficients A, B, L, M of the transformation, they are not closed.
It is, therefore, necessary and sufficient to append the time evolution of these coefficients, i.e., the last four equations, to make the system closed. These four equations are the compatibility conditions of the transformation and are called geometric conservation laws.
4. Grid Movement
Since there are two arbitrary functions, U and V, we prescribe two requirements:
(A) Coordinate lines= const. shall be material
lines of fluid particles, meaning
0 Dt
Dq (6)
Together with the first of (4), we get
(v - V)A = (u - U)B (7) Observations:
(a) Contact lines, being material lines, must coincide with coordinate lines and, therefore, can be resolved sharply.
(b) As the body surface is a material line, condition (7) guarantees that the unified grid is automatically a body-fitted grid at all time. This provides the basis for automatic grid generation.
(c) A material interface (including a free surface) corresponds to = const and thus can be resolved sharply.
(d)(x, y, t) is a level set function; hence there is no need to introduce an extra function when using the level set method.
All these make the unified coordinate approach particularly suitable and useful for problems of multi-fluids flow and flow-solid interactions, etc.
(B) Grid angles and hence grid orthogonality shall be preserved during the
marching computation. Thus we have0 cos
cos
2 2 2 2 1 1
M L B A
BM AL
(8)
After eliminating V from (7) and using the geometric conservation laws, (8) becomes an ordinary differential equation for U
) ,
; ( ) ,
;
(
P U Q
U
(9)
2 2 2 2 2 2
2 2 2 2
,
) ,
; ( ) ,
; (
) ,
; (
B A T M L S
uP
B u A v J L A v B u J T
A Q S
B A A B AJ
L B A
A B J T P S
We can specify any initial data for U at =
const.
5. Computation Procedure
The special case of steady supersonic flow was successfully studied in [3], and here we are concerned with steady flow which may have subsonic regions. Flows of this type are computed by marching in time
until a steady state is reached.
Md Bd Vd dy
Ld Ad Ud dx
d dt
The computation procedure for uniform flow past a body is illustrated by a Mach 0.8 steady air flow (
= 1.4) past a NACA 0012 airfoil as follows:Initialization Stage --- automatic generation of body-fitted grid in a computational window.
Given the grid sizes, x andy , and the number of cells, MN, in the window (we use M = 200 and N = 100 in the example).
Step 1. Begin with a column of N orthogonal cells, representing the given uniform flow in the x-direction (Fig. 1a). This gives the initial vales of (A, B, L, M) = (1, 0, 0, 1). We also take (U, V) = (u, v) initially.
Step 2. Compute the solution to Eq. (5) by marching in timeusing dimensional splitting: splitting into two 1-D systems inand each of them is solved using the standard Godunov/MUSCL scheme with the minmod limiter.
(Details are as follows. To update the solution from time n to time n+1: (a) Solve the first four equations (the physical conservation laws) of Eq. (5) for (ρ, p, u, v) keeping A, B, L, M, U and V at time n level. (b) Use this updated values together with a specified initial data to solve Eq. (9) for U and then Eq. (7) for V at time n+1. (c) Use these updated values of U and V to update (A, B, L, M) at time level n+1 by integrating the geometric conservation laws. At all outer boundaries of the computational region at every time step, we apply the characteristic boundary conditions). After one time step,this column of cells moves to the right by U.
Step 3.After several time steps when the initial column of cells has moved to the right by a distance equal tox, add one new column of cells on the left that is identical to the initial column.
Step 4. Repeat this process of adding cell columns on the left of the computational region until the leading column meets the body surface (Fig. 1b) we then impose the boundary condition of zero normal velocity on the body surface.
Step 5. Continue this process until after the columns of cells cover the whole body surface and further downstream, when we have M columns of cells in the window (Figs. 1c, 1d and 1e).
This completes the initialization stage, and we now have a body-fitted grid (Figs. 1e and 1f) and a flow field around the airfoil in the window.
The computed grid is seen orthogonal, as predicted. It is also fairly uniform in the x-direction, because in solving Eq. (9) we have specified uniform data for U at = const. The associated flow field computed (e.g., surface pressure in Fig. 1g) is, however, only a very rough approximation to the correct one (see, e.g., the potential solution of Hafez et al [4]), partly because it has not reached the steady state and partly because the downstream boundary conditions used at the transient times, e.g., in Figs. 1a to 1d, are obviously incorrect as the computational regions at those times are not the full window.
To progress further, one could use the body-fitted orthogonal grid generated so far to perform an Eulerian computation with the associated flow field as an initial solution. This can be easily done by putting U = V = 0 (without solving Eq. (9)) during the subsequent iterations towards a steady state. In this way, the unified coordinate approach plays the role of grid generation for Eulerian computation.
An alternative and better way is to continue the unified coordinate computation to proceed to the Main Stage --- iteration with flow-adjusted grids.
Step 6. To iterate the solution towards a steady state, whenever we add a new column of cells on the left of the window we also simultaneously delete the right-most column of cells from the computation window, thus keeping the window in the same size. At the same time we improve the solution by using the information of the flow field at every time step, e.g., the surface pressure gradient, to adjust the initial data of U at= const in solving (9) so that the grid is refined in regions of high pressure gradient (Fig. 1h). We note that this flow-adjusted refined grid remains orthogonal. The computed surface pressure distribution is shown in Fig. 1i, which is much better than that in Fig. 1g and is in good agreement with the potential flow computation of Hafez et al [4]. The pressure contours at the same time is shown in Fig 1j.
The flow-adjusted grid in Fig. 1h looks similar to those obtained by the grid re-distribution method in Euelrian computation.
But there are differences: grid re-distribution requires generating another grid at every time step by solving an elliptic equation, and conservation properties in the interpolations of
geometry and flow variables between the two grids must be ensured. These issues do not arise in our flow-adjusted grid approach because we need only one grid; the only modification is to specify the initial data for U at = const by using the pressure gradient information.
In terms of computing time, most of the CPU time is used in solving the Riemann problems for the physical conservation laws, which is the same as in Eulerian computation. However, additional times are needed to solve Eq. (9) and to update the geometric variables (A, B, L, M).
These typically increase the CPU time by 5-10%.
(This can also be reduced if in step 6 after the flow-adjusted grid is well established we put U
= V = 0 so that there is no need to solve Eq. (9) subsequently). On the other hand, grid re-distribution in Eulerian computation by solving an elliptic equation increases CPU time.
What is most important is that Eulerian computation always needs generating a body-fitted grid prior to flow computation, and this can be time-consuming.
6. Examples
Example 1 showing a transonic flow computation is already given above. Example 2 shows a sample computation for a Mach 2.2 steady flow over a NACA 0012 airfoil at an angle of attack of 8 degrees and a horizontal air-SF6 material. Figs. 2a to 2c show the flow-generated grids at different times, whereas Figs. 2d to 2h show the computed contours of density, pressure and entropy at those times. The interface between air and SF6 was initially identified with a particular value of and it remained so, because a contact line coincides with a coordinate line = const in the unified coordinate system. The computation was straight forward, with no special treatment; yet the results are better than the corresponding Eulerian ones (Figs. 3a to 3d), in particular the interface and the slip line behind the airfoil are resolved sharper.
7. Conclusion
This paper demonstrates that coordinate system plays an important role in computational fluid dynamics. It further shows that the unified coordinate system is superior to the Eulerian and the Lagrangian system.
References
[1] W. H. Hui, P. Y. Li, and Z. W. Li., “A unified coordinate system for solving the two-dimensional Euler equations”. J.
Comput. Phys., 153:596–637, 1999.
[2] W. H. Hui, and S. Kudriakov, “A unified coordinate system for solving the two-dimensional Euler equations”. J.
Comput. Phys., 172: 235-260, 2001.
[3] Hui, W.H. and Hu, J.J., “Space-marching gridless computation of steady supersonic/hypersonic flow”, Int’l J.
Comput. Fluid Dynamics, (in print) 2006.
[4] M. M. Hafez, S. Osher,, and W. Whitlow,
“Improved finite different schemes for transonic potential calculations”. AIAA Paper 84-0092, AIAA 22nd Aerospace Sciences Meeting. 1984.
(a) Initial Grid at t=0.0
(b) Flow-Generated Grid at t=4.0
(c) Flow-Generated Grid at t=4.6
(d) Flow-Generated Grid at t=4.9
(e) Flow-Generated Grid at t=10.1
(f) Blow-Up Grid at t=10.1
(g) Surface Pressure at t=10.1
X -CP
0 0.25 0.5 0.75 1
-1 -0.5 0 0.5 1
Pressure at t=10.1 Hafez et al. [4]
(h) Flow-Adjusted Blow-Up Grid at t=50.0
(g) Surface Pressure at t=50.0
X -Cp
0 0.25 0.5 0.75 1
-1 -0.5 0 0.5 1
Present Hafez et al. [4]
(h) Pressure Contours at t=50.0
Figure 1. Computed results for a Mach 0.8 steady flow past a NACA 0012 airfoil, showing flow-generated grids at different times, surface pressure and pressure contours.
= 1.4(b) Flow-Generated Grid, t=1.04 (a) Flow-Generated Grid, t=0.128
(c) Flow-Generated Grid, t=4.0
Figure 2. Sample computation of a Mach 2.2 steady flow past a NACA 0012 airfoil at angle of attack of 8 degrees and a horizontal air-SF6
material interface. The flow-generated grid, and the contours of the density, pressure and entropy are shown at three different times: (a) t = 0.128, (b) t = 1.04, and (c) t = 4.0.
(d) Density Contours, t=0.128
(e) Density Contours, t=1.04
(f) Density Contours, t=4.0
(g) Pressure Contours, t=4.0
(h) Entropy Contours, t=4.0
Figure 3. Eulerian computation of a Mach 2.2 steady flow past a NACA 0012 airfoil at angle of attack of 8 degrees and a horizontal air-SF6
material interface.
(a) Local Eulerian Grid
(b) Density Contours
(c) Pressure Contours
(d) Entropy Contours
