\documentclass[twoside]{article} \usepackage{amssymb} % used for R in Real numbers \usepackage{epsfig} \pagestyle{myheadings} \setcounter{page}{149} \markboth{\hfil Governing Equations of Fluid Mechanics \hfil}% {\hfil Swungho Lee \& Bharat K. Soni \hfil} \begin{document} \title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent {\sc Differential Equations and Computational Simulations III}\newline J. Graef, R. Shivaji, B. Soni, \& J. Zhu (Editors)\newline Electronic Journal of Differential Equations, Conference~01, 1997, pp. 149--157. \newline ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu \newline ftp ejde.math.swt.edu (login: ftp)} \vspace{\bigskipamount} \\ Governing Equations of Fluid Mechanics in Physical Curvilinear Coordinate System \thanks{ {\em 1991 Mathematics Subject Classifications:} 65N30, 35Q30. \hfil\break\indent {\em Key words and phrases:} Physical curvilinear coordinate systems, \hfil\break\indent incompressible flows, transformations. \hfil\break\indent \copyright 1998 Southwest Texas State University and University of North Texas. \hfil\break\indent Published November 12, 1998.} } \date{} \author{Swungho Lee \& Bharat K. Soni} \maketitle \begin{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. \end{abstract} \section*{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.~\cite{pat} note that, \begin{quote}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.\end{quote} 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~\cite{tru}. Demirdzic~et~al.~\cite{dem} derived the physical curvilinear-components form in nonorthogonal coordinates for Reynolds-averaged Navier-Stokes equations and the equations of the $k - \epsilon$ 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.~\cite{tak} 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 two-dimensional 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. This work explores an approach which is different from Demirdzic~et~al~\cite{dem} 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. \section*{Basic transformations} In the following analysis, $x_i$ are Cartesian coordinates, $\xi^i$ are curvilinear coordinates and $\xi^{(i)}$ are physical curvilinear coordinates. First, consider the general transformation law under the two changes of the coordinates system $x_i$ to $\xi^i$ and $\xi^i$ to $\xi^{(i)}$. The relationship between the coordinates $x_i$ and the coordinates $\xi^i$ can be expressed as follows: $$\xi^i=\xi^i(x_j, t)$$ The relationship between the coordinates $\xi^i$ and the coordinates $\xi^{(i)}$ can be expressed as: $$\xi^{(i)} = \xi^{(i)}(\xi^i)\,, \textrm{ where }\Delta\xi^{(i)} = \sqrt{\mathrm{g}_{ii}}\Delta\xi^i$$ $\sqrt{\mathrm{g}_{ii}}$ are evaluated at $\xi^k$ = constant and $k \neq i$. $\xi^{(i)}$ resemble the coordinate stretching in each direction of $\xi^i$. In view of transforming the coordinates from $x_i$ to $\xi^i$ and $\xi^i$ to $\xi^{(i)}$, the vector $d\vec r$ can be written as: \begin{eqnarray} d\vec r & = & {\frac{\partial\vec r}{\partial\xi^i}}d\xi^i = \vec a_id\xi^i\qquad\textrm{(sum on i)} \\ & = & {\frac{\partial\vec r}{\partial\xi^{(i)}}}d\xi^{(i)} = \vec a_{(i)}d\xi^{(i)}\,, \textrm{ where } \vec a_{(i)} = {\frac{1}{\sqrt{\mathrm{g}_{ii}}}}\vec a_i \end{eqnarray} $\vec a_i$ are covariant base vectors in the curvilinear coordinate system and $\vec 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: \begin{eqnarray} \vec a_{(i)} \cdot \vec a_{(j)} = \mathrm{g}_{(ij)} & = & {\frac{1}{{\sqrt{\mathrm{g}_{ii}}}{\sqrt{\mathrm{g}_{jj}}}}} \mathrm{g}_{ij} \,, \textrm{ where } \vec a_{i}\cdot \vec a_{j} = \mathrm{g}_{ij} \\ \vec a^{(i)}\cdot \vec a^{(j)} = \mathrm{g}^{(ij)} & = & \sqrt{\mathrm{g}_{ii}}\sqrt{\mathrm{g}_{jj}}\mathrm{g}^{ij} \,, \textrm{ where } \vec a^{i}\cdot \vec a^{j} = \mathrm{g}^{ij}\,, \quad \vec a^{(i)} = \sqrt{\mathrm{g}_{ii}}\vec a^{i} \end{eqnarray} $\mathrm{g}_{(ij)}$ and $\mathrm{g}^{(ij)}$ are the physical covariant and the physical contravariant metric tensors, respectively. The physical covariant metric tensor $\mathrm{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^\mathrm{th}$ component projected onto the $i^\mathrm{th}$ physical curvilinear coordinate direction, $$u^{(i)} = \vec u \cdot \vec a_{(i)}\,.$$ 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. \begin{eqnarray} \label{six} {\frac{\partial\vec a_{(i)}}{\partial\xi^{(j)}}} & = & [{\frac{\sqrt{\mathrm{g}_{kk}}}{{\sqrt{\mathrm{g}_{ii}}} {\sqrt{\mathrm{g}_{jj}}}}}\Gamma_{ij}^{k}-\delta_{i}^{k} {\frac{\mathrm{g}_{km}}{{\mathrm{g}_{ii}}{\sqrt{\mathrm{g}_{jj}}}}} \Gamma_{ij}^{m}\,]\vec a_{(k)} \quad \textrm{(sum on k and m)} \nonumber \\ & = & \Gamma_{(ij)}^{(k)}\vec a_{(k)} \,\, \textrm{(sum on k)} \,, \end{eqnarray} where $$\Gamma_{(ij)}^{(k)} = {\frac{\sqrt{\mathrm{g}_{kk}}}{{\sqrt{\mathrm{g}_{ii\!}}} {\sqrt{\mathrm{g}_{jj\!}}}}}\Gamma_{ij}^{k}-\delta_{i}^{k} {\frac{\mathrm{g}_{km}}{{\mathrm{g}_{ii\!}}{\sqrt{\mathrm{g}_{jj\!}}}}} \Gamma_{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{\partial \vec a^{(i)}}{\partial \xi^{(j)}}} = -\Gamma_{(kj)}^{(i)}\vec a^{(k)} \quad \textrm{(sum on k)}$$ Here the Christoffel and physical counterparts of the Christoffel symbols have the following properties: $$\Gamma_{jk}^{i} = \Gamma_{kj}^{i} \,\, \textrm{and} \,\, \Gamma_{(jk)}^{(i)} \neq\Gamma_{(kj)}^{(i)}$$ Using the general transformation laws for a scalar $\phi$, the gradient can be written as follows: $$\nabla\phi = {\frac{\partial\phi}{\partial\xi^{(i)}}}\vec a^{(i)} = \mathrm{g}^{(ik)} {\frac{\partial\phi}{\partial\xi^{(k)}}}\vec a_{(i)} \quad \textrm{(sum on i and k)}\,,$$ where $\vec a^{(i)} = \mathrm{g}^{(ik)}\vec a_{(k)}$. Also, the gradient of a vector $\vec u$ can be expressed using equation as: \begin{eqnarray*} \nabla\vec u & = & {\frac{\partial\vec u}{\partial\xi^{(i)}}}\vec a^{(i)} = {\frac{\partial(u^{(k)}\vec a_{(k)})}{\partial\xi^{(i)}}}\vec a^{(i)} \quad \textrm{(sum on i and k)} \\ & = & u_{(,i)}^{(k)}\vec a_{(k)}\vec a^{(i)} = \mathrm{g}^{(ij)}u_{(,i)}^{(k)}\vec a_{(k)}\vec a_{(j)} \,, \end{eqnarray*} where $$u_{(,i)}^{(k)} = {\frac{\partial u^{(k)}}{\partial\xi^{(i)}}} + u^{(m)}\Gamma_{(im)}^{(k)}\,.$$ The quantity $u_{(,j)}^{(k)}$ is called the covariant derivative of the physical curvilinear components of a vector $\vec u$. One can easily evaluate the divergence of a vector $\vec u$ using equation (\ref{six}) as: \begin{eqnarray} \nabla\cdot\vec u & = & {\frac{\partial\vec u}{\partial\xi^{(i)}}}\cdot\vec a^{(i)} = {\frac{\partial u^{(k)}\vec a_{(k)}}{\partial\xi^{(i)}}}\cdot\vec a^{(i)} \quad \textrm{(sum on i and k)} \nonumber \\ & = & u_{(,i)}^{(i)} = {\frac{\sqrt{\mathrm{g}_{ii}}}{{J}}} {\frac{\partial}{\partial\xi^{(i)}}} ({\frac{{J}}{\sqrt{\mathrm{g}_{ii}}}}u^{(i)}) \quad \textrm{(sum on i)} \end{eqnarray} The Laplacian can be evaluated by the divergence of the gradient of a vector $\vec u$, \begin{eqnarray*} \nabla^{2} \vec u &\! = &\! \nabla \cdot \hat{T} = {\frac{\partial \hat{T}}{\partial\xi^{(j)}}} \cdot \vec{a}^{(j)} \,, \textrm{ where } \quad \hat{T} = \nabla_{\vec u} = u_{(i,)}^{(k)} \vec a_{(k)} \vec a^{(i)} \\ &\! = &\! \textrm{g}^{(jk)}{} [u_{(,j)}^{(m)} \Gamma_{(m\!k)}^{(i)} -{} u_{(,m)}^{(i)} \Gamma_{(j\!k)}^{(m)} + {\frac{\partial u_{(,j)}^{(i)}}{\partial \xi^{(k)}}}] \vec a_{(i)} \, \textrm{(sum on i, j, k and m)} \end{eqnarray*} The time derivative of a scalar or a vector F is given as \mbox{Warsi}~\cite{war}: \begin{eqnarray} \label{13a} {\frac{\partial F}{\partial t}} |_{x_i} & = & [{\frac{\partial F}{\partial \tau}} + {\frac{\partial F}{\partial \xi^{(i)}}} {\frac{\partial \xi^{(i)}}{\partial t}}]|_{\xi^(i)} \\ \label{13b} {\frac{\partial F}{\partial \tau}} |_{\xi^(i)} & = & [{\frac{\partial F}{\partial t}} + {\frac{\partial F}{\partial x_{i}}} {\frac{\partial x_{i}}{\partial \tau}}] |_{x_i} \end{eqnarray} Here $\tau$ represents the time in an unsteady physical curvilinear coordinate system. The grid speed can be easily evaluated by replacing $F$ with $\xi^{(i)}$ in equation (\ref{13b}). $$w^{(i)} = {\frac{\partial \xi^{(i)}}{\partial t}} = - {\frac{\partial \xi^{(i)}}{\partial x_{l}}} {\frac{\partial x_{l}}{\partial \tau}} = - {\frac{{\sqrt{\mathrm {g}}_{ii}}}{J}} b_{l}^{i} {\frac{\partial x_{l}}{\partial \tau}} \quad \textrm{(sum on l)}$$ where $$b_{l}^{i} = {\frac{\partial x_{m}}{\partial \xi^{j}}} {\frac{\partial x_{n}} {\partial \xi^{k}}}- {\frac{\partial x_{n}}{\partial \xi_{j}}} {\frac{\partial x_{m}}{\partial \xi^{k}}}$$ with $i, j, k$, and $l, m, n$ in each cyclic order. The divergence of the grid speed vector $\vec w$ is written as: $$\nabla \cdot \vec w = {\frac{{\sqrt{\mathrm {g}}_{ii}}}{J}} {\frac{\partial({\frac{J} {\sqrt{\mathrm {g}}_{ii}}} w^{(i)})}{\partial \xi^{(i)}}} = - {\frac{l}{J}} {\frac{\partial J}{\partial \tau}}$$ \section*{Governing equations in the unsteady \\ physical curvilinear coordinate system} By replacing $F$ with $\vec u$ in equation (\ref{13a}), 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 (\ref{16}). $$\label{16} {\frac{\partial \vec u}{\partial \tau}} + (\nabla \vec u) \cdot \vec v = - \nabla P + \nu_{E} \nabla^{2} \vec u + [\nabla \vec u + (\nabla \vec u)^{T}] \cdot \nabla \nu_{E}\,,$$ where $$\vec v = \vec u + \vec w,{} P = \rho + {\frac{z}{Fn^{2}}} + {\frac{2}{3}}k \, \textrm{and }\, \nu_{E} = {\frac{l}{R_{\mathrm{eff}}}} = {\frac{l}{R_{\mathrm{e}}}} + \nu_{t}\,.$$ Here $\vec w$ is a grid speed vector, $Fn$ is a Froude number, $P$ is total pressure, and $\rho$ is a static pressure. $\nu_{t}$ and $R_{\mathrm{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: \begin{eqnarray} &&{\frac{\partial}{\partial \xi^{i}}} ({\frac{J}{\sqrt\mathrm{g}_{ii}}} u^{(i)}) \ = \ 0 \nonumber \\ \lefteqn{ R_{\mathrm{eff}} {\frac{\partial u^{(i)}}{\partial\tau}} + R_{\mathrm{eff}} v^{(j)} ({\frac{\partial u^{{(i)}}}{\partial \xi^{(j)}}} + u^{(k)} \Gamma_{(kj)}^{(i)} ) } \label{17b} \\ &=& - R_{\mathrm{eff}}\, [\mathrm{g}^{(ij)} {\frac{\partial P}{\partial \xi^{(j)}}} - {\frac{\partial \nu_{E}}{\partial \xi^{(j)}}} [\mathrm{g}^{(jk)} {\frac{\partial u^{(i)}}{\partial\xi^{k}}} + \mathrm{g}^{(jk)} \Gamma_{(lk)}^{(i)} u^{(l)} + \mathrm{g}^{(ik)} {\frac{\partial u^{(i)}}{\partial\xi^{k}}} \nonumber \\ &&+ \mathrm{g}^{(ik)} \Gamma_{(lk)}^{(j)} u^{(l)}]\,] + \mathrm{g}^{(jk)} [{\frac{\partial^{2} u^{(i)}}{\partial \xi^{(k)} \partial \xi^{(j)}}} + {\frac{\partial (u^{(l)} \Gamma_{(lj)}^{(i)})}{\partial \xi^{(k)}}} + \Gamma_{(lk)}^{(i)} u_{,(j)}^{(l)} - \Gamma_{(jk)}^{(l)} u_{,(l)}^{(i)}] \nonumber \end{eqnarray} Equation (\ref{17b}), which has a nonconservative form, is rearranged into the standard form for the use of the finite analytic method \cite{swun}. Using the stretched coordinates $\xi^{i}\ast$, the 12-point finite analytic discretization scheme based on the local nonuniform grid spacing \cite{swun} is employed. The stretched coordinates $\xi^{i}\ast$ are defined as: \begin{displaymath} \xi^{i} = {\sqrt{\mathrm{g}^{ii}}} \xi^{i}\ast \textrm{ or } \xi^{i}\ast = {\frac{l}{\sqrt{\mathrm{g}^{ii}}}} \xi^{i} \end{displaymath} The first equation shows the relation between the curvilinear coordinates and the stretched coordinates. The $\sqrt{\mathrm{g}^{ii}}$ are evaluated at $\xi^{k}$ = constant and $k\neq i$. \section*{Results and discussions} An unsteady three-dimensional incompressible flow solver based on the physical curvilinear coordinate system has been developed \cite{swun}. 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 \cite{swun}. The results of the following two test cases to validate the pertinent numerical simulations are presented here. \subsection*{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 \mbox{$16\!\times\!16\!\times\!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 \cite{bal}. 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. \begin{figure}%[hbtp] \epsfig{file=fig1.eps, width=4.5in} \caption{Convergence Histories for Different Time Steps} \end{figure} \begin{figure}%[hbtp] \epsfig{file=fig2.eps, width=4.5in} \caption{Comparison of the Centerline u-Velocity Profile} \end{figure} 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 $\Delta t = 1$. \mbox{Peric~\cite{per}} has reported that if the angle between the two coordinate lines is greater than $135^\circ$ or less than $45^\circ$, then the pressure correction equation does not converge at all, or the convergence rate is too slow. \mbox{Cho and Chung~\cite{cho}} 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 the angle, the narrower the region of relaxation factor. In the present research, the computations are performed at several inclined angles \\ ($90^\circ, 60^\circ, 45^\circ, 30^\circ, 15^\circ,10^\circ, \textrm{and} 5^\circ$) to check the influence on the rate of the convergence due to the grid skewness. As mentioned, all computations of the three-dimensional 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. \begin{figure}%[hbt] \epsfig{file=fig3.eps,width=4.5in} \caption{Convergence Histories in Various Inclined Angles (Degree) } \end{figure} \subsection*{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. \begin{table} [!hbp] \caption{Grid Dependence Tests for the Wigley Hull} \label{table 1} \begin{center} \begin{tabular}{|lccc|} \hline & I & II & III \\ \hline Grid Points & $125\times35\times34$ & $125\times40\times40$ & $125\times50\times48$ \\ 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\times10^6$ & $1.0\times10^6$ & $1.0\times10^6$ \\ Froude Number & 0.289 & 0.289 & 0.289 \\ \hline \end{tabular} \end{center} \end{table} Table~\ref{table 1} shows the information for these computations. The computations were performed with three different grids under the same conditions. The Baldwin-Lomax turbulence model \cite{bal} was used to calculate the eddy viscosity in the turbulent flow. The Froude number and the Reynolds number used in the experiment \cite{coop} are 0.289 and $3.3 \times 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 ($\Delta \textrm{x}$) in a region of relatively high gradient. The residue remains around $10^{-4}$after $t = 0.5$. \begin{figure}%[hbtp] \epsfig{file=fig4.eps,width=4.5in} \caption{The Wave Elevations on the Hull Surface for different grid sizes and a Perspective view of the wave elevation} \end{figure} 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. \begin{thebibliography}{99} \bibitem{pat} 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. \bibitem{tru} Truesdell, C., The Physical components of vectors and Tensors,''Journal of Applied Math. Mech., Vol~33, pp~345--356, 1953. \bibitem{dem} 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. \bibitem{tak} 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. \bibitem{war} Warsi, Z.~U.~A., Conservation Form of the Navier-Stokes Equations in General Nonsteady Coordinates,'' AIAA Journal, Vol.~19, No.~2, 1981. \bibitem{swun} Swungho Lee, Three-Dimensional Incompressible Viscous Solutions Based on the Physical Curvilinear Coordinate System,'' Dissertation, Department of Computational Engineering, Mississippi State University, December 1997. \bibitem{per} Peric, M., Analysis of Pressure Velocity Coupling on Non-Staggered Grids,'' Numerical Heat Transfer, Part~B, Vol.~17, pp~63--82, 1990. \bibitem{cho} 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. \bibitem{bal} Baldwin, B.~S. and Lomax, H., Thin Layer Approximation and Algebraic Model for Separated Turbulent Flows,'' AIAA 78--257, 1978. \bibitem{coop} Cooperative Experiments on Wigley Parabolic Models in Japan,'' 17th ITTC Resistance Committee Report, 1983. \bibitem{chen} 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. \bibitem{ros} 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. \bibitem{kar} 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. \end{thebibliography} \bigskip \noindent{\sc Swungho Lee} \\ Hydrodynamics Research Dept. \\ Hyundai Maritime Research Institute \\ 1, Cheonha-Dong, Dong-Ku, Ulsan, Korea \\ E-mail address: slee@hhi.co.kr \medskip \noindent{\sc Bharat K. Soni} \\ Department of Aerospace Engineering \\ Mississippi State University \\ Mississippi State, MS 39210 USA \\ E-mail address: bsoni@erc.msstate.edu \end{document}