Differential Equations and Computational Simulations III J. Graef, R. Shivaji, B. Soni, & J. Zhu (Editors) Electronic Journal of Differential Equations, Conference 01, 1997, pp. 149–157. ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu ftp ejde.math.swt.edu (login: ftp)

Governing Equations of Fluid Mechanics in Physical Curvilinear Coordinate System ^{∗}

Swungho Lee & Bharat K. Soni

Abstract This paper presents the development of unsteady three-dimensional incompressible Navier-Stokes and Reynolds-averaged Navier-Stokes equations in an unsteady physical curvilinear coordinate system. It is demonstrated that the numerical simulations based on governing equations in a physical curvilinear coordinate system are less mesh sensitive than other schemes.

Introduction

The ultimate goal of the current research is to develop a numerical flow simulation system applicable to unsteady three dimensional incompressible Navier-Stokes equations that is accurate and efficient in view of CPU time taken for convergence. Patel et al. [1] note that,

If the geometry is highly curved and skewness of angles between the velocity components and the faces of the computational cells are large, an approach that transforms only the independent coordinate variables in the equations representing the transport of mass and momentum may lead to increased numerical diffusion.

This fact provides the motivation for the development of governing equations in physical curvilinear component form. In this coordinate system, components of the velocity have the same direction as that of the coordinate lines and have physical values. The physical curvilinear-components form of the velocity was first introduced by Truesdell [2]. Demirdzic et al. [3] derived the physical curvilinear-components form in nonorthogonal coordinates for Reynoldsaveraged Navier-Stokes equations and the equations of the k − ϵ turbulence model. In their derivation, the equations of the Cartesian tensor forms were transformed directly into physical curvilinear-component forms by a two step procedure. Takizawa et al. [4] used this form to simulate a two-dimensional free-surface problem using a different concept for the connection coefficients. However, in this approach, practical applications have been limited to twodimensional problems because of the large storage requirements of geometric tensors and connection coefficients, as well as the numerical error associated with the evaluation of these geometric tensors and connection coefficients.

∗1991 Mathematics Subject Classifications: 65N30, 35Q30.

Key words and phrases: Physical curvilinear coordinate systems,

incompressible flows, transformations.

\begin{matrix}c \\ ◯ \end{matrix}1998 Southwest Texas State University and University of North Texas.

Published November 12, 1998.

This work explores an approach which is different from Demirdzic et al [3] in that the partial differential equations with the coordinate-free vector form are transformed into the physical curvilinear coordinate system using general transformation laws. The resulting unsteady three-dimensional incompressible viscous equations based on the unsteady physical curvilinear coordinate system derived in this research are validated by numerically simulating the laminar three-dimensional lid-driven cavity and the free surface turbulent flows.It is demonstrated that these numerical simulations are less mesh sensitive.

Basic transformations

In the following analysis, x_{i} are Cartesian coordinates, ξ^{i} are curvilinear coordinates and ξ^{(i)} are physical curvilinear coordinates. First, consider the general transformation law under the two changes of the coordinates system x_{i} to ξ^{i} and ξ^{i} to ξ^{(i)}. The relationship between the coordinates x_{i} and the coordinates ξ^{i} can be expressed as follows:

ξ^{i} = ξ^{i}(x_{j}, t)  (1)

The relationship between the coordinates ξ^{i} and the coordinates ξ^{(i)} can be expressed as:

ξ^{(i)} = ξ^{(i)}(ξ^{i}) , where Δξ^{(i)} = \sqrt{g_{ii}}Δξ^{i}  (2)

\sqrt{g_{ii}} are evaluated at ξ^{k} = constant and k ≠ i. ξ^{(i)} resemble the coordinate stretching in each direction of ξ^{i}. In view of transforming the coordinates from x_{i} to ξ^{i} and ξ^{i} to ξ^{(i)}, the vector dr⃗ can be written as:

dr⃗ = \frac{∂r⃗}{∂ξ^{i}} dξ^{i} = a⃗_{i}dξ^{i}  (sum on i)  (3)
= \frac{∂r⃗}{∂ξ^{(i)}} dξ^{(i)} = a⃗_{(i)}dξ^{(i)} , where a⃗_{(i)} = \frac{1}{\sqrt{g_{ii}}} a⃗_{i}  (4)

a⃗_{i} are covariant base vectors in the curvilinear coordinate system and a⃗_{(j)} are covariant base vectors in the physical curvilinear coordinate system. The relationships between each coordinate system for each covariant and contravariant metric tensors are written as:

a⃗_{(i)} · a⃗_{(j)} = g_{(ij)} = \frac{1}{\sqrt{g_{ii}}\sqrt{g_{jj}}} g_{ij} , where a⃗_{i} · a⃗_{j} = g_{ij}  (5)
a^{(i)} · a^{(j)} = g^{(ij)} = \sqrt{g_{ii}}\sqrt{g_{jj}}g^{ij} , where a^{i} · a^{j} = g^{ij} , a^{(i)} = \sqrt{g_{ii}}a^{i} (6)

g_{(ij)} and g^{(ij)} are the physical covariant and the physical contravariant metric tensors, respectively. The physical covariant metric tensor g_{(ij)} are equal to one if the subscripts i and j are the same. To obtain the divergence, gradient and Laplacian operators of a vector in the physical curvilinear coordinate system, one starts from the covariant and the contravariant derivatives of the base vectors. Before obtaining the divergence, however, one needs to define the physical curvilinear components of the velocity vector. These components can be defined as the magnitude of the i^{th} component projected onto the i^{th} physical curvilinear coordinate direction,

u^{(i)} = u⃗ · a⃗_{(i)} .  (7)

Here u^{(i)} represents the physical curvilinear component of the velocity vector and has physical values. In the Cartesian coordinate system, u^{(i)} are identical to the physical components of the velocity u(i).

The derivatives of the covariant base vector in the physical curvilinear coordinate system are obtained by taking the derivative of the covariant base vector in the curvilinear coordinate system.

\frac{∂a⃗_{(i)}}{∂ξ^{(j)}} = [\begin{matrix}√ & \overline{g_{kk}} \\ \sqrt{g_{ii}} & \sqrt{g_{jj}} \end{matrix} Γ_{ij}^{k} − δ_{i}^{k} \frac{g_{km}}{g_{ii}\sqrt{g_{jj}}} Γ_{ij}^{m}]a⃗_{(k)} (sum on k and m)
= Γ_{(ij)}^{(k)}a⃗_{(k)} (sum on k) ,  (8)

where

Γ_{(ij)}^{(k)} = \begin{matrix}√ & \overline{g_{kk}} \\ \sqrt{g_{ii}} & \sqrt{g_{jj}} \end{matrix}Γ_{ij}^{k} − δ_{i}^{k} \frac{g_{km}}{g_{ii}\sqrt{g_{jj}}}Γ_{ij}^{m} .

Similarly, the derivatives of the contravariant base vector in the physical curvilinear coordinate system are obtained by taking the derivative of the contravariant base vector.

\frac{∂a^{(i)}}{∂ξ^{(j)}} = −Γ_{(kj)}^{(i)}a^{(k)} (sum on k)  (9)

Here the Christoffel and physical counterparts of the Christoffel symbols have the following properties:

Γ_{jk}^{i} = Γ_{kj}^{i} and Γ_{(jk)}^{(i)} ≠ Γ_{(kj)}^{(i)}  (10)

Using the general transformation laws for a scalar ϕ, the gradient can be written as follows:

∇ϕ = \frac{∂ϕ}{∂ξ^{(i)}} a^{(i)} = g^{(ik)} \frac{∂ϕ}{∂ξ^{(k)}} a⃗_{(i)} (sum on i and k) ,  (11)

where a^{(i)} = g^{(ik)}a⃗_{(k)}. Also, the gradient of a vector u⃗ can be expressed using equation as:

∇u⃗ = \frac{∂u⃗}{∂ξ^{(i)}} a^{(i)} = \frac{∂(u^{(k)}a⃗_{(k)})}{∂ξ^{(i)}}a^{(i)} (sum on i and k)
= u_{(,i)}^{(k)}a⃗_{(k)}a^{(i)} = g^{(ij)}u_{(,i)}^{(k)}a⃗_{(k)}a⃗_{(j)} ,

where

u_{(,i)}^{(k)} = \frac{∂u^{(k)}}{∂ξ^{(i)}} + u^{(m)}Γ_{(im)}^{(k)} .

The quantity u_{(,j)}^{(k)} is called the covariant derivative of the physical curvilinear components of a vector u⃗. One can easily evaluate the divergence of a vector u⃗ using equation (8) as:

∇ · u⃗ = \frac{∂u⃗}{∂ξ^{(i)}} · a^{(i)} = \frac{∂u^{(k)}a⃗_{(k)}}{∂ξ^{(i)}} · a^{(i)} (sum on i and k)
= u_{(,i)}^{(i)} = \frac{\sqrt{g_{ii}}}{J} \frac{∂}{∂ξ^{(i)}} (\frac{J}{\sqrt{g_{ii}}} u^{(i)}) (sum on i)  (12)

The Laplacian can be evaluated by the divergence of the gradient of a vector u⃗,

∇^{2}u⃗ = ∇ · T̂ = \frac{∂T̂}{∂ξ^{(j)}} · a^{(j)} , where T̂ = ∇u⃗ = u_{(i,)}^{(k)}a⃗_{(k)}a^{(i)}

= g^{(jk)}[u_{(,j)}^{(m)}Γ_{(mk)}^{(i)} − u_{(,m)}^{(i)}Γ_{(jk)}^{(m)} + \frac{∂u_{(,j)}^{(i)}}{∂ξ^{(k)}} ]a⃗_{(i)} (sum on i, j, k and m)

The time derivative of a scalar or a vector F is given as Warsi [5]:

\frac{∂F}{∂t} |_{x_{i}} = [\frac{∂F}{∂τ} + \frac{∂F}{∂ξ^{(i)}} \frac{∂ξ^{(i)}}{∂t}]|_{ξ^{(}i)}  (13)
\frac{∂F}{∂τ} |_{ξ^{(}i)} = [\frac{∂F}{∂t} + \frac{∂F}{∂x_{i}} \frac{∂x_{i}}{∂τ}]|_{x_{i}}  (14)

Here τ represents the time in an unsteady physical curvilinear coordinate system. The grid speed can be easily evaluated by replacing F with ξ^{(i)} in equation (14).

w^{(i)} = \frac{∂ξ^{(i)}}{∂t} = −\frac{∂ξ^{(i)}}{∂x_{l}} \frac{∂x_{l}}{∂τ} = −\frac{\sqrt{g}_{ii}}{J} b_{l}^{i} \frac{∂x_{l}}{∂τ} (sum on l)  (15)

where

b_{l}^{i} = \frac{∂x_{m}}{∂ξ^{j}} \frac{∂x_{n}}{∂ξ^{k}} − \frac{∂x_{n}}{∂ξ_{j}} \frac{∂x_{m}}{∂ξ^{k}}

with i, j, k, and l, m, n in each cyclic order.

The divergence of the grid speed vector w⃗ is written as:

∇ · w⃗ = \frac{\sqrt{g}_{ii}}{J} \frac{∂(\frac{J}{\sqrt{g}_{ii}} w^{(i)})}{∂ξ^{(i)}} = −\frac{l}{J} \frac{∂J}{∂τ}

Governing equations in the unsteady physical curvilinear coordinate system

By replacing F with u⃗ in equation (13), one can get the equations into an unsteady coordinate system. The vector form of incompressible Reynolds-averaged Navier-Stokes equations, with the body force in unsteady coordinates system, is given by equation (16).

\frac{∂u⃗}{∂τ} + (∇u⃗) · v⃗ = −∇P + ν_{E}∇^{2}u⃗ + [∇u⃗ + (∇u⃗)^{T}] · ∇ν_{E} ,  (16)

where

v⃗ = u⃗ + w⃗, P = ρ + \frac{z}{F n^{2}} + \frac{2}{3}k and ν_{E} = \frac{l}{R_{eff}} = \frac{l}{R_{e}} + ν_{t} .

Here w⃗ is a grid speed vector, F n is a Froude number, P is total pressure, and ρ is a static pressure. ν_{t} and R_{eff} represent the eddy viscosity and the effective Reynolds number, respectively.

The procedure for the transformation of the incompressible Reynolds-averaged Navier-Stokes equations, based on an unsteady physical curvilinear coordinate system, is now introduced using the derivations for the gradient, divergence operator, Laplacian, and time derivative. An unsteady physical curvilinear component form of the continuity and the Reynolds-averaged Navier-Stokes equations can be presented as:

\frac{∂}{∂ξ^{i}} (\frac{J}{\sqrt{g}_{ii}} u^{(i)}) = 0

R_{eff} \frac{∂u^{(i)}}{∂τ} + R_{eff}v^{(j)}(\frac{∂u^{(i)}}{∂ξ^{(j)}} + u^{(k)}Γ_{(kj)}^{(i)})  (17)
= −R_{eff} [g^{(ij)} \frac{∂P}{∂ξ^{(j)}} − \frac{∂ν_{E}}{∂ξ^{(j)}} [g^{(jk)} \frac{∂u^{(i)}}{∂ξ^{k}} + g^{(jk)}Γ_{(lk)}^{(i)}u^{(l)} + g^{(ik)} \frac{∂u^{(i)}}{∂ξ^{k}}

+g^{(ik)}Γ_{(lk)}^{(j)}u^{(l)}] ] + g^{(jk)}[\frac{∂^{2}u^{(i)}}{∂ξ^{(k)}∂ξ^{(j)}} + \frac{∂(u^{(l)}Γ_{(lj)}^{(i)})}{∂ξ^{(k)}} + Γ_{(lk)}^{(i)}u_{,(j)}^{(l)} − Γ_{(jk)}^{(l)}u_{,(l)}^{(i)}]

Equation (17), which has a nonconservative form, is rearranged into the standard form for the use of the finite analytic method [6]. Using the stretched coordinates ξ^{i}∗, the 12-point finite analytic discretization scheme based on the local nonuniform grid spacing [6] is employed. The stretched coordinates ξ^{i}∗ are defined as:

ξ^{i} = \sqrt{g^{ii}}ξ^{i} ∗ or ξ^{i}∗ = \frac{l}{\sqrt{g^{ii}}} ξ^{i}

The first equation shows the relation between the curvilinear coordinates and the stretched coordinates. The \sqrt{g^{ii}} are evaluated at ξ^{k} = constant and k ≠ i.

Results and discussions

An unsteady three-dimensional incompressible flow solver based on the physical curvilinear coordinate system has been developed [6]. The 12-point finite analytic scheme with enhanced kinematic boundary condition and numerical approach was utilized in this development. The detailed discussions on the numerical scheme can be found in [6]. The results of the following two test cases to validate the pertinent numerical simulations are presented here.

Figure 1: Convergence Histories for Different Time Steps

Figure 2: Comparison of the Centerline u-Velocity Profile

Lid-driven three-dimensional cavity flow

The velocity components on the wall are zero, except on the moving wall with a velocity of 1. The computations are performed on a grid consisting of 16×16×16 grid points. The dimensionless time step is taken from 0.01 to 2. The matrix that consists of the coefficients resulting from the finite analytic method was solved using the GMRES (Generalized Minimal RESidual) method [9]. To obtain the solution for the steady state, only one iteration per time step is used. All computations are performed using a relaxation factor of 1. Figure 1 shows that the rate of the convergence depends on the size of the time step in the range from 0.01 to 2. The computations lead to a fully converged solution within fewer than 200 iterations.

Figure 2 shows the comparison of the center line u-velocity profiles at the different Reynolds numbers. The computations were performed on both uniform and nonuniform grids. The time increment is set to Δt = 1.

Peric [7] has reported that if the angle between the two coordinate lines is greater than 135^{∘} or less than 45^{∘}, then the pressure correction equation does not converge at all, or the convergence rate is too slow. Cho and Chung [8] used a new treatment method for nonorthogonal terms in the pressure-correction equation in order to enlarge the ranges for convergence and found that the smaller

Figure 3: Convergence Histories in Various Inclined Angles (Degree)

the angle, the narrower the region of relaxation factor. In the present research, the computations are performed at several inclined angles (90^{∘}, 60^{∘}, 45^{∘}, 30^{∘}, 15^{∘}, 10^{∘}, and5^{∘}) to check the influence on the rate of the convergence due to the grid skewness. As mentioned, all computations of the threedimensional cavity flow are performed using a relaxation factor of 1. The solution always converged, even for very small inclined angles, but more iterations were required for convergence for small inclined angles. Figure 3 shows the convergence histories in various inclined angles from 90 to 5 degrees.

Ship flow with the free surface

It has been shown that the present code is less mesh sensitive and converges well even at the large grid skewness for the three-dimensional cavity flow. For the next case, three-dimensional moving free surface turbulent flow was simulated. The upper boundaries move arbitrarily with the flow, and the grid in the computational domain is generated at every time step until the solution of the steady state is obtained.

Table 1: Grid Dependence Tests for the Wigley Hull

I  II  III
Grid Points  125 × 35 × 34 125 × 40 × 40 125 × 50 × 48
Total Nodes  148,750  200,000  300,000
Time Increment  0.005  0.005  0.005
Total time steps  400  400  400
Total CPU (hours)  34.91  46.56  60.68
Reynolds Number  1.0 × 10^{6}  1.0 × 10^{6}  1.0 × 10^{6}
Froude Number  0.289  0.289  0.289

Table 1 shows the information for these computations. The computations Figure 4: The Wave Elevations on the Hull Surface for different grid sizes and a Perspective view of the wave elevation

were performed with three different grids under the same conditions. The Baldwin-Lomax turbulence model [9] was used to calculate the eddy viscosity in the turbulent flow. The Froude number and the Reynolds number used in the experiment [10] are 0.289 and 3.3 × 10^{6}, respectively. The wave elevations on the hull surface and a perspective view of the wave elevation is shown in Figure 4. A deviation of wave profiles is observed in the bow region, while better agreement is seen toward the stern. The bow peak is not captured properly due to the large spacing of (Δx) in a region of relatively high gradient. The residue remains around 10^{−4}after t = 0.5.

The comparison of these numerical simulations with the results reported in open literature have shown very good agreement. It is demonstrated, especially in the cavity flow simulation, that the numerical simulations involving a physical curvilinear coordinate system are less mesh sensitive.

References

[1] Patel, V. C., Chen, H. C. and Ju, S., “Solutions of the Fully-Elliptic Reynolds-Averaged Navier-Stokes Equations and Comparisons with Experiments,” IIHR Report Number 323, The University of Iowa, 1988.

[2] Truesdell, C., “The Physical components of vectors and Tensors,”Journal of Applied Math. Mech., Vol 33, pp 345–356, 1953.

[3] Demirdzic, I., Gosman, A. D., Issa, R. I., and Peric, M., “A Calculation Procedure for Turbulent Flow in Complex Geometries,” Computer & Fluids Vol. 15, No. 3, pp 251–273, 1987.

[4] Takizawa, A., Koshizuka, S. and Kondo, S., “Generalization of Physical Component Boundary Fitted Co-ordinate (PCBFC) Method for the Analysis of Free-Surface Flow,” International Journal for Numerical Methods in Fluids, Vol. 15, pp 1213-1237, 1992.

[5] Warsi, Z. U. A., “Conservation Form of the Navier-Stokes Equations in General Nonsteady Coordinates,” AIAA Journal, Vol. 19, No. 2, 1981.

[6] Swungho Lee, “Three-Dimensional Incompressible Viscous Solutions Based on the Physical Curvilinear Coordinate System,” Dissertation, Department of Computational Engineering, Mississippi State University, December 1997.

[7] Peric, M., “Analysis of Pressure Velocity Coupling on Non-Staggered Grids,” Numerical Heat Transfer, Part B, Vol. 17, pp 63–82, 1990.

[8] Cho, M. J. and Chung M. K., “New Treatment of Nonorthogonal Terms in the Pressure-Correction Equation,” Numerical Heat Transfer, Part B, Vol. 26, pp 133–145, 1994.

[9] Baldwin, B. S. and Lomax, H., “Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows,” AIAA 78–257, 1978.

[10] “Cooperative Experiments on Wigley Parabolic Models in Japan,” 17th ITTC Resistance Committee Report, 1983.

[11] Chen, C. J., Bravo, R. H., Chen, H. C. and Xu, Z., “Accurate Discretization of Incompressible Three-Dimensional Navier-Stokes Equations,” Numerical Heat Transfer, Part B, Vol. 27, pp 371–392, 1995.

[12] Rosenfeld, M. and Kwak, D., “A Fractional Step Solution Method for the Unsteady Incompressible Navier-Stokes Equations in Generalized Coordinate Systems,” Journal of Computational Physics, Vol. 94, pp 102–137, 1991.

[13] Karki, K. C., Sathyamurthy, P. S. and Patankar, S. V., “Performance of a Multigrid Method with an Improved Discretization Scheme for Three-Dimensional Fluid Flow Calculations,” Numerical Hear Transfer, Part B, Vol. 29, pp 275–288, 1996.

Swungho Lee Hydrodynamics Research Dept. Hyundai Maritime Research Institute 1, Cheonha-Dong, Dong-Ku, Ulsan, Korea E-mail address: slee@hhi.co.kr

Bharat K. Soni Department of Aerospace Engineering Mississippi State University Mississippi State, MS 39210 USA E-mail address: bsoni@erc.msstate.edu