\documentclass[reqno]{amsart}
\usepackage{hyperref}
\usepackage{graphicx}

\AtBeginDocument{{\noindent\small 
Seventh Mississippi State  - UAB Conference on Differential Equations 
and Computational Simulations, 
{\em Electronic Journal of Differential Equations},
Conf. 17 (2009),  pp. 149--157.\newline 
ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or 
http://ejde.math.unt.edu
\newline ftp ejde.math.txstate.edu}
\thanks{\copyright 2009 Texas State University - San Marcos.}
\vspace{9mm}}

\begin{document} \setcounter{page}{149}
\title[\hfilneg EJDE-2009/Conf/17\hfil Spectral element simulations]
{Spectral element simulations of flow past an ellipsoid at
different Reynolds numbers}

\author[D. Liu, G. Karniadakis, M. Maxey\hfil EJDE/Conf/17 \hfilneg]
{Don Liu, George Karniadakis, Martin Maxey}  

\address{Don Liu \newline 
Mathematics and Statistics, Louisiana Tech University, Ruston, LA
71270, USA} 
\email{donliu@LaTech.edu Tel: 318 257 4670}

\address{George Karniadakis \newline
Division of Applied Mathematics, Brown University, Providence, RI
02912, USA} 
\email{George\_Karniadakis@brown.edu}

\address{Martin Maxey \newline
Division of Applied Mathematics, Brown University, Providence, RI
02912, USA} 
\email{Maxey@dam.brown.edu}


\thanks{Published April 15, 2009.}
\subjclass[2000]{74S25, 74S05, 76T99}
\keywords{Two phase flow; spectral element method; ellipsoid;
\hfill\break\indent
direct numerical simulation}

\begin{abstract}
 A new method is developed to simulate the fully coupled motion
 involving an  ellipsoidal particle and the ambient fluid.
 This method essentially reduces a two-phase flow problem into a
 single-phase fluid flow problem. The momentum exchange at the
 interface between the particle and the fluid is computed in a
 spectral element method. To validate this method and demonstrate
 the accuracy of the results, theoretical data as well as results
 from the direct numerical simulations (DNS) via a spectral element
 method are obtained to provide reference and comparison.
\end{abstract}

\maketitle
\numberwithin{equation}{section}

\section{Introduction}
\label{chap:Introduction}

Many industrial, mechanical and biomedical phenomena involve  two
phase flows. The mutual interactions between particles and a
carrier fluid have been studied  experimentally and numerically,
for example \cite{Hu96,MaxeyPCW97,LiuMK02,LiuMK04}. However, most
of the relevant research papers deal with spherical particles.
Two-phase flows with non-spherical particles are quite commonly
met but not as well studied as spherical particles because of the
extra complexity in the ellipsoidal interface. Therefore two-phase
flow involving non-spherical particles becomes interesting to
researchers. This paper presents a spectral/hp element simulation
of an ellipsoidal particle in confined domains. A novel method to
describe the two-way coupled motion between the fluid and the
ellipsoid is addressed here. The novelty of this method lies in
the fact that it fully utilizes the high order accuracy of a
spectral element method and also essentially reduces a two-phase
flow problem into a single-phase flow problem.


In this paper, a mathematical description of this method is
presented,  especially about how the two-way coupled motion is
handled, and how an ellipsoid is represented in the computational
domain. This method is related to the low order multipole
expansion method, based on the Force-Coupling Method (FCM) which
was developed by Maxey et al. \cite{MaxeyP01}. The presence of the
ellipsoidal particle is implicitly included in the mathematical
formulation by a set of terms called force monopole and  torque
dipole terms. The forces and torques involved in the mutually
coupled motion are distributed finitely via tailored Gaussian
distribution functions with the mean being coordinates of the
center of the ellipsoid. The envelopes of the Gaussian functions
are selected according to the orientation and the size of the
semi-axes of the ellipsoid. A spectral/hp element method is used
to obtain the numerical solutions of this two phase flow problem
with the immersed ellipsoid.


Among a series of systematic simulation results, two benchmark
simulation examples, a Stokes flow case and a finite Reynolds
number flow case are selected and presented in this paper. The
induced drag, lift and torque are verified and compared with the
analytic data from Happel et al. \cite{HappelB86} as well as
results from the direct numerical simulation (DNS) on a
body-fitted mesh. The associated computational errors are
tabulated in detail.

\section{Mathematical Formulation and Simulation Methods}\label{chap:Formulation}

The two-way coupling between the particle and fluid is the core
of two-phase flow problems. To describe the particle and fluid
motion efficiently, artificial source terms are introduced in
Navier-Stokes equations to account for the mutual effects between
the ellipsoid and the ambient fluid. The governing equations of
the fluid momentum in terms of the primitive variables ${\bf u}$
and p are given as:
\begin{gather}
\rho\frac{D{\bf u}}{D t}=-\nabla
   p+\mu\nabla^{2}{\bf u}+{\bf f}({\bf x}, t), \label{NS}\\
\nabla\cdot{\bf u}=0. \label{Continue}
\end{gather}
Where $\rho$, $p$ and $\mu$ are the fluid density, pressure and
viscosity, respectively. The source term ${\bf f}({\bf x},t)$
describes the two-way mutual interactions between the fluid and
the ellipsoid centered at ${\bf Y}(t)$, and is the sum of the
force monopole (the first term) and the dipole (the second term)
in (\ref{source}):
\begin{equation}
{\bf f}({\bf x},t)={\bf F}\Delta(\vec{\sigma}_m,{\bf x}-{\bf
Y}(t)) + {\bf G} \nabla \big[\Delta(\vec{\sigma}_d,{\bf x}-{\bf
Y}(t))\big], \label{source}
\end{equation}
where $\Delta(\vec{\sigma}_m,{\bf x})$ and $\Delta(\vec{\sigma}_d,{\bf x})$
are the monopole and dipole Gaussian distribution functions with the monopole
envelope vector $\vec{\sigma}_m=(\sigma_{m1},\sigma_{m2},\sigma_{m3})$ and
the dipole envelope $\vec{\sigma}_d=(\sigma_{d1},\sigma_{d2},\sigma_{d3})$
respectively. For spherical particles \cite{LiuMK04}, one Gaussian envelope
is used in all directions since spheres are isotropic in space. For an
ellipsoid, however, the semi-axes are different in three directions.
Therefore different envelopes are needed in different directions for both the
monopole and dipole Gaussian functions. For example, the monopole Gaussian
distribution function for an ellipsoid with the principal axes aligned with
the fixed Cartesian coordinate axes $(x_1, x_2, x_3)$ can be factorized as:
\begin{equation}
\Delta(\vec{\sigma}_m,{\bf x}-{\bf Y})=\Delta_1(\sigma_{m1},x_1-Y_1)
\Delta_2(\sigma_{m2},x_2-Y_2)\Delta_{3}(\sigma_{m3},x_3-Y_3),
\label{ELL1}
\end{equation}
where $\Delta_1$, $\Delta_2$ and $\Delta_3$ can be computed in the
same  way as below with the corresponding envelope: $\sigma_{m1}$,
$\sigma_{m2}$ or $\sigma_{m3}$, respectively:
\begin{equation}
\Delta_1(\sigma_{m1},x_1-Y_1)=
\frac{1}{\sqrt{2\pi}\sigma_{m1}}e^{-\frac{(x_1-Y_1)^2}{2\sigma^2_{m1}}}.
\label{ELL2}
\end{equation}
These envelopes are related to the semi-axes of the ellipsoid in the following way:
\begin{gather}
\sigma_{m1}=\frac{a_1}{\sqrt{\pi}},\quad
\sigma_{m2}=\frac{a_2}{\sqrt{\pi}},\quad
\sigma_{m3}=\frac{a_3}{\sqrt{\pi}}, \label{envelop1}
\\
\sigma_{d1}=\frac{a_1}{\sqrt[3]{6\sqrt{\pi}}},\quad
\sigma_{d2}=\frac{a_2}{\sqrt[3]{6\sqrt{\pi}}},\quad
\sigma_{d3}=\frac{a_3}{\sqrt[3]{6\sqrt{\pi}}},\quad .
\label{envelop2}
\end{gather}

The monopole strength $\bf F$ is a force vector and equals the sum of the
forces on this particle. The components of the dipole ${\bf G}$ in
(\ref{source}) has the dimension of a torque. It is a 3 by 3 matrix
for each particle and is  determined via the Dipole Iteration
Scheme \cite{LiuPhD04}.

The particle phase is described in Lagrangian approach.
The particle velocity is filtered out as
\begin{equation}
{\bf V}(t)=\iiint_{\rm Vol}
         {\bf u}({\bf x},t)\Delta(\vec{\sigma}_m,{\bf x}-{\bf Y}(t))
\,d{\rm Vol}.
\label{PartVelo}
\end{equation}
The particle angular velocity $\bf \Omega$ is calculated similarly from the fluid
vorticity and dipole Gaussian distribution function. The particle position $Y(t)$
is then updated by integrating the following equation
\begin{equation}
{\bf V}(t)=\frac{d{\bf Y}(t)}{dt}.
\label{PartPosit}
\end{equation}
Therefore the particle moves to the new position at the next time
step. The details of the Force-Coupling method are given by Maxey
et al. \cite{MaxeyP01} and Lomholt \cite{Lomholt01}.


Most problems in reality are {\it mobility} problems,  where the
forces on the moving particles are known and the motion is to be
found. Another category of problems is called {\it resistance}
problems, where the motion is known but the force and torque to
maintain the known motion is to be determined. To simulate this
resistance problem, a novel method, which is related to the
penalty methods, is developed here to calculate the unknown force
${\bf F}(t)$ and torque $T_{ij}(t)$ on the ellipsoid:
\begin{equation}
\frac{d{\bf F}(t)}{dt}=\lambda[{\bf V}^T(t)-{\bf V}(t)],
\label{PenaltyF}
\end{equation}
and
\begin{equation}
\frac{dT_{ij}(t)}{dt}=\lambda\epsilon_{ijk}[\Omega^T_k(t)-\Omega_k(t)].
\label{PenaltyT}
\end{equation}
Where the target velocity ${\bf V}^T$ and the target angular
velocity  $\Omega^T_k$ are zero because the ellipsoid is fixed
still in this resistance problem. In Eqs. (\ref{PenaltyF}) and
(\ref{PenaltyT}), $\lambda$ is the penalty parameter;
$\epsilon_{ijk}$ is the permutation symbol. A spectral element
method using both the {\it h} and {\it p} type refinements is used
to solve the above problem in terms of the primitive variables.
Details of the spectral element methods can be found in
Karniadakis \cite{KarniadakisS1999}.

\section{Simulation Results}\label{chap:Results}

Flow past an ellipsoid is studied here with the proposed method.
Two benchmark example problems - a Stokes flow and a low Reynolds
number flow are presented in this paper. Simulation results are
compared with the analytic data from Happel et al.
\cite{HappelB86} and also with the results from direct numerical
simulation. As a reference for comparison, the analytic data is
reliable and accurate. The ellipsoid studied here has the
semi-axes $a_1=2$, $a_2=a_3=1$, which are aligned with the fixed
coordinate axes. The ellipsoid is fixed in an off-centered
location inside a channel, which is the computational domain. The
dimensions of the channel are $-80\leq x_1/a_2 \leq 80$ in the
streamwise direction, $-10/3\leq  x_2/a_2\leq 10$ in the lateral
direction, and $-20\leq x_3/a_2\leq 20$ in the spanwise direction.
The origin of the coordinate system is set at the center of the
ellipsoid. A pressure gradient is imposed so that the ellipsoid is
immersed in a fully developed Poiseuille flow. The density and
dynamic viscosity of the fluid are non-dimensionalized into 1.
Periodic boundaries are specified in the streamwise $x_1$ and
spanwise $x_3$ directions.  Non-slip boundary conditions are
specified in the lateral $x_2$ direction.

The computational domain for the simulation with the FCM is the
entire  channel including the volume which is supposed to be
occupied by the ellipsoid since a virtual particle is imposed in
this channel. To compare the flow details, a direct numerical
simulation with a spectral element method is conducted under the
same boundary conditions. The computational domain is inside the
channel and outside the ellipsoid. Non-uniform structured meshes
were created to discretize both the computational domains.

\section{Stokes Flow} \label{chap:StokesFlow}

We consider a resistance problem in a Stokes flow. A dimensionless
pressure gradient of 0.06075 is imposed to drive the flow over a
fixed ellipsoid. The dimensionless viscosity is $\mu=1$. Upon
convergence, the centerline fluid velocity at the entrance is
$u_{cl}=1.366$, and the approach velocity (Happel et al.
\cite{HappelB86}, p.333) is $u_{o}=1.0245$. The approach velocity
is $3/4$ of the centerline velocity due to its location. Although
the ellipsoid is off-centered, no lift force is induced on the
ellipsoid because of the nature of a Stokes flow. However the drag
and torque acting on this ellipsoid are present and are to be
determined. A penalty method, which is new in this paper, is
developed in order to determine the {\it restoring force} and the
{\it restoring torque}, which balance the drag and torque from the
fluid, see Eqs. (\ref{PenaltyF}) and (\ref{PenaltyT}). From the
simulation result, the scaled restoring force in $x_1$ is
$F_1/(\mu a_2 u_o)=29.265$ and the scaled restoring torque in
$x_3$ is $T_{12}/(\mu a^2_2 u_o)=4.376$. The particle velocity and
angular velocity  approach to zero with small errors as shown in
table \ref{tab:AlignedSpheroidStokes-1}.


From the analytic data - table 7-5.1 in Happel et al.
\cite{HappelB86},  the drag normalized by $\mu a_2 u_o$ is 28.816.
Compared with this analytic value, the restoring force has an
error of $1.56\%$ in table \ref{tab:AlignedSpheroidStokes-1}. The
difference in the scaled drag between the DNS result and this
analytic value is only $0.76\%$. The value of the analytic torque
is unavailable for comparison; therefore the error in the torque
given in table \ref{tab:AlignedSpheroidStokes-1} is relative to
the DNS result. To see the periodic effect in both the streamwise
and the spanwise boundary conditions, different channel sizes are
tested in the subsequent parametrical studies in in Liu
\cite{LiuPhD04}. For example, a smaller geometry $-10\leq
x_3/a_2\leq 10$ increases the drag to $5.48\%$ of the analytic
value. When the size of the channel is further reduced in both
$x_1$ and $x_3$ directions to $-20\leq x_1/a_2\leq 20$ and $-5\leq
x_3/a_2 \leq 5$, an increased drag is anticipated. The simulation
result in \cite{LiuPhD04} confirms that the drag actually
increases up to $11\%$.

\begin{table}[ht]
\begin{center}
\begin{tabular}{||c||c|c|c|c||}\hline\hline
Source&Drag/$(\mu a_2 u_o)$
&Torque/$(\mu a^2_2 u_{o})$& $V_{1}/u_o$& $\Omega_3 a/u_o$\\\hline\hline
\cite{HappelB86} & $28.816$  &    N/A     & 0        & 0     \\\hline
Simul.& 29.265(1.56\%)&4.376(3.5\%)&0.000973(0.1\%)&0.000763(0.1\%)\\\hline
DNS        & 28.598(0.76\%) & 4.231 & 0    & 0\\\hline
\end{tabular}
\end{center}
\caption{Comparisons of the scaled results from the analytic data, the DNS, and the
simulation with the FCM for an aligned ellipsoidal particle in a Stokes flow. The numbers
in percentage show the relative errors against the analytic data or DNS.}
\label{tab:AlignedSpheroidStokes-1}
\end{table}


Further validation with the analytic results \cite{HappelB86} was
conducted by setting the semi-axes of the ellipsoid to be
$a_1=a_2=a_3=1$. The measures of the channel are $-80\leq
x_1/a_2\leq 80$, $-10\leq x_3/a_2 \leq 10$, and the dimension in
$x_2$ remains the same. The pressure gradient is set to be 0.05
and the dynamic viscosity is still 1. The centerline fluid
velocity at the entrance becomes 1.152 and the approach velocity
is 0.864. From table 7-5.1, Happel et al. \cite{HappelB86},
$a/c=1$ and $c/l=0.3$, the error in drag is $4.2\%$. If the
measure in $x_3$ is expanded to $-20\leq x_3/a_2\leq 20$ and the
pressure gradient is set to 0.06075, the centerline fluid velocity
at the entrance becomes 1.395 and the approach velocity is
1.04625. From table 7-5.1 in \cite{HappelB86}, the error in the
drag drops down to $1.1\%$ of the analytic value. Further details
are documented in Liu \cite{LiuPhD04}.


To verify the details of the flow field, a set of DNS with a
spectral  element method on different meshes were conducted. The
computational domain is the volume inside the channel of the same
size as before and outside the ellipsoid of semi-axes $a_2=2$,
$a_2=a_3=1$. Similar boundary conditions and pressure gradient are
imposed in the DNS. From the Gauss-Legendre-Lobatto integrations,
the non-dimensional values for the drag and the torque on the
ellipsoid are found to be 28.598 and 4.23, respectively. Compared
with the analytic value from Happel et al. \cite{HappelB86}, the
drag from the DNS has an error of $0.76\%$ as shown in table
\ref{tab:AlignedSpheroidStokes-1}. The analytic data for the
torque is unavailable; therefore the difference of $3.5\%$, which
is between the DNS result and the simulation with the proposed
method in this paper, is listed in table
\ref{tab:AlignedSpheroidStokes-1}.

\begin{figure}[htbp]  
\begin{center}
\includegraphics[width=.45\textwidth]{fig1a} \quad %SpheroidStokes-u1x1.eps
\includegraphics[width=.45\textwidth]{fig1b} \\ SpheroidStokes-u2x1.eps
(a) \hspace{55mm} (b)
\end{center}
\caption{Fluid velocity comparisons between the results from the DNS and
from the simulation with the FCM for a Stokes flow: (a) $u_1$ versus $x_1$;
(b) $u_2$ versus  $x_1$.}
\label{Stokes-u1u2-x1}
\end{figure}

\begin{figure}[htbp] %fig2
\begin{center}
\includegraphics[width=.45\textwidth]{fig2a} \quad %SpheroidStokes-u1x2.eps
\includegraphics[width=.45\textwidth]{fig2b} \\ %SpheroidStokes-u2x2.eps
(a) \hspace{55mm} (b)
\end{center}
\caption{Fluid velocity comparisons between the results from the DNS and
from the simulation with the FCM for a Stokes flow: (a) $u_1$ versus $x_2$;
(b) $u_2$ versus  $x_2$.}
\label{Stokes-u1u2-x2}
\end{figure}


Comparisons of the velocity distributions between the DNS and the
simulation using the proposed method (denoted as FCM in the
following figures) are made near the ellipsoid. Figure
\ref{Stokes-u1u2-x1} shows the fluid velocity $u_1$ and $u_2$
versus the streamwise distance in $x_1$ with $x_2=x_3=0$. Figure
\ref{Stokes-u1u2-x2} presents $u_1$ and $u_2$ versus the lateral
distance in $x_2$ with $x_3=0$. Since a virtue ellipsoid is
actually placed in the simulation, there is no DNS result inside
this ellipsoid in the range $-2\leq x_1/a_2\leq2$. Outside this
ellipsoid, good agreement is achieved.


\subsection{Finte Reynolds Number Flow}
\label{chap:FiniteReynoldsFlow}

At a finite Reynolds number the flow is convective, not only does
the off-centered ellipsoid subject to a drag, but also it induces
a lateral lift force. It is demonstrated in this paper that the
penalty method described in  (\ref{PenaltyF}) and
(\ref{PenaltyT}) is able to predict the sensitive lateral lift
force as well as the torque on the ellispod. The computed
restoring forces balance the drag ($D$) and the lift ($L$) on the
ellipsoid to keep it from translating; the computed restoring
torque balances the torque ($\tau$) induced from the fluid on the
ellipsoid to prevent it from rotating. The ellipsoid of the same
size is placed at the origin of a channel of a smaller size
$-14\leq x_1/a_2\leq 6$, $-10/3\leq x_2/a_2\leq 10$, and $-5\leq
x_3/a_2\leq 5$. A parabolic velocity profile
$u_1=(1+3x_2/10)(1-x_2/10)$ is specified at the inlet $x_1=-14$.
Therefore the approach velocity is $u_o=1$. Based on this approach
velocity, the length scale $2a_2$, and the kinematic viscosity
$\nu=1$, the particle Reynolds number is 2. Periodic boundary
conditions are specified at $x_3=\pm5$. Non-slip boundary
conditions are imposed at $x_2=-10/3$ and $x_2=10$. Since there is
no analytic solution for this problem, the comparison was made
between the results from the DNS and the simulation with the FCM.
The computational meshes used in both the FCM simulation and the
DNS are conforming, non-uniform and structured, with 1920 and 2220
hexahedral elements respectively. The scaled drag, lift, torque,
and their associated relative errors are tabulated in table
\ref{tab:AlignedSpheroid_NS}. The relatively larger error in the
lift is because the lift is proportional to the shear stress and
it is very sensitive to the spatial resolution.

\begin{table}
\begin{center}
\begin{tabular}{||c|c|c|c|c|c||}\hline\hline
Source&$D/(\mu a_2 u_o)$&$L/(\mu a_2 u_o)$&$\tau/(\mu a^2_2
u_{o})$ &$V_{1}/u_o$& $\Omega_3 a_2/u_o$ \\\hline\hline Simul.
&34.371& 3.078 &4.102 &$4.61\times10^{-4}$&
$1.47\times10^{-3}$\\\hline DNS    &34.231& 2.898 & 4.12 & 0 &
0\\\hline Error& 0.4\%& 6.2\%& 3.5\%& 0.05\%& 0.15\%\\\hline
\end{tabular}
\end{center}
\caption{Comparisons of the scaled drag, lift and torque from the
DNS  and the simulation with the FCM for an ellipsoid in a channel
flow at Reynolds number 2. The numbers in percentage show the
relative errors against the values from the DNS.}
\label{tab:AlignedSpheroid_NS}
\end{table}

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.45\textwidth]{fig3a} \quad % SpheroidNS-u1x1.eps
\includegraphics[width=.45\textwidth]{fig3b} \\ %SpheroidNS-px1.eps
(a) \hspace{55mm} (b)
\end{center}
\caption{Comparisons between the results from the DNS and from the
simulation with the FCM at particle $Re=2$: (a) streamwise fluid velocity
$u_1$ versus $x_1$; (b) normalized pressure $p$ versus  $x_1$.}
\label{NS-u1p-x1}
\end{figure}

The comparisons of the flow field details between the DNS and the
simulation with the FCM are presented in the following figures.
Figure \ref{NS-u1p-x1}(a) shows the streamwise fluid velocity
distribution along $x_1$ through the center of the ellipsoid.
Figure \ref{NS-u1p-x1}(b) demonstrates the pressure distribution
along $x_1$. Both plots are along the line $x_2=x_3=0$. Near the
tips of the ellipsoid where $x_1/a_2=\pm 2$, there is a slight
difference between the two results.


\begin{figure}[htbp]
\begin{center}
\includegraphics[width=.45\textwidth]{fig4a} \quad %SpheroidNS-u1x2.eps
\includegraphics[width=.45\textwidth]{fig4b} \\ %SpheroidNS-u2x2.eps
(a) \hspace{55mm} (b)
\end{center}
\caption{Fluid velocity comparisons between the results from the DNS and
from the simulation with the FCM at particle $Re=2$: (a) $u_1$ versus $x_2$;
(b) $u_2$ versus  $x_2$.}
\label{NS-u1u2-x2}
\end{figure}

Figure \ref{NS-u1u2-x2}(a) presents the streamwise fluid velocity
profiles $u_1$ versus the lateral distance $x_2$ between the two
side walls. The locations of the profiles are at the center of the
ellipsoid, tangential to it in both the upstream and the
downstream directions, and further away from the ellipsoid. All
the velocity profiles are within the plane $x_3=0$. Figure
\ref{NS-u1u2-x2}(b) plots the lateral velocity ($u_2$) profiles at
the same locations as in Figure \ref{NS-u1u2-x2}(a). There are
small errors around the surface of the ellipsoid. However, in the
locations outside the volume of the  ellipsoid, the DNS result and
the simulation with the FCM agree well. But inside the ellipsoid,
the comparison can not be made, because the simulation involves a
virtue ellipsoid which is actually filled with fluid while in the
DNS a real ellipsoid is placed in the channel.


\begin{figure}[ht]
\begin{center}
\includegraphics[width=.45\textwidth]{fig5a} \quad % SpheroidNS-contours.eps
\includegraphics[width=.45\textwidth]{fig5b} % SpheroidNS-meshes.eps
(a)\hspace{55mm} (b)
\end{center}
\caption{Comparisons in the plane $x_3=0$ through the center of the
ellipsoid: (a) contour lines of the streamwise fluid velocity ; 
(b)  computational meshes; Top - simulation with the FCM; Bottom - DNS.}
\label{NS-ContourMesh}
\end{figure}

To further compare the overall flow field characteristics between
the  simulation with the FCM and the DNS, figure
\ref{NS-ContourMesh}(a) presents the contour lines of the
streamwise fluid velocity field. The top plot is the simulation
with the FCM while the bottom one is the DNS. In the top plot the
fluid fills the volume nominally occupied by the ellipsoid and
forms a virtual ellipsoid. In the DNS plot, the non-slip boundary
condition is imposed on the surface of the ellipsoid, which is
impermeable to the fluid. Therefore contour lines are outside the
ellipsoid. It is shown in this figure that outside the ellipsoid,
the flow field from the simulations with the FCM is similar and
comparable to that from the DNS. However, the simulation with the
FCM is computationally less expensive than DNS.


Figure \ref{NS-ContourMesh}(b) shows the computational meshes used in the
simulation with the FCM and the DNS. To enhance resolution while maintaining
efficiency with less elements, finer elements are fitted together around this
ellipsoid in the DNS mesh. On the other hand, the top plot does not need to
have mesh refinement around the position where the ellipsoid is placed.
Therefore computational cost is much less than the DNS.


\subsection*{Conclusions}
%\label{section:Conclusion}

This paper proposed an efficient method to simulate the  flow with
a two-way force interaction between the ellisoid and the ambient
fluid. The merit of this method lies in the fact that the presence
of an immersed object in the flow field is described implicitly in
the formulation. There is no need to enhance the resolution around
the solid-fluid interface. Therefore this helps to significantly
reduce the size of the unknowns in the computation domain and
lower the computational cost. The proposed method is capable of
providing good accuracy and remain robust in low Reynolds number
flow Regime due to the fact that the envelopes of the Gaussian
functions are determined \cite{LiuPhD04}. The spectral DNS on a
different mesh in this paper verifies this method. However, inside
the immersed ellipsoid, the flow field is ficticious and
comparisons can not be made. This method is good for predicting
both the hydrodynamical parameters and the flow field detail at a
relative less cpu time.

This method can be easily applied to the situation involving many immersed
spheres and/or ellipsoids while the computational overhead is only a small
fraction (depends on the number of immersed objects) of the cost for
computing the fluid phase.


\subsection*{Acknowledgements}
The authors would like to acknowledge NSF and DARPA for their
funding support on this research. The acknowledgement also goes to
Dr. Sune Lomholt at Danske Bank Copenhagen Area, Denmark, for
providing valuable suggestion and advice.

\begin{thebibliography}{0}

\bibitem{HappelB86}
J.~Happel and H.~Brenner.
\newblock {\em Low Reynolds Number Hydrodynamics}.
\newblock Martinus Nijhoff publishers, 1986.

\bibitem{Hu96}
H.H. Hu.
\newblock {Direct simulation of flows of solid-liquid mixtures}.
\newblock {\em Int. J. Multiphase Flow}, 22:335--352, 1996.

\bibitem{KarniadakisS1999}
G.E. Karniadakis and S.J. Sherwin.
\newblock {\em Spectral/hp Element Methods for CFD}.
\newblock Oxford University Press, 1999.

\bibitem{LiuPhD04}
D.~Liu.
\newblock {\em {Spectral Element$/$Force Coupling Method: Application to
  Colloidal Micro-Devices and Self-Assembled Particle Structures in 3D
  Domains.}}
\newblock PhD thesis, Brown University, 2004.

\bibitem{LiuMK02}
D.~Liu, M.R. Maxey, and G.E. Karniadakis.
\newblock A fast method for particulate microflows.
\newblock {\em J. Microelectromechanical Systems}, 11:691--702, 2002.

\bibitem{LiuMK04}
D.~Liu, M.R. Maxey, and G.E. Karniadakis.
\newblock A fast method for particulate microflows.
\newblock {\em J. Micromechanics and Microengineering}, 14:567--575, 2004.

\bibitem{Lomholt01}
S.~Lomholt.
\newblock {\em {Numerical investigations of macroscopic particle dynamics in
  microflows}}.
\newblock PhD thesis, Technical University of Denmark, 2001.

\bibitem{MaxeyP01}
M.R. Maxey and B.K. Patel.
\newblock Localized force representations for particles sedimenting in stokes
  flow.
\newblock {\em Int. J. Multiphase Flow}, 27:1603--1626, 2001.

\bibitem{MaxeyPCW97}
M.R. Maxey, B.K. Patel, E.J. Chang, and L.P. Wang.
\newblock {Simulation of Dispersed Turbulent Multiphase Flow}.
\newblock {\em Fluid Dynamic Research}, 20:143--156, 1997.

\end{thebibliography}


\end{document}