\documentclass[twoside]{article}
\usepackage{amsmath, graphicx}
\pagestyle{myheadings}
\markboth{\hfil Simulation of incompressible flows \hfil EJDE--2003/Conf/10}
{EJDE--2003/Conf/10 \hfil Jalal Abedi \& Shahrouz Aliabadi \hfil}
\begin{document}
\setcounter{page}{1}
\title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent
Fifth Mississippi State Conference on Differential Equations and
Computational Simulations, \newline
Electronic Journal of Differential Equations,
Conference 10, 2003, pp 1--11. \newline
http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu (login: ftp)}
\vspace{\bigskipamount} \\
%
Simulation of incompressible flows with heat and mass transfer using
parallel finite element method%
%
\thanks{ {\em Mathematics Subject Classifications:} 65C20, 81T80.
\hfil\break\indent
{\em Key words:} Finite element method, parallel simulation,
contaminant dispersion, \hfil\break\indent free convection.
\hfil\break\indent
\copyright 2003 Southwest Texas State University. \hfil\break\indent
Published February 28, 2003. } }
\date{}
\author{Jalal Abedi \& Shahrouz Aliabadi}
\maketitle
\begin{abstract}
The stabilized finite element formulations based on the SUPG
(Stream\-line-Upwind/Petrov-Galerkin) and PSPG (Pressure-Stabilization/Petrov-Galerkin) methods
are developed and applied to solve buoyancy-driven incompressible flows with heat and
mass transfer. The SUPG stabilization term allows us to solve flow problems at high
speeds (advection dominant flows) and the PSPG term eliminates instabilities associated
with the use of equal order interpolation functions for both pressure and velocity. The
finite element formulations are implemented in parallel using MPI. In parallel
computations, the finite element mesh is partitioned into contiguous subdomains using
METIS, which are then assigned to individual processors. To ensure a balanced load, the
number of elements assigned to each processor is approximately equal. To solve
nonlinear systems in large-scale applications, we developed a matrix-free GMRES
iterative solver. Here we totally eliminate a need to form any matrices, even at the
element levels. To measure the accuracy of the method, we solve 2D and 3D example of
natural convection flows at moderate to high Rayleigh numbers.
\end{abstract}
\allowdisplaybreaks
\numberwithin{equation}{section}
\section{Introduction}
Many geophysical and astrophysical heat transport problems have to do with convection
in a fluid layer with internal energy sources. Within this fluid layer, the flow field is
developed due to the buoyancy forces generated by temperature gradients. These same
forces are also responsible for transporting pollutants, contaminants and chemicals. In
order to adequately simulate this complex flow, the incompressible Navier-Stokes
equations are coupled with heat and mass transfer equations. This coupling results in an
extremely nonlinear system of equations, especially at high Rayleigh numbers.
Historically, many upwind finite difference schemes have been developed to relax these
nonlinearities \cite{b1,f1,r1}, but in most cases these algorithms lack stability and contain
inconsistencies. Thus, stable and consistent numerical algorithms are needed to solve
these problems with a reasonable degree of accuracy.
In this article, we apply stabilized finite element formulations to solve
buoyancy driven
incompressible flows with heat and mass transfer. The stabilizations are based on the
SUPG (Streamline-Upwind/Petrov-Galerkin) and PSPG (Pressure-Stabilization/Petrov-
Galerkin) methods \cite{a1,a2,a3}. The SUPG stabilization term allows us to solve flow problems
at high speeds (advection dominant flows) and the PSPG term eliminates instabilities
associated with the use of equal order interpolation functions for both pressure and
velocity. Prior to the computation, a finite element mesh is partitioned into contiguous
subdomains, which are then assigned to individual processors. To ensure a balanced load,
the number of elements assigned to each processor is approximately equal.
Discretizations of finite element formulations results in coupled, nonlinear systems of
equations which need to be solved at every time step. For complex 3D systems as
encountered in flows affected by heat transfer, this involves millions of unknowns that
require special attention to solve. Primarily, the computational requirements are
immense, requiring parallel supercomputers with hundreds of processors such as the
CRAY T3E to reduce the computing time. A message passing computing paradigm is
needed to allow cross-processor communication and parallel implementation. This is
accomplished using MPI (Message Passing Interface) libraries. The nonlinear systems are
linearized and solved iteratively using the Newton-Raphson method \cite{a3}. The solutions of
the linear systems are also obtained iteratively using an advanced linear solver algorithm
\cite{s1}. For very large systems of equations, a matrix-free iteration strategy is used to obtain
the solution of the nonlinear system.
The 3D applications in buoyancy driven flows are large-scale. Parallel supercomputers with
hundreds of fast processors, such as the CRAY T3E and IBM SP are used to reduce the
computing time. In the parallel implementation, we use a message-passing computing
paradigm, making the cross-processor communication explicit. This is accomplished by using
the MPI (message passing interface) libraries. Prior to the computation, the finite element mesh
is partitioned into contiguous subdomains, and these subdomains are assigned to individual
processors. To ensure load balancing for each processor, each subdomain contains
approximately the same number of elements. Element-level computations are carried out
independently for each processor.
The governing equations for buoyancy-driven will be introduced in Section 2.
The finite element formulations are presented in Section 3. In section 4,
the iterative solution strategy is outlined and finally, the numerical
results are shown in Section 5.
\section{Governing Equations}
The governing equations are the incompressible Navier-Stokes equations coupled with
the heat and mass transfer equations written over the spatial domain $\Omega$
with boundary of $\Gamma$. In non-dimensional form, these equations can
be written as follows:
\begin{gather}
\nabla \cdot \mathbf{u} =0 \label{e1}\\
\rho \big(\frac{\partial \mathbf{u}}{\partial t} + \mathbf{u}\cdot \nabla \mathbf{u}\big)=
\rho \big(\mathbf{g}-\frac{G_r}{R_e^2} \theta \mathbf{n}_g\big)+\nabla\cdot\sigma,\label{e2}\\
\sigma = -p\mathbf{I} +\frac 2{R_e} \epsilon (\mathbf{u}), \label{e3} \\
\epsilon = \frac12 \big(\nabla\mathbf{u}+\nabla\mathbf{u} ^T\big) \label{e4}\\
\frac{\partial \theta}{\partial t}+\mathbf{u}\cdot\nabla\theta
= \frac 1{R_eP_r}\nabla \cdot \nabla \theta +\dot q, \label{e5}\\
\frac{\partial \rho_c}{\partial t} +\mathbf{u}\cdot\nabla\theta\rho_c
= \frac1{R_eP_rL_e} \nabla\cdot \rho_c+ \dot n\,, \label{e6}
\end{gather}
where, $\rho$, $\mathbf{u}$, $p$, $\theta$, $\rho_c$, $\mathbf{g}$, and $\mathbf{n}_g$ are the total density,
velocity, pressure, temperature,
density of contaminant, gravitational force, and the unit vector in the
direction of the gravity, respectively.
The strain tensor is denoted by $\epsilon$ and $\mathbf{I}$ represents the
identity tensor. The Reynolds, Prandtl, Lewis, and Grashoff numbers are
denoted by $R_e$, $P_r$, $L_e$, and $G_r$, respectively.
The sources of heat and mass generation are denoted by $\dot q$ and $\dot n$,
respectively. Here, the total density is the density of air plus the density
of contaminants.
In non-dimensional form, the total density is defined as:
\begin{equation} \label{e7}
\rho = 1+\rho_c
\end{equation}
Since the density of the air is assumed to be much larger that the density
of the contaminant, we assume that the total density is constant and equal
to the density of the air.
\section{Finite element formulations} %sec.3
In the finite element formulations we first define appropriate sets of
trail solution spaces
$S_u^h$, $S_p^h$, $S_\theta^h$, and $S_c^h$
also the weighting function spaces
$V_u^h$, $V_p^h$, $V_\theta^h$, and $V_c^h$
for the velocity, pressure, temperature and the density of contaminant.
The stabilized finite element formulation of
Equations \eqref{e1}--\eqref{e6} can then be written as follows: for all
$\mathbf{w}^h\in V_u^h$, $q^h\in V_p^h$, $\Psi^h\in V_\theta^h$, and
$\Phi^h\in V_c^h$, find $\mathbf{u}^h\in S_u^h$, $p^h\in S_p^h$,
$\theta^h\in S_\theta^h$, and $\rho_c^h\in S_c^h$, such that:
\begin{equation} \label{e8}
\begin{aligned}
&\int_\Omega \mathbf{w}^h\cdot\rho\big(\frac{\partial u^h}{\partial t} +
\mathbf{u}^h\cdot\nabla \mathbf{u}^h-\mathbf{g} + \frac{G_r}{R_e^2}\theta^h\mathbf{n}_g\big)
\,d\Omega
+\int_\Omega \epsilon(\mathbf{w}^h:\sigma(p^h,\mathbf{u}^h)\,d\Omega\\
&+\int_\Omega q^h\nabla\cdot \mathbf{u}^h\,d\Omega
+\sum_{e=1}^{ne}\int_{\Omega^e} \tau_s\nabla\cdot\mathbf{w}^h \rho^h\nabla
\mathbf{u}^h\,d\Omega
\\
&+\sum_{e=1}^{ne}\int_{\Omega^e} \frac{\tau_m}{\rho} \big[\mathbf{u}^h\cdot\nabla\mathbf{w}^h
-\nabla\cdot\sigma(q^h,\mathbf{w}^h)\big]\\
&\times\Big[\rho\big(\frac{\partial u^h}{\partial t}
\mathbf{u}^h\cdot\nabla \mathbf{u}^h-\mathbf{g} + \frac{G_r}{R_e^2}\theta^h\mathbf{n}_g\big)
-\nabla\cdot\sigma(q^h,\mathbf{w}^h)\Big]\, d\Omega \\
&=\int_{\Gamma_{h_u}} \mathbf{w}^h\cdot h\,\Gamma\,,
\end{aligned}
\end{equation}
%
\begin{equation} \label{e9}
\begin{aligned}
&\int_\Omega \Psi^h\big(\frac{\partial \theta^h}{\partial t} +
\mathbf{u}^h\cdot\nabla \theta^h-\dot q\big)\,d\Omega
+\int_\Omega \frac1{R_eP_r}\nabla\Psi^h\cdot \nabla\theta^h\, d\Omega\\
&+ \sum_{e=1}^{ne}\int_{\Omega^e} \tau_\theta\mathbf{u}^h\cdot\Psi^h
\big( \frac{\partial \theta^h}{\partial t} +
\mathbf{u}^h\cdot\nabla \theta^h-\frac1{R_eP_r} \nabla\cdot\nabla\theta^h-\dot q\big)\,d\Omega\\
&= \int_{\Gamma_{h_\theta}} \Psi^h\cdot h_\theta \,d\Gamma\,,
\end{aligned}
\end{equation}
and
\begin{equation} \label{e10}
\begin{aligned}
&\int_\Omega \Phi^h\big(\frac{\partial \rho_c^h}{\partial t} +
\mathbf{u}^h\cdot\nabla \rho_c^h-\dot n\big)\,d\Omega
+\int_\Omega \frac1{R_eP_rL_e}\nabla\Phi^h\cdot \nabla\rho_c^h\, d\Omega\\
&+ \sum_{e=1}^{ne}\int_{\Omega^e} \tau_c \mathbf{u}^h\cdot\nabla \Phi^h
\big( \frac{\partial \rho_c^h}{\partial t} +
\mathbf{u}^h\cdot\nabla \rho_c^h-\frac1{R_eP_rL_e} \nabla\cdot\nabla\rho_c^h-\dot n\big)\,d\Omega\\
&= \int_{\Gamma_{h_\theta}} \Phi^h\cdot h_c\,d\Gamma\,,
\end{aligned}
\end{equation}
where $h$, $h_\theta$, and $h_c$ are the natural boundary conditions
for the velocity, temperature and the density of contaminant, respectively.
The stabilization parameters are defined as follows:
\begin{gather}
\tau_m=\Big[\big(\frac{2\|\mathbf{u}\|}{h}\big)^2+
\big(\frac4{R_eh^2} \big)^2\Big]^{1/2}, \label{e11}
\\
\tau_s=\frac h2\|\mathbf{u}\| z, \quad\mbox{with }
z=\begin{cases} R_u/3 &\mbox{if }R_u\leq 3\\
1 & \mbox{if } 3< R_u, \end{cases} \label{e12}
\\
\tau_\theta=\Big[\big(\frac{2\|\mathbf{u}\|}{h}\big)^2+
\big(\frac4{R_eP_r h^2} \big)^2\Big]^{1/2}, \label{e13}
\\
\tau_c=\Big[\big(\frac{2\|\mathbf{u}\|}{h}\big)^2+
\big(\frac4{R_eP_rL_eh^2} \big)^2\Big]^{1/2}, \label{e14}
\end{gather}
\paragraph{Remark:}
In Equation \eqref{e8}, the first three integrals together with the right hand side term are the
Galerkin formulation of the Navier-Stokes equations. The fourth term (first element level
integral) is the least-square stabilization of the continuity equation. This term enhances
stabilization at high Reynolds number flows. The fifth term (second element level integral)
includes the SUPG and the PSPG stabilizations. The SUPG stabilization eliminated numerical
oscillation due to advection dominancy and the PSPG stabilization allow us to use equal order
interpolation functions for both velocity and pressure.
\paragraph{Remark:}
In Equation \eqref{e9}, the first two integrals together with the right hand side term are the
Galerkin formulation of the heat transfer equation. In this equation, the element level integral is
the SUPG stabilization term.
\paragraph{Remark:}
In Equation \eqref{e10}, the first two integrals together with the right hand side term are the
Galerkin formulation of the mass transfer equation. In this Equation, the element level integral is
the SUPG stabilization term.
The finite element formulations are discretized using first order polynomials for both the
unknowns and the weighting functions \cite{a3}.
\section{Iterative solution strategy}
The discretization of the finite element formulations results in a series of coupled,
nonlinear systems of equations that need to be solved at every time step. The nonlinear
system of equations in vector form can be written as:
\begin{equation} \label{e15}
\mathbf{F}(\mathbf{\dot s}, \mathbf{s})= \mathbf{L}
\end{equation}
where the vector $\mathbf{F}$ is a function of the nodal unknowns
$\mathbf{s}$ and its time derivative $\mathbf{\dot s}$. Here
$\mathbf{L}$ is the known right hand side vector. After linearization using Newton-Raphson
algorithm, at every nonlinear iteration $\tau$, we need to solve the first order linear
differential equation system in this form:
\begin{equation} \label{e16}
\mathbf{M}\Delta \mathbf{\dot s}^{\tau+1}+ \mathbf{N}\Delta \mathbf{s}^{\tau+1}
= \mathbf{L} - \mathbf{F}^\tau\,,
\end{equation}
where
\[ %17-19
\mathbf{M}=\frac{\partial \mathbf{F}}{\partial \mathbf{\dot s}}\,, \quad
\mathbf{N}=\frac{\partial \mathbf{F}}{\partial \mathbf{s}}\,.
\]
Integrating in time, \eqref{e16} results in
\begin{equation} \label{e24}
(\mathbf{M}^{\tau+1}+\xi\Delta t \mathbf{N}^{\tau+1})\Delta \mathbf{s}^{n+1}
= \mathbf{R}^\tau\,,
\end{equation}
where $n$ is the time step counter, $\delta t$
is the time increment and $\mathbf{R}$ is the residual.
Here $\xi$ is
the scheme selector ( $xi=0$ : fully explicit, $\xi=1$: fully implicit,
$\xi= 0.5$: time accurate) \cite{a4}.
In this work, is set $\xi=0.5$ for second order accuracy in time
integration \cite{a4}. The
linear system of equations in \eqref{e24} is also solved iteratively using the GMRES
update algorithm \cite{a4,s1}. For very large systems of equations, we use a matrix-free
iteration strategy \cite{a4,k1}. This element-vector-based computation totally eliminates the
need to form any matrices, even at the element level.
\section{Numerical examples} %sec. 5
\paragraph{Dispersion of Smoke from a Chimney}
In this test problem, smoke is released into the air from
a chimney as shown in Figure 1. The temperature of the outside air is
$15^\circ$C and the smoke has the
temperature of $56^\circ$C. The density of smoke is 10 precent of the air
density.
The simulations are carried out using an unstructured mesh made of 339,513 elements and 57,266
nodes. The square base of the chimney is 1m ×××× 1m. The simulations are carried out for three wind
velocities, 0.25, 0.125 and 0.0625 m/s. The series of pictures in Figure 1 show the iso-surface of
smoke at 5\% of air density. In this Figure, the left, middle and right pictures correspond to the
wind velocity of 0.25, 0.125 and 0.0625 m/s, respectively. The simulations are carried out on
IBM SP-2 with 16 processors.
At a wind velocity of 0.25 m/s, the wind inertia is relatively greater than bouyancy driven inertia
in the vertical direction. Therefore, the dispersion is in the form of plume. This characteristic
changes as the wind velocity diminishes. At the wind velocity of 0.0625 m/s the flow induced by
the bouyancy force is much stronger than the convection flow due to the wind.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.9\textwidth, clip=]{fig1.eps}
\end{center}
\caption{Dispersion of smoke from a chimney. Iso-surface of smoke at 5\%
of air density. The left, middle and right pictures correspond to
wind velocity of 0.25, 0.125, and 0.0625 m/s}
\end{figure}
\paragraph{2-D simulation of natural convection at various Rayleigh numbers}
To measure the accuracy
of the finite element method, we simulate the natural convection of flow with constant heat
source. Scientifically, there are extensive experimental results for
comparison \cite{k2,k3}.
Here, the computational domain is a square with one unit length. The nondimensional heat source
is 2. The Prandtl number is fixed at 6.5 and the Reynolds number is selected
such that $Re P_r=1$.
The computational domain consisted of 22,500 ($150 \times 150$) quadrilateral elements, and
simulations were carried out for Rayleigh numbers of $10^4$ to $10^8$.
These simulations were
carried out on the CRAY T3E with 64 processors. The mass transfer equation is not
solved in this problem; thus the initial conditions included zero pressure, velocity and
temperature everywhere on the domain. A zero velocity boundary condition was imposed
on all four edges, which were also traction-free. In addition, zero temperature was
imposed on the top edge, and no heat flux on the remaining three edges. The comparison
was of the computed Nusselt number (based on the reference or average temperature
inside the computational domain) for the top edge versus measured Nusselt numbers. The
steady state solution is graphically depicted at Rayleigh numbers of
$10^4$--$10^6$ in Figure 2.
Figure 2 shows the time history of the Nusselt number for Rayleigh numbers
of $10^4$, $10^5$, and $10^6$.
In these three cases, the steady-state solution was obtained.
For Rayleigh numbers more than $10^7$,
strong convection makes the flow unstable. Instabilities grow as we increase the Rayleigh
number. However, these instabilities, we believe are physical, since we could not find any
numerical oscillations in the solution. Figure 3 shows the temperature distribution for various
Rayleigh numbers. As can be seen, the symmetric solution is obtained for Rayleigh numbers up
to $10^7$. At Rayleigh number of $10^8$, the instabilities result in strong
nonsymmetric convection.
In experiment \cite{k2,k3}, the Nusselt number can be expressed in the form of
the power law as:
\begin{equation} \label{e15b}
N_u=0.306R_a^{0.227} \quad\mbox{for } P_r=6.5\mbox{ and } 10^4\leq R_a\leq10^8
\end{equation}
In Figure 4, we compare the computed Nusselt numbers with the expression
in \eqref{e15b}.
The comparison is excellent for Rayleigh numbers $10^4$ and $10^8$ and
very comparable for Rayleigh numbers $10^5$, $10^6$, and $10^7$.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.9\textwidth, clip=]{fig2.eps}
\end{center}
\caption{Time history of the Nusselt number}
\end{figure}
\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.9\textwidth, clip=]{fig3.eps}
\end{center}
\caption{Temperature distribution for $Ra=10^4-10^8$ (from left to right)}
\end{figure}
\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.8\textwidth, clip= ]{fig4.eps}
\end{center}
\caption{Comparison between computed and mesasured Nusselt numbers}
\end{figure}
\paragraph{3-D simulation of natural convection at $Ra = 10^7$}
Since the steady state solution for $Ra \geq 10^7$ was not achieved for 2-D
(i.e., strong convection results in unstable physical
fluctuations), the simulation was carried out for 3-D using same initial and boundary
conditions. The computational domain consisted of 3,375,000
($150 \times 150 \times 150$)
quadrilateral elements, and simulations were carried out for Rayleigh number
of $10^7$.
Some results of this simulation are shown in Figure 5, which depicts iso-surfaces of the
temperature after release of heat at different time.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.9\textwidth, clip=]{fig5.eps}
\end{center}
\caption{Iso-surface illustration for temperature}
\end{figure}
Figure 6 shows the time history of the Nusselt number for Rayleigh numbers of
$10^7$. This Figure
compares the result of simulation for 2-D and 3-D. As can be seen, the steady-state solution was
obtained only for 3-D.
\begin{figure}[tb]
\begin{center}
\includegraphics[width=0.8\textwidth,clip=]{fig6.eps}
\end{center}
\caption{Time history of the Nusselt number for 2-D and 3-D}
\end{figure}
\paragraph{Conclusion}
This project demonstrates a highly accurate SUPG finite element method, developed for
buoyancy-driven incompressible flows with heat and mass transfer. The finite element
method was implemented in parallel using MPI libraries, and the accuracy of the method
was demonstrated for flows at moderate to high Rayleigh numbers. The computed results
showed good correlation with experimental or measured data, and show the applicability
of the method to solving 3D applications. Instabilities due to numerical oscillations were
eliminated for 3D model.
\paragraph{Acknowledgement/disclaimer}
We thank the DOD High Performance Computing Modernization Program, ERDC Major
Shared Resource Center through Programming Environment and Training (PET) for
funding this work. We also thank the Army HPC Research Center for partial support and
providing computer time. Views, opinions, and/or findings contained in this report are
those of the authors and should not be constructed as an official Department of Defense
position, policy, or decision unless so designated by other official documentation.
\begin{thebibliography}{00} \frenchspacing
\bibitem{a1} S. K. Aliabadi and T.E. Tezduyar, ``Space-time Finite Element Computation of
Compressible Flows Involving Moving Boundaries and Interfaces", Computer
Methods in Applied Mechanics and Engineering, 107, 209-223 (1993).
\bibitem{a2} S. K. Aliabadi, S. E. Ray and T. E. Tezduyar, ``SUPG Finite Element
Computation of Viscous Compressible Flows Based on the Conservation and
Entropy Variables Formulations", Computational Mechanics, 11, 300-312 (1993).
\bibitem{a3} S. K. Aliabadi and T. E. Tezduyar, ``Parallel Fluid Dynamics Computations in
Aerospace Applications", International Journal for Numerical Methods in
Fluids, 21, 783-805 (1995).
\bibitem{a4} S. Aliabadi and T. Tezduyar, ``Parallel Fluid Dynamics Computations in
Aerospace Applications", International Journal for the Numerical Methods in
Fluids, 21:783-805 (1995).
\bibitem{b1}
Z. H. Barakat and A. J. Clark, ``Analytical and Experimental Study of Transient
Laminar Natural Convection Flows in Partially Filled Liquid Containers",
Proceedings of Third International heat transfer Conference, Chicago, Illinois,
Vol. 2, 152-162 (1966).
\bibitem{f1} E. J. Formm,``The Time Dependent Flow of an Incompressible Viscous Fluid",
Methods in Computational Physics, Vol. 3, 345-382 (1964).
\bibitem{k1} G. Karypis and V. Kumar, ``Parallel Multilevel k-Way Partitioning Scheme for
Irregular Graphs", SIAM Review, 41 (2), 278-300 (1999).
\bibitem{k2} M. Keyhani and F. A. Kulacki, ``Experiments on Transient Thermal Convection
with Internal Heating Large Time Results", ASME Journal of Heat Transfer,
Vol. 5, No. 5, 261-266 (1988).
\bibitem{k3} F. A. Kulacki and A. A. Emara, ``Steady and Transient Thermal Convection in a
Fluid Layer with Uniform Volumetric Energry Sources", Journal of Fluid
Mechanics, Vol. 55, Part 2, 271-287 (1977).
\bibitem{r1} J. P. Roache, ``Computational Fluid Dynamics", Hermosa Publishes,
Albuquerque, N.M.
\bibitem{s1} Y. Saad and M. Schultz, ``GMRES: Generalized Minimal Residual Algorithm for
Solving Nonsymmetric linear Systems", SIAM Journal of Scientific and Statistical
Computing, 7, 856-896 (1986).
\end{thebibliography}
\noindent\textsc{Jalal Abedi} (e-mail: jabedi@cau.edu Phone 404-880-6938)\\
\textsc{Shahrouz Aliabadi} (e-mail: aliabadi@cau.edu)\\[2pt]
Department of Engineering, Clark Atlanta University\\
223 James P. Brawley Dr. S. W.\\
Atlanta, Georgia 30314, USA
\end{document}