\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}