\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. 13--31.\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}{13}
\title[\hfilneg EJDE-2009/Conf/17/\hfil A computational domain decompositio]
{A computational domain decomposition approach for solving
coupled flow-structure-thermal interaction problems}
\author[E. Aulisa, S. Manservisi, P. Seshaiyer\hfil EJDE/Conf/17 \hfilneg]
{Eugenio Aulisa, Sandro Manservisi, Padmanabhan Seshaiyer} % in alphabetical order
\address{Eugenio Aulisa \newline
Mathematics and Statistics, Texas Tech University, Lubbock,
TX 79409, USA}
\email{eugenio.aulisa@ttu.edu}
\address{Sandro Manservisi \newline
DIENCA-Lab. di Montecuccolino, Via dei Colli 16,
40136 Bologna, Italy}
\email{sandro.manservisi@mail.ing.unibo.it}
\address{Padmanabhan Seshaiyer \newline
Mathematical Sciences, George Mason University, Fairfax, VA 22030, USA}
\email{pseshaiy@gmu.edu}
\thanks{Published April 15, 2009.}
\thanks{Supported by grants DMS 0813825 from the
National Science Foundation, and \hfill\break\indent
ARP 0212-44-C399 from the Texas Higher Education Coordinating Board.}
\subjclass[2000]{65N30, 65N15}
\keywords{Fluid-structure-thermal interaction;
domain decomposition; \hfill\break\indent multigrid solver}
\begin{abstract}
Solving complex coupled processes involving
fluid-structure-ther\-mal interactions is a challenging problem in
computational sciences and engineering. Currently there exist
numerous public-domain and commercial codes available in the area
of Computational Fluid Dynamics (CFD), Computational Structural
Dynamics (CSD) and Computational Thermodynamics (CTD). Different
groups specializing in modelling individual process such as CSD,
CFD, CTD often come together to solve a complex coupled
application. Direct numerical simulation of the non-linear
equations for even the most simplified fluid-structure-thermal
interaction (FSTI) model depends on the convergence of iterative
solvers which in turn rely heavily on the properties of the
coupled system. The purpose of this paper is to introduce a
flexible multilevel algorithm with finite elements that can be used
to study a coupled FSTI. The method relies on decomposing
the complex global domain, into several local sub-domains,
solving smaller problems over these sub-domains and then gluing
back the local solution in an efficient and accurate fashion to
yield the global solution. Our numerical results suggest that
the proposed solution methodology is robust and reliable.
\end{abstract}
\maketitle
\numberwithin{equation}{section}
\allowdisplaybreaks
\section{Introduction}
Engineering analysis is constantly evolving with a goal to develop
novel techniques to solve coupled processes that arise in
multi-physics applications. The efficient solution of a complex
coupled system which involves FSTI is still a challenging problem
in computational mathematical sciences. The solution of the
coupled system provides predictive capability in studying complex
nonlinear interactions that arise in several applications. Some
examples include a hypersonic flight, where the structural
deformation due to the aerodynamics and thermal loads leads to a
significant flow field variation or MAVs (Micro Air Vehicles)
where geometric changes possibly due to thermal effects may lead
to a transient phase in which the structure and the flow field
interact in a highly non-linear fashion.
The direct numerical simulation of this highly non-linear system,
governing even the most simplified FSTI, depends on the
convergence of iterative solvers which in turn relies on the
characteristics of the coupled system. Domain decomposition
techniques with non-matching grids have become increasingly
popular in this regard for obtaining fast and accurate solutions
of problems involving coupled processes. The mortar finite element
method \cite{BMP1,Belgacem} has been considered to be a viable
domain decomposition technique that allows coupling of different
subdomains with nonmatching grids and different discretization
techniques. The method has been shown to be stable mathematically
and has been successfully applied to a variety of engineering
applications \cite{SS2,BCP}. The basic idea is to replace the
strong continuity condition at the interfaces between the
different subdomains by a weaker one to solve the problem in a
coupled fashion. In the last few years, mortar finite element
methods have also been developed in conjunction with multigrid
techniques, \cite{AMS1,AMS2,BDW,GP}. One of the great advantages
of the multigrid approach is in the grid generation process
wherein the corresponding refinements are already available and no
new mesh structures are required. Also, the multigrid method
relies only on local relaxation over elements and the solution on
different domains can be easily implemented over parallel
architectures.
The purpose of this paper is to introduce a flexible multigrid
algorithm that can be used to study different physical processes
over different subdomains involving non-matching grids with less
computational effort. In particular, we develop the method for a
model problem that involves Fluid-Structure-Thermal interaction (FSTI).
In section 2, the equations of the coupled model are discretized via
the finite element discretization. In section 3, the multigrid
domain decomposition algorithm to solve the discrete problem is discussed.
Finally in section 4 the method is applied to a two-dimensional
FSTI application.
\section{Model and governing equations}
Let the computational domain $\Omega \subset \Re^2 $ be an open set
with boundary $\Gamma$.
Let the fluid subdomain $\Omega_f$
and the solid subdomain $\Omega_s$
be two disjoint open sets
with boundary $\Gamma_f$ and $\Gamma_s$,
respectively and let $\Omega=\bar{\Omega}_f \cup \bar{\Omega}_s$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig1a}
\includegraphics[width=0.45\textwidth]{fig1b}
\end{center}
\caption{Domain $\Omega=\Omega_f \cup \Omega_s $ in two different
configurations.} \label{fig1}
\end{figure}
Figure \ref{fig1} presents illustrations of two sample
computational domains. $\Gamma_{sf}$ is the interior boundary
between $\Omega_f$ and $\Omega_s$, $\Gamma_f^e=\Gamma \cap
\Gamma_f$ is the fluid exterior boundary and $\Gamma_s^e=\Gamma
\cap \Gamma_s$ is the solid exterior boundary. For simplicity let
us assume that the only boundary which can change in time is the
interior boundary $\Gamma_{sf}$. In agreement with this assumption
both subdomains $\Omega_f$ and $\Omega_s$ are time dependent and
constrained by
\begin{equation}
\bar{\Omega}_f(t) \cup
\bar{\Omega}_s(t)=\bar{\Omega} \,. \label{DOMEGA_DT}
\end{equation}
In the model problem, we employ the unsteady Navier-Stokes equations
describing the flow of a fluid in a region $\Omega_f$ given by:
\begin{gather*}
\rho_f \frac{\partial \vec {u}}{\partial t} - \mu_f \Delta \vec{u}
+ \rho_f (\vec{u} \cdot
\nabla) \vec{u} + \nabla p = \vec{f} \quad\text{in }\Omega_f \times
(0, T) \\
\nabla \cdot \vec{u} = 0 \quad\text{in }\Omega_f \times (0, T) \\
\vec{u} = \vec{U} \quad\text{on }\Gamma_1
\end{gather*}
where $\rho_f$ and $\mu_f$ are the density and the viscosity and
$ \vec{f}$ is the body force.
This is coupled with the energy equation given by,
\begin{gather*}
\rho c_p \frac{\partial T}{\partial t} - k \Delta T
+ \rho c_p (\vec{u} \cdot \nabla T) = 0 \quad\text{in }\Omega \times (0, T) \\
T = \Theta \quad\text{on }\Gamma_2
\end{gather*}
that is solved over the whole domain
$\Omega $.
In the solid region $\Omega_s $
the approximate Euler-Bernoulli beam equation is considered.
In this approximation
plane cross sections perpendicular to the axis of the beam are
assumed to remain plane and perpendicular
to the axis after deformation \cite{REDDY} and
under these hypotheses only a mono-dimensional
model is required for the normal transverse deflection field $w$.
We will denote by $\Lambda$ the beam axis
and by $(\xi,\eta)$ a local reference system
oriented with the $\xi$-axis parallel to $\Lambda$.
\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.6\textwidth]{fig2}
\end{center}
\caption{Domain notation for the beam domain $\Omega_s $.}
\label{fig2}
\end{figure}
As shown in Figure \ref{fig2}, variables
$\delta$ and $ L $ are the thickness and the length of the beam
respectively,
the interior boundary $\Gamma_{sf}$
is in $(\xi,\pm\delta /2) $ for $ 0 \le \xi \le L $ and in
$ (L,\eta) $ for $-\delta/2 \le \eta \le \delta / 2 $.
In $\Gamma_1 \subset \Gamma_f^e$,
Dirichlet boundary conditions are imposed for the velocity field $\vec{u}$;
Neumann homogenous boundary conditions are considered on the
remaining part,
$\Gamma_f^e \setminus \Gamma_1$.
Similarly, on $\Gamma_2 \subset \Gamma$
Dirichlet boundary conditions are imposed for the temperature $T$,
while Neumann homogenous boundary conditions are considered on $\Gamma \setminus \Gamma_2$.
In $\xi=0$ Dirichlet zero boundary conditions are imposed for the solid displacement $w$ and its derivatives.
Conditions of displacement compatibility and force equilibrium
along the structure-fluid interface $\Gamma_{sf}$ are satisfied.
Let $ \vec{U} \in \mathbf{H}^{1/2}(\Gamma_1) $
be the prescribed boundary velocity over $ \Gamma_1 $,
satisfying the compatibility condition, and
$\Theta \in {H}^{1/2}(\Gamma_2) $
be the prescribed temperature over $ \Gamma_2 $.
We are using ${H}^k(\Omega)$ to denote the space of functions with $k$
generalized derivatives. We set $L^2(\Omega)={H}^0(\Omega)$ and note that
the derivation of these spaces can be extended to non-integer values $k$
by interpolation.
The velocity, the pressure, the temperature and the beam deflection
$(\vec{u},p,T,w) \in \mathbf{H}^1(\Omega_f) \times L^2(\Omega_f)\times
{H}^{1}(\Omega)\times H^2(\Lambda)$
satisfy the weak variational form of the unsteady fully coupled system given by
the Navier-Stokes system over $\Omega_f$,
\begin{gather}
\begin{aligned}
&\langle \rho_f\,\frac{\partial \vec{u}}{\partial t},\vec{v}_f\rangle+
a_f\big(\vec{u},\vec{v}_f\big)+
b_f\big(\vec{v}_f,p\big)+
c_f\big(\vec{u};\vec{u},\vec{v}_f\big) \\
&=-\langle \rho_f\,\vec{g}\,\beta(T-T_0),\vec{v}_f\rangle \quad
\forall \vec{v}_f \in \mathbf{H}^1(\Omega_f)\,,
\end{aligned} \label{NS}
\\
b_f\big(\vec{u},r_f\big) = 0 \quad \forall r_f \in L^2(\Omega_f)\,,
\label{CNT} \\
\langle \vec{u} - \vec{U},\vec{s}_f\rangle _{\Gamma_{1}}=0 \quad
\forall \vec{s}_f \in \mathbf{H}^{-\frac{1}{2}}(\Gamma_1)\,, \label{BC1}
\\
\langle \vec{u} - \dot{w} \cdot \hat{n}_{sf},\vec{s}_{sf}\rangle_{\Gamma_{sf}}=0
\quad \forall \vec{s}_{sf} \in \mathbf{H}^{-\frac{1}{2}}(\Gamma_{sf})\,,
\label{BC4}
\end{gather}
the energy equation over $\Omega$
\begin{gather}
\langle \rho\, c_p\,\frac{\partial T}{\partial t},v\rangle +
a\big(T,v\big)+
c\big(\vec{u};T,v\big)= 0 \quad \forall v \in {H}^{1}(\Omega)\,,
\label{TMP} \\
\langle T - \Theta,s\rangle _{\Gamma_{2}}=0 \quad
\forall s \in {H}^{-\frac{1}{2}}(\Gamma_2)\,, \label{BC2}
\end{gather}
and the Euler-Bernoulli beam equation over $\Omega_s$,
\begin{gather}
\langle \rho_s\, \delta\,\ddot{w},v_s\rangle +
a_s\big(w,v_s\big)=
\langle p(\vec{x}(\xi,-\frac{\delta}{2}),t)
-p(\vec{x}(\xi,\frac{\delta}{2}),t),v_s \rangle \quad
\forall v_s, \in H^2(\Lambda)\,,\label{DSP} \\
w(0,t)=0\,,\quad \frac{\partial w(0,t)}{\partial \xi}=0\,, \quad
\dot{w}(0,t)=0\,. \label{BC3}
\end{gather}
In (\ref{NS})-(\ref{CNT}) the continuous bilinear forms are defined as
\begin{gather}\label{NSFORM1}
a_f ( \vec{u}, \vec{v} )=\int_{\Omega_f}2\mu_f\,D( \vec{u} )
:D(\vec{v} )\,d \vec{x}
\quad \forall \,\vec{u} ,\vec{v} \in \mathbf{H}^1(\Omega_f),
\\ \label{NSFORM2}
b_f(\vec{v},r)=-\int_{\Omega_f} r \,\nabla \cdot \vec{v} \,d \vec{x}
\quad \forall\, r\in L^2(\Omega_f) \,,\; \forall\, \vec{v} \in\mathbf{H}^1(\Omega_f)
\end{gather}
and the trilinear form as
\begin{equation}\label{NSFORM3}
c_f(\vec{w};\vec{u},\vec{v} )
=\int_{ \Omega_f} \rho_f\,(\vec w\cdot\nabla) \vec u\cdot\vec v\,d\vec{x}
\quad \forall\, \vec{w}, \vec{u} ,\vec{v} \in \mathbf{H}^1(\Omega_f)\,,
\end{equation}
where $\rho_f$ and $\mu_f$ are the density and the viscosity of
the fluid. The distributed force in Eq. (\ref{NS})
is the Boussinesq approximation of the buoyancy force,
where $\vec{g}$ is the gravity acceleration, $\beta$ the
volumetric expansion coefficient of the fluid and $T_0$ a reference
temperature. For $T>T_0$ the fluid expands then
the density decreases and the buoyancy force points in the direction
opposite to the gravity.
When $T