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