\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<T_0$ both the buoyancy force and the gravity point in the same direction. In
 (\ref{TMP}) the bilinear form is defined as
\begin{equation}\label{TMPFORM1}
a(T,v)=\int_{\Omega} k\, \nabla T \cdot \nabla v \,d \vec{x}
\quad \forall \,T ,v \in { H}^1(\Omega)
\end{equation}
and the trilinear form
\begin{equation}\label{TMPFORM2}
c(\vec{u};T,v)
=\int_{\Omega}\rho\,c_p (\vec u\cdot\nabla T) v\,d\vec{x}
\quad \forall\, \vec{u} \in \mathbf{H}^1(\Omega), \; T , v
\in {H}^1(\Omega)\,,
\end{equation}
where $\rho$, $c_p$, and $k$ are the density, the heat capacity
and the  heat conductivity, respectively. If the integral is over
the subdomain $\Omega_f$, the fluid physical properties $\rho_f$,
$c_{p_f}$, and $k_f$ are used, otherwise over $\Omega_s$ the solid
properties $\rho_s$, $c_{p_s}$, and $k_s$ are used. Furthermore in
the solid region the trilinear form $c(\vec{u},T,v)$ is
identically zero, since the velocity $\vec{u}$ is zero. In
(\ref{DSP}) the bilinear form is defined as
\begin{equation}\label{DSPFORM1}
a_s(w,v)=
\int_{\Lambda} EI\,\frac{\partial^2 w}{\partial \xi^2}
\frac{\partial^2 v}{\partial \xi^2}\,d \vec{x}
\quad \forall \,w ,v \in { H}^2(\Lambda)\,,
\end{equation}
where $EI$ is the beam stiffness coefficient for unitary deepness.
In the right hand side of Eq. (\ref{DSP}) the load due to the
pressure difference between the two sides of the beam is given.
(\ref{DSP}) represents the force equilibrium constraint between
the two subdomains $\Omega_f$ and $\Omega_s$ on the common
boundary $\Gamma_{sf}$. For details concerning the function
spaces, the bilinear and the trilinear forms and their properties,
one may consult \cite{REDDY,GR1,TM}. Equations (\ref{BC1}),
(\ref{BC2}) and (\ref{BC3}) represent the exterior Dirichlet
boundary condition for the velocity, the temperature and the
displacement, respectively. Eq. (\ref{BC4}) represents the
compatibility constraints between the velocity field $\vec{u}$ and
the time derivative of the beam deflection $w$ on $\Gamma_{sf}$.
The unitary normal $\hat{n}_{sf}$ points in the same direction as
the local $\eta$-axis. According to Eqs. (\ref{DOMEGA_DT}) and
(\ref{BC4}), changes in  the fluid and solid subdomains $\Omega_f$
and $\Omega_s$  should be also considered.

\subsection{Domain decomposition}

Let us now introduce a non-conforming formulation of the problem.
Let the domain $\Omega$ be partitioned into $m$ non-overlapping
sub-domains $\{\Omega^i\}_{i=1}^{m}$ such that
$\partial\Omega^i\cap\partial\Omega^j $ ($i\ne j$) is either
empty, a vertex, or a collection of edges of $\Omega^i$ and
$\Omega^j$. In the latter case, we denote this interface by
$\Gamma^{ij}$ which consists of individual common edges from the
domains  $\Omega^i$ and $\Omega^j$. Let now the fluid domain
$\Omega_f$ be partitioned into $m$ non-overlapping sub-domains
$\{\Omega^i_f\}_{i=1}^{m}$, where $\Omega_f^i$ is given by
$\Omega_i \cap \Omega_f$. The fluid partition
$\{\Omega^i_f\}_{i=1}^{m}$ is obtained from the domain partition
$\{\Omega^i\}_{i=1}^{m}$, subtracting the solid region from each
subdomain $\Omega^i$, then $\Omega_f^i$ is an empty region if
$\Omega^i$ is a subset of $\Omega_s$. The common boundary between
two subregions $\Omega_f^i$ and $\Omega_f^j$  is denoted by
$\Gamma_f^{ij}$. The velocity, the pressure,  the stress vector,
the temperature, the heat flux and the displacement $ (\vec{u}^{
i}, p^{i},\vec{\tau}^{ij}, T^{i}, q^{ij}, w) \in
\mathbf{H}^1(\Omega_f^{i}) \times L^2(\Omega_f^{i}) \times
\mathbf{H}^{-1/2}(\Gamma_f^{ij}) \times H^1(\Omega^{i}) \times
H^{-1/2}(\Gamma^{ij}) \times  H^2(\Lambda)$ satisfy the following
system of equations
\begin{gather}
\begin{aligned}
&\langle \rho_f\,\frac{\partial \vec{u}^{i}}{\partial
t},\vec{v}^{i}_f\rangle + a_f\big(\vec{u}^{i},\vec{v}^{i}_f\big)+
b_f\big(\vec{v}^{i}_f,p^{i}\big)+
c_f\big(\vec{u}^{i};\vec{u}^{i},\vec{v}^{i}_f\big) \\
&+\langle \vec{\tau}^{ij} , \vec{v}_f^{i}\rangle_{\Gamma_f^{ij}} = -\langle
\rho_f\,\vec{g}\,\beta(T^{i}-T_0),\vec{v}_f^{i}\rangle
 \quad \forall \vec{v}_f^{i}\in \mathbf{H}^{1}(\Omega^i_f )\,,
 \end{aligned}  \label{SYS2A}\\
b_f\big(\vec{u}^{i},r_f^{i}\big) =  0 \quad
\forall r^{i}_f \in L^2 (\Omega^i_f )\,,\\
\langle \vec{u}^{i} - \vec{U},\vec{s}^{i}_f\rangle _{\Gamma^{i}_1}=0 \quad
\forall \vec{s}^{i}_f    \in \mathbf{H}^{-1/2}(\Gamma_1^i)\,,
\\
\langle \vec{u}^{i}-  \vec{u}^{j},\vec{s}^{ij}_f \rangle _{\Gamma^{ij}_f} =
0 \quad \forall \vec{s}^{ij}_f    \in \mathbf{H}^{-1/2}(\Gamma^{ij}_f)\,,  \\
\langle \vec{u}^{i} - \dot{w} \cdot
\hat{n}_{sf},\vec{s}^{i}_{sf}\rangle _{\Gamma^{i}_{sf}}=0  \quad \forall
\vec{s}^{i}_{sf}  \in \mathbf{H}^{-1/2}(\Gamma^{i}_{sf})\,,
\\
\langle \rho\, c_p\,\frac{\partial T^{i}}{\partial t},v^{i}\rangle +
a\big(T^{i},v^{i}\big)+
c\big(\vec{u}^{i};T^{i},v^{i}\big) % \\
+\langle \vec{q}^{ij} , \vec{v}^{i}\rangle_{\Gamma^{ij}}= 0   \quad
\forall v^{i}    \in {H}^{1}(\Omega^i )\,, \label{SYS2B}
\\
\langle T^{i} - \Theta,s^{i}\rangle_{\Gamma_2^{i}}=0   \quad
\forall s^{i}    \in {H}^{-1/2}(\Gamma_2^i)\,,
\\
\langle T^{i}-  T^{j},\vec{s}^{ij} \rangle_{\Gamma^{ij}} \,  = 0 \quad
 \forall s^{ij}   \in {H}^{-1/2}(\Gamma^{ij})\,,
\\
\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{SYS2C}
\\
w(0,t)=0\,,\quad \frac{\partial w(0,t)}{\partial \xi}=0\,,
\quad \dot{w}(0,t)=0\,,   \label{AAAA}
\end{gather}
for  $i=1,2,\dots ,m$  where
$\Gamma_1^i=\Gamma_1 \cap \partial{\Omega^{i}_f}$,
 $\Gamma_2^i=\Gamma_2 \cap \partial{\Omega^{i}}$,
the stress vector $ \vec{\tau}^{ij}=-\mu \nabla \vec{u}^{i} \cdot
\hat{n}^{ij}_f + p^{i}\,\hat{n}^{ij}_f$ and the heat flux $
q^{ij}= -k \nabla T^{i} \cdot \hat{n}^{ij}$, with ${n}^{ij}_f$ and
${n}^{ij}$  the unitary external vectors normal to  the subdomains
$\Omega^{i}_f$ and $\Omega^{i}$, respectively.

Note that in the continuous case, on the boundary $\Gamma^{ij}_f$,
the velocity vectors, $\vec{u}^{i}$ and $\vec{u}^{j}$ are in the
same space $\mathbf{H}^{1/2}( \Gamma^{ij}_f)$, namely   we have
$\vec{u}^{i}=\vec{u}^{j}$ pointwise. The same remark can be
applied to
 the stress vectors $\vec{\tau}^{ij}$ and $\vec{\tau}^{ji}$,
the temperatures $T^{i}$ and $T^{j}$ and the heat fluxes $q^{ij}$
and $q^{ji}$. In the rest of the paper, the variables $\vec{u}$,
$p$ and $T$, without the label $^i$, should be considered as a
collection of all the local variables $\vec{u}^{i}$, $p^{i}$ and
$T^{i}$, where
\begin{gather*}
\vec{u}(\vec{x},t)= \vec{u}^{i}(\vec{x},t) \quad
 \forall\;\vec{x}\in\Omega_f^{i}\,, \\
p(\vec{x},t)= p^{i}(\vec{x},t) \quad \forall\;\vec{x}\in\Omega_f^{i}\,,  \\
T(\vec{x},t)= T^{i}(\vec{x},t) \quad \forall\;\vec{x}\in\Omega^{i}\,,
\end{gather*}
for $i=1,2,\dots ,m$. It is straightforward to prove that the system
(\ref{SYS2A})-(\ref{AAAA}) for the local state variables $
(\vec{u}^{i}, p^{i},\vec{\tau}^{ij}, T^{i}, q^{ij}, w) $, for
$i=1,2,\dots ,m$,  implies the  system (\ref{NS})-(\ref{BC3}) for the
global state variables $ (\vec{u}, p,\vec{\tau}, T, q, w) $ .

\subsection{ALE formulation and time discretization}

In order to account for the changing nature of the fluid and solid
subdomains, we wish to define a dynamic mesh for the space
discretization. However, to avoid extreme distortion, we choose to
move the mesh independently of the fluid velocity in the interior
of $\Omega_f$. Such a scheme, called arbitrary Lagrangian-Eulerian
(ALE) formulation, is commonly applied when studying
fluid-structure interaction \cite{SP,DGH,GGM,HLZ}. Inside the
solid region each point is moving according to the time derivative
of the displacement $w$. The grid velocity $\vec{u}_g$  can be any
velocity satisfying the following constraints
\begin{gather}
\vec{u}_g=0 \quad\text{on } \Gamma\,, \label{LAEU1}\\
\vec{u}_g=\vec{u}\quad\text{on } \Gamma_{sf}\,, \label{LAEU2}\\
\vec{u}_g(\vec{x}(\xi,\eta),t)=\dot{w}(\xi,t) \, \hat{n}_{sf}
\label{LAEU3}\quad\text{in } \Omega_s \,.
\end{gather}
If the grid velocity is known as a function of time, the
trajectory inside the domain $\Omega$ of a generic point of
coordinate $\vec{x}(t)$ can be traced by solving the integral
equation
\begin{equation}
\vec{x}(t_1)=\vec{x}(t_0)+\int_{t_0}^{t_1} \vec{u}_g(\vec{x}(t'),t')\, dt'\,.
\end{equation}
The Lagrangian derivative of the velocity field
$\vec{u}(\vec{x}(t),t)$, evaluated along the point trajectory, is
given by
\begin{equation}
\begin{aligned}
\frac{D\,\vec{u}(\vec{x}(t),t)}{D\,t}
&=\frac{\partial\vec{u}}{\partial t}+
\frac{\partial\vec{u}}{\partial x_1}\,\frac{d{x}_1}{dt}+
\frac{\partial\vec{u}}{\partial x_2}\,\frac{d{x}_2}{dt} \\
&=\frac{\partial\vec{u}}{\partial t}+
\frac{\partial\vec{u}}{\partial x_1}\,u_{g_1}+
\frac{\partial\vec{u}}{\partial x_2}\,u_{g_2}\\
&= \frac{\partial\vec{u}}{\partial t}+
(\vec{u}_g\cdot \nabla)\, \vec{u}\,.
\end{aligned} \label{LAGU}
\end{equation}
Then the Eulerian derivative can be written as the difference between
the Lagrangian derivative and the corresponding grid velocity advection term
\begin{equation}
\frac{\partial\vec{u}}{\partial t}=
\frac{D\,\vec{u}(\vec{x}(t),t)}{D\,t}-
(\vec{u}_g\cdot \nabla)\, \vec{u}\,. \label{LAGU1}
\end{equation}
In the same way the Eulerian derivative of the temperature $T$ can be expressed as
the difference between the Lagrangian derivative and the corresponding grid velocity advection term
\begin{equation}
\frac{\partial T}{\partial t}=
\frac{D\,T(\vec{x}(t),t)}{D\,t}-
(\vec{u_{g}}\cdot \nabla) T\,.\label{LAGT}
\end{equation}
Using  (\ref{LAGU1})-(\ref{LAGT}) in system (\ref{SYS2A})-(\ref{AAAA}),
the Navier-Stokes and energy conservation equations become
\begin{gather}
\begin{aligned}
&\langle \rho_f\,\frac{D\,\vec{u}^{i}(\vec{x}(t),t)}{D\,t},\vec{v}^{i}_f\rangle
+ a_f\big(\vec{u}^{i},\vec{v}^{i}_f\big)+
b_f\big(\vec{v}^{i}_f,p^{i}\big)
+ c_f\big(\vec{u}^{i}-\vec{u}^{i}_g;\vec{u}^{i},\vec{v}^{i}_f\big)
+\langle \vec{\tau}^{ij} , \vec{v}_f^{i}\rangle_{\Gamma_f^{ij}}\\
&= -\langle \rho_f\,\vec{g}\,\beta(T^{i}-T_0),\vec{v}_f^{i}\rangle
\end{aligned}
\\
\langle \rho\, c_p\,\frac{D\,T^{i}(\vec{x}(t),t)}{D\,t},v^{i}\rangle
+a\big(T^{i},v^{i}\big)+
c\big(\vec{u}^{i}-\vec{u}^{i}_g;T^{i},v^{i}\big) +\rangle \vec{q}^{ij} ,
\vec{v}^{i}_f\rangle_{\Gamma^{ij}}= 0 \,. \nonumber
\end{gather}
Given the initial velocity field $\vec{u}_0$ the Lagrangian derivatives
can be  discretized in time using a simple first order integration method.

The structural equation is discretized in time by a using Newmark
integration scheme \cite{REDDY}. In this method the displacement and its time
derivative are approximated according to
\begin{gather}
w_{t+1}=w_t+\Delta t\, \dot{w}_t +\frac{1}{2}\Delta t^2\,
\ddot{w}_{\,t+\gamma}\,, \label{NM1}\\
\dot{w}_{t+1}=\dot{w}_t+\Delta t\,\ddot{w}_{\,t+\alpha}\,, \label{NM2}
\end{gather}
where
\begin{equation}
\ddot{w}_{\,t+\theta} = (1-\theta)\ddot{w}_t+\theta\ddot{w}_{t+1}\,,
 \label{NM3}
\end{equation}
and $\alpha$ and $\gamma$ are parameters that determine the
stability  and accuracy of the method. We chose
$\alpha=\gamma=0.5$, which is known as the constant-average
acceleration method. This technique is stable for each time step
$\Delta t$ and conserves energy for free vibration problem. The
use of  (\ref{NM1})-(\ref{NM3}) in  (\ref{SYS2C}) gives the
following approximation
\begin{equation}
\begin{gathered}
\begin{aligned}
&\langle \rho_s\, \delta\,(a_3\,w_{t+1}),v_s\rangle +
a_s\big(w_{t+1},v_s\big)
 -\langle (p_{t+1}(\vec{x}(\xi,-\frac{\delta}{2}),t)-
    p_{t+1}(\vec{x}(\xi,\frac{\delta}{2}),t),v_s \rangle \\
&=\langle \rho_s\, \delta  (a_3\,w_t + a_4\,\dot{w}_t
+a_5\,\ddot{w}_t),v_s\rangle \,,
\end{aligned} \label{DSP2}
\\
a_3=\frac{2}{\gamma\,\Delta t^2},\quad
a_4=\frac{2}{\gamma\,\Delta t},\quad
a_5=\frac{1}{\gamma}-1\,,
\end{gathered}
\end{equation}
where all the terms in the right hand side, $w_t$, $\dot{w}_t$ and
$\ddot{w}_t$, are evaluated at the previous time step. Note that
the calculation of the right hand side requires knowledge of the
initial conditions $w_0$, $\dot{w}_0$ and $\ddot{w}_0$.
In practice, one does not have $\ddot{w}_0$.
As an approximation, it can be calculated from  (\ref{SYS2C})
\[
\langle \rho_s\, \delta\,\ddot{w}_0,v_s\rangle
=\langle p_0(\vec{x}(\xi,-\frac{\delta}{2}),t)-p_0(\vec{x}
(\xi,\frac{\delta}{2}),t),v_s \rangle
-a_s\big(w_0,v_s\big)\,,
\]
for given initial displacement $w_0$ and initial pressure $p_0$.
At the end of each time step, the new velocity
$\dot{w}_{t+1}$ and acceleration $\ddot{w}_{t+1}$ are computed using
\begin{equation} \label{ACCVEL}
\begin{aligned}
 \ddot{w}_{t+1}=a_3(w_{t+1}-w_s)-a_4\,\dot{w}_t -a_5\,\ddot{w}_t \,,\\
 \dot{w}_{t+1} = \dot{w}_t + a_2\,\ddot{w}_t + a_1\,\ddot{w}_{t_1}\,,\\
 a_1=\alpha\,\Delta t,\; a_2=(1-\alpha)\,\Delta t\;.
\end{aligned}
\end{equation}

\section{System discretization}
\subsection{Multilevel domain decomposition}

In this section a multilevel domain decomposition methodology
 will be described for the
whole domain $\Omega$ and the corresponding variable $T$.
 The same consideration can be directly
extended to the fluid domain $\Omega_f $ and the corresponding
variables $\vec{u}$ and $p$.

Let us introduce a finite element discretization in each subdomain
$\Omega^i$ through the mesh parameter $h$ which tends to zero. Let
$\{\Omega_h^{i}\}_{i=1}^m$ be the partition of the discretized
domain $ \Omega_h $. Now, by starting  at the multigrid coarse
level $ l=0 $, we subdivide each $ \Omega^i_h $ and consequently $
\Omega_h $ into triangles or rectangles by families of meshes $
\mathcal{T}^{i, 0}_h $. Based on a simple element midpoint
refinement different multigrid levels can be built to reach the
finite element meshes $ \mathcal{T}^{i, n}_h$ at the top finest
multigrid level $l=n$. At the coarse level, as at the generic
multigrid level $l$, the triangulation over two adjacent
subdomains, $ \Omega^i_h $ and  $ \Omega^j_h $, obeys the finite
element compatibility constraints along the common interfaces
$\Gamma^{ij}_h$. For details on multigrid levels and their
construction one may consult \cite{BR,ST,T1}.

By using this methodology  we have constructed
a sequence of meshes for each
multigrid level in a standard finite element fashion
with compatibility  enforced across all the element
interfaces built over midpoint refinements.
In every subdomain $ \Omega^i_h $ the energy equations
can be solved over a different level mesh,
generating a global solution over $ \Omega_h $,
consisting mesh solutions at different levels over different subdomains.
Let  $ \Omega^{i, l}_h $  be the subdomain $i$ where the
solution will be computed at the multigrid level  $l$.
It should be noted that the multigrid levels at which the solution is
computed over adjacent subdomains,
 $\Omega_h^{i, l}$ and  $\Omega_h^{j, k}, $ may be
different from each other ($l\ne k$), with no compatibility
enforced across the common interface $ \Gamma^{ij}_h$.


\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig3}
\end{center}
\caption{Domain (a), domain partitioning (b),conforming coarse level mesh (c),
midpoint refinement generation of two finer conforming level meshes (d),(e),
non-conforming mesh (f).} \label{fig3}
\end{figure}

In Figure \ref{fig3} an example of non-conforming mesh partitioning
is given for a square domain $\Omega_h$ (Figure \ref{fig3}.a).
The square is divided in
four subdomains (Figure \ref{fig3}.b) and 3 conforming
level meshes are considered ($l=0,1,2$).
The coarse mesh (Figure \ref{fig3}.c) is given, while the other
two (Figures 3.d-e) are generated
by successive  midpoint refinements. In Figure \ref{fig3}.f one possible
non-conforming mesh configuration is illustrated.
Note that in Figure \ref{fig3}.f only across the common boundary
$\Gamma^{1,2}_h$ finite element compatibility is
obtained since the two subdomain meshes $\Omega^{1,0}_h$
and $\Omega^{2,0}_h$ are chosen at the same level ($l=0$).
On the other three common boundaries,
$\Gamma^{1,3}_h$, $\Gamma^{2,,4}_h$ and
$\Gamma^{3,4}_h$,  compatibility is not enforced
since different level meshes are considered between
 two adjacent subdomains.
It should be noticed that on the
common boundaries the nodes of the coarsest grid
are always included in the nodes of
the finest one. This is  the key point for
discretization and resolution
of  (\ref{SYS2A})-(\ref{AAAA}) and
it is always true
for different level meshes generated
by using midpoint refinements.

Finite element approximation spaces can be generated regularly,
as function of the characteristic length $h$ over each multigrid level $l$
resulting in different approximation spaces over the solution meshes  $ \Omega^{i, l}_h $.
In the rest of the paper we denote with labels
$  i,l $ the solution over the corresponding
subdomains, i.e., for the temperature $T^{i,l}$,
and with no labels  the extended solution over
$\Omega $, i.e. $ T_h $ for the extended temperature.
Note that the temperature $T^{i,l}$  is computed over each
$\Omega^{i, l}_{h}$ at the  corresponding level $l$,
but the extended solution $T_h$ on the top level $n$ is
defined over all $ \Omega_h $ in a standard and regular way.
There may be parts of the domain where the solution is not computed
at the top level but a projection operator $\mathcal{I}_l^n$ from the
coarse level $l$ to the top level $n$ can always be used to
interpolate the solution over $\Omega_h$.
The extended temperature is therefore defined by
\[
T_h(\vec{x},t)=\mathcal{I}_l^n\,T^{i,l}(\vec{x},t)\,,
\]
for all $\vec{x} \in \Omega^{i}_h$, $i=1,2,\dots ,m$. We can easily
generalize the notations and evaluate
the extended velocity and pressure
by using the same operator  $ \mathcal{I}_l^n $ as
\begin{gather}
\vec{u}_h(\vec{x},t)=\mathcal{I}_l^n\,\vec{u}^{i,l}(\vec{x},t),\\
 p_h(\vec{x},t)=\mathcal{I}_l^n\,p^{i, l}(\vec{x},t)\,,
\end{gather}
for all $\vec{x} \in \Omega^{i}_{h_f}$, $i=1,2,..,m$,
where $\{\Omega^{i}_{h_f}\}_{i=1}^m$ is the partition
over the discrete subdomains of $\Omega_{h_f}$.
These extended functions take  the same values over
the  coarse and the top mesh
 at those nodes included in both meshes.



\subsection{Non-conforming finite element discretization}

Let $\mathbf{X}^l_{h}(\Omega_f) \subset \mathbf{H}^1(\Omega_f)$,
$S^l_{h}(\Omega_f) \subset L^2(\Omega_f)$,
${\bf R}^l_{h}({\Gamma_f}) = \mathbf{X}^l_{h}|_{\Gamma_f}\subset\mathbf{H}^{1/2}(\Gamma_f)$,
${X}^l_{h}(\Omega) \subset {H}^1(\Omega)$ and
${R}^l_{h}({\Gamma}) = {X}^l_{h}|_{\Gamma}\subset{H}^{1/2}(\Gamma)$
be the approximation spaces.
At each level mesh $l$ we chose the families of finite element spaces
to satisfy appropriate  stability and approximation properties \cite{GH1,GM2}
that will allow us building a regular conforming approximation.
We indicate with $\mathbf{P}^l_{h}({\Gamma_f}) \subset\mathbf{H}^{-1/2}(\Gamma_f)$
the dual space of ${\bf R}^l_{h}({\Gamma_f})$ and with
${P}^l_{h}({\Gamma}) \subset{H}^{-1/2}(\Gamma)$ the dual space of ${R}^l_{h}({\Gamma})$.
The beam space discretization is obtain by choosing
the family of Hermitian elements,
$Z_h(\Lambda)\subset H^2(\Lambda)$,
for which the interpolation functions are continuous with non-zero derivatives
up to order two.

Let $J^i$ be the set of the $j$-indices of all the neighboring
regions $\Omega^j$ surrounding the subdomain $\Omega^{i}$ and
$J_f^i$ be the set of the $j$-indices of all the neighboring
regions $\Omega^j_f$ surrounding the subdomain $\Omega^{i}_f$. Let
$ (\vec{u}^{i,l}, p^{i,l},\vec{\tau}^{ij,l}, T^{i,l}, q^{ij,l},w)$
be in $ \mathbf{X}^l_h(\Omega_f^{i}) \times S_h^l(\Omega_f^{i}) \times
\mathbf{P}^l_h(\Gamma_f^{ij}) \times X^l_h(\Omega^{i}) \times
P^l_h(\Gamma^{ij}) \times H^2(\Lambda)$
 be  the velocity, the pressure, the stress vector,
the temperature, the heat flux and the beam displacement over the
corresponding subdomains. The   variable state $ (\vec{u}^{i,l},
p^{i,l},\vec{\tau}^{ij,l}, T^{i,l}, q^{ij,l},w) $ satisfies the
Navier -Stokes discrete  system
\begin{gather}
\begin{aligned}
&\langle \rho_f\,\frac{D\,\vec{u}^{i,l}(\vec{x}_h,t)}{D\,t},\vec{v}^{i,l}_{f}
\rangle
+ a_f\big(\vec{u}^{i,l},\vec{v}^{i,l}_{f}\big)
+ b_f\big(\vec{v}^{i,l}_{f},p^{i,l}\big)\\
&+ c_f\big(\vec{u}^{i,l}-\vec{u}_{g}^{i,l};
\vec{u}^{i,l},\vec{v}^{i,l}_{f}\big)
+\langle \vec{\tau}^{ij,l}, \vec{v}_{f}^{i,l}\rangle_{\Gamma_{f_h}^{ij}} \\
&=
-\langle \rho_f\,\vec{g}\,\beta(T^{i,l}-T_0),\vec{v}_{f}^{i,l}\rangle
\quad \forall \vec{v}_{f}^{i,l}\in \mathbf{X}_h^l(\Omega^i_{f} )\,,
 \label{SYS3A}
\end{aligned}\\
b_f\big(\vec{u}^{i,l},r_f^{i,l}\big) =  0  \quad
\forall r^{i,l}_f \in L^2 (\Omega^i_f )\,,
\\
\langle \vec{u}^{i,l} - \vec{U},\vec{s}^{i,l}_{f}\rangle_{\Gamma^{i}_{1_h}}=0
\quad \forall \vec{s}^{i,l}_{f}     \in \mathbf{P}_h^l(\Gamma_1^i)\,,
\\
\langle \vec{u}^{i,l}-  \vec{u}^{j,k},\vec{s}^{ij,l}_{f}
\rangle_{\Gamma^{ij}_{f_h}}   = 0
\quad \forall \vec{s}^{ij,l}_{f}  \in \mathbf{P}_h^l(\Gamma^{ij}_f)\,,
\label{BBBB}
\\
\langle \vec{u}^{i,l} - \dot{w}_h \cdot
\hat{n}_{sf},\vec{s}^{i,l}_{sf}\rangle_{\Gamma^{i}_{sf_h}}=0
\quad  \forall \vec{s}^{i,l}_{sf}  \in \mathbf{P}_h^l(\Gamma^{i}_{sf})\,,
 \label{EEEE}
\end{gather}
over $\Omega_f^i $ for all $j\in J_f^i$ and $i=1,2,\dots ,m$, the energy
equation
\begin{gather}
\begin{aligned}
&\langle \rho\, c_p\,\frac{D\, T^{i,l}(\vec{x}_h,t)}{D\, t},v^{i,l}\rangle
+ a\big(T^{i,l},v^{i,l}\big)\\
&+ c\big(\vec{u}^{i,l} -\vec{u}_{g}^{i,l};T^{i,l},v^{i,l}\big)
+\langle \vec{q}^{ij,l} , \vec{v}^{i,l}\rangle_{\Gamma^{ij}_h}= 0
\quad \forall v^{i,l} \in {X}^l_h(\Omega^i )\,,
\end{aligned}\label{SYS3B}
\\
\langle T^{i,l} - \Theta,s^{i,l}\rangle_{\Gamma_{2_h}^{i}}=0
\quad \forall s^{i,l}    \in P_h^l(\Gamma_2^i)\,,
\\
\langle T^{i,l}-  T^{j,k},\vec{s}^{ij,l} \rangle_{\Gamma^{ij}_h}  = 0
\quad \forall s^{ij,l}   \in P_h^l(\Gamma^{ij})\,,
\label{CCCC}
\end{gather}
over $\Omega^i $  for all $j\in J_f^i$ and $i=1,2,\dots ,m$, and
 the discrete Euler-Bernoulli beam equation
\begin{gather}
\begin{aligned}
&\langle \rho_s\, \delta\,\ddot{w}_h,v_{s_h}\rangle
+ a_s\big(w_h,v_{s_h}\big)\\
& = \langle p_h(\vec{x}_h(\xi,-\frac{\delta}{2}),t)-
 p_h(\vec{x}_h(\xi,\frac{\delta}{2}),t),v_{s_h}\rangle
 \quad v_{s_h} \in Z(\Lambda )\,, \end{aligned}
 \label{SYS3C} \\
w_h(0,t)=0\,,\quad \frac{\partial w_h(0,t)}{\partial \xi}=0\,,
\quad \dot{w}_h(0,t)=0\,. \label{DDDD}
\end{gather}
over $\Lambda $. On the shared boundaries $\Gamma^{ij}_{f_h}$ the
stress vectors, $\vec{\tau}^{ij, l}$ and $\vec{\tau}^{ji, k}$,
belong to the two different spaces $\mathbf{P}^{l}_h $ and
$\mathbf{P}^{k}_h $, and their equivalence should be considered
in the weak form by
\begin{eqnarray}\label{WEAKTAU}
\langle \vec{\tau}^{ij, l},\vec{v}^{i, l}\rangle_{\Gamma^{ij}_{f_h}}=
 -\langle \vec{\tau}^{ji, k},\vec{v}^{i, l}\rangle_{\Gamma^{ij}_{f_h}}\,.
\end{eqnarray}
Similarly $q^{ij, l}$ and $q^{ji, k}$ belong to the two
different spaces ${ P}^{l}_h $ and ${ P}^{k}_h $, and
\begin{equation} \label{WEAKHEAT}
\langle q^{ij, l},{v}^{i, l}\rangle_{\Gamma^{ij}_h}=
-\langle q^{ji, k},{v}^{i, l}\rangle_{\Gamma^{ij}_h}\,.
\end{equation}
Since different level meshes are generated by using midpoint
refinement the boundary vector spaces are nested, with ${\bf
R}^{l}_h \subseteq {\bf R}^{k}_h $ and ${R}^{l}_h \subseteq
{R}^{k}_h$, assuming $l \le k$. Then the traces of the test
functions $\vec{v}^{i,l}$ and ${v}^{i,l}$ on the coarsest boundary
can be decomposed as a linear combination of traces of test
functions $\vec{v}^{j,k}_a$ and ${v}^{j,k}_a$ on the finest
boundary. Under these hypotheses,  Eqs. (\ref{WEAKTAU}) and
(\ref{WEAKHEAT}) become
\begin{gather}\label{WEAKTAU2}
\langle \vec{\tau}^{ij, l},\vec{v}^{i, l}\rangle_{\Gamma^{ij}_{f}}=
-\langle \vec{\tau}^{ji, k},\sum_a w_a\,\vec{v}^{j, k}_a\rangle_{\Gamma^{
ji}_{f_h}}:=
-\langle \vec{\tau}^{ji, k},\mathcal{R}_k^{l}\vec{v}^{j, k}
\rangle_{\Gamma^{ji}_{f_h}}\,, \\
\langle q^{ij, l},{v}^{i, l}\rangle_{\Gamma^{ij}_h}= -<q^{ji, k},\sum_a w_a
\,{v}^{j, k}_a\rangle_{\Gamma^{ji}_h}:=
-\langle q^{ji, k},\mathcal{R}_k^{l}{v}^{j, k}\rangle_{\Gamma^{
ji}_h}\,,\label{WEAKHEAT2}
\end{gather}
which means that in  (\ref{SYS3A}) and (\ref{SYS3B}) the
Lagrange multipliers
$\langle \vec{\tau}^{ij, l},\vec{v}^{i, l}\rangle_{\Gamma^{ij}_{f}}$ and
$\langle q^{ij, l},{v}^{i, l}\rangle_{\Gamma^{ij}_h}$ can be always
discretized and computed on the finest grid available on the
boundary and then restricted over the coarsest one. Eqs.
(\ref{WEAKTAU2}) and (\ref{WEAKHEAT2}), define the restriction
operator $\mathcal{R}_k^{l}$ from the finest to the coarsest
vector space.

In  (\ref{BBBB}) and (\ref{CCCC}) the equalities between
$\vec{u}^{i,l}$ and  $\vec{u}^{j,k}$ on the boundary
$\Gamma^{ij}_{f_h}$, and between  $T^{i,l}$ and  $T^{j,k}$ on the
boundary $\Gamma^{ij}_h$, must be  considered in their weak form
since the  velocities and the temperatures belong to the two
corresponding different vector spaces ${\bf R}^{l}_h $, ${\bf
R}^{k}_h $ and ${R}^{l}_h $, ${R}^{k}_h$ respectively. Assuming
$l\le k$,  the spaces are nested with ${\bf R}^{l}_h \subseteq
{\bf R}^{k}_h $ and ${R}^{l}_h \subseteq {R}^{k}_h $. The weak
equations may turn into pointwise  equations if the velocities or
the temperatures belong to the coarsest spaces ${\bf R}^{l}_h $ or
${R}^{l}_h $, respectively.

\subsection{Solution strategy}

Although the  equation system (\ref{SYS3A})-(\ref{DDDD})
is fully coupled, its solution
is achieved with an iterative strategy, where the
three systems of equations are
solved separately and  in succession,
always using the latest information, until convergence
is reached. An iterative multigrid solver
is used for both the
Navier-Stokes and the energy equation systems since
the number of unknowns could be quite large. For the
solution of the beam equation a direct LU decomposition is used.

At each iteration, the linearized Navier-Stokes
system is assembled, using the
latest updated value of the temperature $T$ and
the latest updated value of the grid velocity $\vec{u}_{g}$
in the nonlinear term $c_f\big(\vec{u}-\vec{u}_{g};\vec{u},\vec{v}_f\big)$.
In the nonlinear term, the first of the two velocity $\vec{u}$
is  considered  explicitly.
On the boundary $\Gamma_{sf}$ Dirichlet boundary conditions
are imposed according to the latest updated value
of the beam displacement time derivative $\dot{w}$.
A V-cycle  multigrid algorithm is used to obtain a new update
solution for the pressure $p$ and the velocity $\vec{u}$.
Then the  energy equation system is assembled,
using  the previously evaluated velocity and grid velocity
in the advection term
$c\big(\vec{u}-\vec{u}_{g};T,{v}\big)$.
A multigrid V-cycle is solved and updated values of the
temperature $T$ are found.
Finally the beam equation system is built,
 where the load field is
computed using the previous evaluated pressure
$p$. Since the number of the subdomain unknowns
is limited an direct LU decomposition solver can be used for computing
the new displacement $w$ and its time derivatives.
The grid velocity is then computed according to Eqs. (\ref{LAEU1})-(\ref{LAEU3})
and the grid nodes are advected
along the corresponding characteristic lines.
The whole procedure is repeated until convergence
is finally reached.

The Navier-Stokes system (\ref{SYS3A})-(\ref{EEEE}) is solved using
a fully coupled iterative multigrid solver \cite{VA}
 with a Vanka-type smoother.
Multigrid solvers for coupled velocity/pressure system compute
simultaneously the solution for both the pressure and the velocity
field, and they are known to be one of the best class of solvers
for laminar Navier-Stokes equations (see for examples
\cite{ST,T1}). An iterative coupled solution for the linearized
discretized Navier-Stokes system requires the solution of a large
number of sparse saddle point problems. In order to optimally
solve the equation system (\ref{SYS3A})-(\ref{EEEE}), involving
the unknown stress vector $ \vec{\tau}^{ij}$, we use the block
Gauss-Seidel method, where each block consists of a small number
of degrees of freedom (for details see \cite{T1,VA,JO2,JT,PB}).
The characteristic feature of this type of smoother is that in
each smoothing step a large number of small linear systems of
equations has to be solved. Each block of equations corresponds to
all the degrees of freedom which are connected to few elements.
For example, for conforming finite elements, the block may consist
of all the elements, containing some pressure vertices. Thus, a
smoothing step consists of a loop over all the blocks, solving
only the equations involving the unknowns inside the elements that
are around the considered pressure vertices. The velocity and
pressure variables are updated many times in one smoothing step.

The Vanka smoother employed in our multigrid solver involves
the solution of a small number of degrees of freedom given by
the conforming Taylor-Hood finite element discretization used.
For this kind of element the pressure is computed only at the vertices
while the velocity field is computed also at the midpoints.
Over the internal part of the generic subregion $\Omega^{i}_{f_h}$,
where there are no boundary elements, our Vanka-block consists
of an element and all its neighboring elements.
We solve for all the degrees of freedom
inside the block, with boundary condition taken on the external boundaries.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.9\textwidth]{fig4}
\end{center}
\caption{Unknowns (black circles) and boundary conditions
(white circles) are shown for both
the velocity field ${\vec{ u}_h}$ (a) and the pressure $p_h$ (b) for our
particular Vanka-block smoother.} \label{fig4}
\end{figure}

For example, in Figure \ref{fig4}, our block consists of four vertex points
and 12 midpoints to be solved, for a total of 36 unknowns.
We have also used different blocks with different
performances but we have found this particular block
to be very robust and reliable even at high Reynolds numbers.
Examples of computations with this kind of solver can be found
in \cite{AMS1,AMS2,ST,T1}.

The fact that the solution is searched locally allows us to solve
for $ \vec{\tau}^{ij} $ only near the common boundary
$\Gamma^{ij}_{f_h}$, where different level meshes are coupled, and
to solve inside the subdomains by using Vanka smoother. When
solving block elements near the common boundary, distinction
should be made if the block element is inside the coarse or the
fine subregion. Let us consider the two adjacent subdomains
$\Omega^{i, l}_{f_h}$ and $\Omega^{j, k}_{f_h}$ with $l\le k$.
We have already seen, that in our particular discretization, where
the nodes of the coarsest mesh are included in the finest one, the
traces of the two velocities $u^{i, l}$ and $u^{j, k}$  are
equal and belong to the coarsest vector spaces $\mathbf{R}_h^l$.
Thus, when solving boundary blocks over the finest subregion
$\Omega^{j, k}_{f_h}$, we can use the standard solution technique
imposing the trace of the coarsest velocity $u^{i, l}$ as
Dirichlet boundary condition for the block equation system. We
have also seen that each Lagrange multiplier in the coarsest
system can be represented has a linear combination of Lagrange
multipliers of the finest system, thus when solving over boundary
blocks on the coarsest subregion $\Omega^{i, l}_{f_h}$, the
Lagrange multiplier term in  (\ref{SYS3A}) is replaced with
(\ref{WEAKTAU2}), where the stress tensor vector
$\vec{\tau}^{ji, k}$ is evaluated over the finest subdomain, by
using the pressure, $p^{j, k}$, and velocity field,
$\vec{v}^{j, k}$. In this way a constant and reciprocal exchange
of information between the two subdomains is ensured.

In order to increase the convergence rate, the considered
Vanka-type smoother has been coupled with a standard V-cycle
multigrid algorithm. The multigrid does not change the nature of
the solver, but allows the information to travel faster among
different parts of the domain. A rough global solution is
evaluated on the coarsest mesh $l=0$ and projected on the finer
grid $l=1$, where Vanka-loops are performed improving its details.
The updated solution is then projected on the mesh level $l=2$ and
improved. The procedure is repeated until the finest mesh is
reached. Solving the equation system in fine meshes improves
solution details, but at the same time reduces the communication
speed over the domain. However, this does not affect the global
convergence rate since a considerable information exchange among
different parts of the domain has been already done when solving
in coarser mesh levels. The analysis of the convergence of the
Vanka type multigrid solvers and the associated invertibility of
the discrete system for the Navier-Stokes can be found in
\cite{SM, AMS1}. All these considerations can be directly extended
to the energy equation  solver, where the same element block is
considered. For an interior block, as in Figure \ref{fig4}-a, the only
unknowns are the values of the temperature at the vertices and in
the midpoints, for a total of 12 variables, while for a boundary
block the heat flux $q^{ij,l}$ and $q^{ji,k}$ should be also
solved  over the common boundary $\Gamma^{ij}_h$.

\section{Numerical Experiments}

In this section we test the FSTI non-conforming multilevel
formulation and
its solution strategy.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig5}
\end{center}
\caption{Computational domain.} \label{fig5}
\end{figure}

As shown in Figure \ref{fig5}
let  the  rectangular region  $\Omega = [4\,m]\times [2\,m] $
be   the computational domain with boundary  $\Gamma$
and
$\Omega_f$ and $\Omega_s$ be the fluid  and the solid
subdomain respectively.
The solid
region $\Omega_s$ consists of
a  beam, clamped at the point $(1\,m,0)$,  with  length equal to $0.5$m
and thickness equal to $0.04$m.
 The fluid
and the solid boundaries, $\Gamma_f$ and $\Gamma_s$
are the contours
of the two shaded regions and their intersection
is labeled by $\Gamma_{sf}$.

On the left side of the domain inflow boundary
conditions are imposed for the velocity field $\vec{u}=(u_1,u_2)$
with parabolic profile $ u_1=0.1\,y\,(2-y)$m/s and  $ u_2=0 $.
On the right side of the domain outflow boundary
conditions are imposed
while on the remaining part of the boundary
non-slip conditions are set.
The temperature is set equal to zero in the inlet region
and to $100^o$C
on the solid boundary where the beam is clamped.
Adiabatic conditions
are imposed on the rest of the domain.
The initial conditions
for both the temperature and the velocity field are zero.

The fluid and the solid properties are chosen in
order to produce a
large deformation of the beam.
This choice implies strong interactions
among all the parts of the system and test the reliability
of the solver in challenging situations.
In the Navier-Stokes system, the fluid density $\rho_f$, the
viscosity   $\mu_f$, the
volumetric expansion coefficient $\beta$ and
the reference temperature $T_0$
 are equal to
$100$kg/m$^3$,  $0.01K$g/m\,s, $ 0.01\,K^{-1} $
and $0^o$C, respectively. In the temperature equation
the solid density $\rho_s$ is $200$kg/m$^3$,
 while the heat capacity
$c_p$ and the heat conductivity $K$,
are $100\,J/Kg\,K $ and $10\,W/m\,K $ in the fluid region,
and $10\,J/Kg\,K $
and $400\,W/m\,K$ in the solid region respectively.
The stiffness for unitary  length of the beam is
equal to 1kg\,m$^2$/s$^2$.

In  all the simulations the same time step
 $\Delta T=0.01\,s$ is used,
for a total of $500$ time steps ($5$ seconds).
Only the four level meshes,
$l_0$,$l_1$,$l_2$ and $l_3$,
are considered and
in Figure \ref{fig6}, the two different level meshes,
$l_0$ and $l_3$, are shown.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig6a}\quad
\includegraphics[width=0.45\textwidth]{fig6b}
\end{center}
\caption{The coarse mesh level $l_0$ (left) and finest level mesh
$l_3$ (right).} \label{fig6}
\end{figure}

The coarse mesh level $l_0$ has  $207$ elements, while the mesh
level $l_3$ obtained after three  consecutive midpoint refinements
has $13248$ elements. The one-dimensional mesh on the beam axis
follows the same midpoint refinement algorithm used for the
two-dimensional computational domain $\Omega$. On the coarse level
$l_0$ three elements are available, while on the fine grid $l_3$
after 3 refinements, the number of elements becomes $24$. Since
the number of unknowns is quite small, $(24+1)\times 2=50$, the
solution of the beam equation is always evaluated on the finest
mesh using a direct LU decomposition.

The results obtained with our coupled model (case B) are first
compared with the results obtained for the same geometry with a
rigid beam, $EI=\infty$, and with zero buoyancy force, $\beta=0$
(case A). All the computations are done at the time $ t=5\,s$ and
over the finest level $l_3$.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig7a} \quad
\includegraphics[width=0.45\textwidth]{fig7b}
\end{center}
\caption{Beam bending and grid distortion (left); velocity field map (right).}
\label{fig7}
\end{figure}

In Figure \ref{fig7} on the left, the beam bending and the
corresponding grid deformation  are displayed, showing the strong
influence of the pressure load on the beam shape. Figure
\ref{fig7} on the right shows the velocity field map and  clearly
indicates that the stationary solution is not reached since  new
vortices are constantly created and advected towards  the outflow
region.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig8a} \quad
\includegraphics[width=0.45\textwidth]{fig8b}\\
\includegraphics[width=0.45\textwidth]{fig8c}\quad
\includegraphics[width=0.45\textwidth]{fig8d}
\end{center}
\caption{Velocity field (top), pressure (bottom-left) and temperature
(bottom right) profiles evaluated along the section $y=0.5$
for $0\le x \le 4$ at the time $t=5$.} \label{fig8}
\end{figure}

In Figure \ref{fig8} the velocity field (top), the pressure
(bottom-left) and the temperature (bottom right) profiles,
evaluated over the section $y=0.5$ for $0\le x \le 4$, are shown
for the two cases $EI=\infty$, $\beta=0$ (case A) and  $EI=1$,
$\beta=0.01$ (case B). The combined effect of the beam deflection
and the buoyancy force modify considerably all the profiles,
pointing out how sensitive is the interaction among all the parts
of the system.

The number of  unknowns (velocity field, pressure, temperature and
displacement), involved in the computation at the mesh level $l_3$
is quite large, approximatively $94000$. However Figure \ref{fig7}
on the right shows that the only part of the system, subjected to
high vorticity is the region downstream of the beam. In the region
upstream of the beam and in the upper part of the domain, the
velocity field is almost stationary.

For this reason we remark that  more efficient computations can be
obtained if the solution is evaluated at the mesh levels $l_2$ or
$l_1$ in parts of the domain where the mesh level $l_3$ is not
needed. The domain $\Omega$ is splitted in  three  subdomains
$\Omega^1$, $\Omega^2$, $\Omega^3$, and three different
non-conforming meshes are built. The subdomains and the three
different non-conforming partitions are shown in Figure
\ref{fig9}.

\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.45\textwidth]{fig9a} \quad
\includegraphics[width=0.45\textwidth]{fig9b}\\
\includegraphics[width=0.45\textwidth]{fig9c} \quad
\includegraphics[width=0.45\textwidth]{fig9d}
\end{center}
\caption{Domain decomposition (top-left) and the different
configurations $P_1$ (top-right),
$P_2$ (bottom-left) and  $P_3$ (bottom right).} \label{fig9}
\end{figure}


\begin{figure}[ht]
\begin{center}
\includegraphics[width=0.7\textwidth]{fig10}
\end{center}
\caption{Oscillations of the beam extreme point for different
conforming and non-conforming meshes.} \label{fig10}
\end{figure}

In the subdomain $\Omega^3$ which is the solid domain $\Omega_s$,
the mesh level $l_1$ is always used. In the first configuration,
$P_1$ (top-right), the mesh levels, $l_2$ and $l_1$, are
considered for the subregions $\Omega^2$ and $\Omega^1$,
respectively. The different couplings of level meshes, $l_3$-$l_2$
and $l_3$-$l_1$ are used in the same subregions for the second
configuration $P_2$ (bottom-left), and the third configuration
$P_3$ (bottom-right). The numbers of nodes is greatly reduced for
all the three non-conforming configurations. In particular
approximatively $11000$, $39000$ and $28000$ are the new numbers
of unknowns for the new configurations $P_1$ $P_2$ and $P_3$,
respectively. In Figure \ref{fig10}, the oscillation of the beam
extreme point is compared for the three conforming meshes $l_0$,
$l_1$ and $l_3$, and for the 3 non-conforming meshes $P_1$, $P_2$
and $P_3$. The results show clear advantages of the non-conforming
discretizations over the conforming ones. Obviously the path
obtained with the finest mesh $l_3$ can be  considered the most
accurate. The $l_2$ path is very close to the $l_3$ in the first
second but differences appear as soon as the time increases. The
$l_1$ path is always below the $l_3$, showing too much stiffness
in the beam response. The beam oscillation obtained with the
non-conforming configuration $P_1$ perfectly overlaps the result
obtained with the conforming mesh $l_2$, and the result obtained
with the configuration $P_2$ perfectly overlaps the result in
$l_3$. It is possible to find very small differences between the
path in $l_3$ and the path in $P_3$, where there are two mesh
levels between the two adjacent regions $\Omega^1$ and $\Omega^2$.
These results clearly indicate how one can use the non-conforming
multilevel partitioning to preserve the same accuracy in regions
of interest, reducing at the same time the degree of freedom in
other parts of the domain. It should be noticed that the
comparisons have been tested on a beam displacement oscillation,
which is indirectly related to the multilevel domain partitioning.
The sensitivity of the beam's response to the different domain
decompositions points out again that the system is fully coupled.

\subsection*{Conclusions}
In this paper  a multigrid approach
applied to the domain decomposition for solving
Fluid-Structure-Thermal interaction problems has been presented
and tested. Our computational results clearly indicate that the
domain decomposition herein used in conjunction with the multigrid
method is a robust and reliable technique for solving FSTI
systems. The method developed in this paper leads to a fast and
flexible algorithm to compute solutions in different domains, over
different multigrid levels efficiently. Moreover in \cite{AMS1},
we had tested the performance of the computational methodology
presented in a parallel environment. The results are very
encouraging and the method has prompted us to investigate more
open and challenging questions involving coupled problems with
fluid-structure-thermal interaction that naturally arise in
engineering and science.

\begin{thebibliography}{00}

\bibitem{BMP1}{C.~Bernardi, Y.~Maday and A.~Patera};
{Domain decomposition by the mortar element method},
Asymptotic and numerical methods for partial differential
equation with critical parameters. H.K. et al. eds,
Reidel, Dordecht, pp. 269-286, 1993.

\bibitem{Belgacem}
F.~Ben Belgacem; The mortar finite element method with Lagrange Multipliers,
{Numer. Math.}, vol. 84(2), pp. 173-197, 1999.

\bibitem{SS2}
{P.~Seshaiyer and M.~Suri}; $hp$ submeshing via non-conforming finite
element methods,
{Comput. Meth. Appl. Mech. Eng.}, vol. 189, pp. 1011-1030, 2000.

\bibitem{BCP}
F.~Ben Belgacem, L.~K.~Chilton and P.~Seshaiyer;
The $hp$-Mortar Finite Element Method for Mixed elasticity and
Stokes Problems, {Comput. Math. Appl.}, vol. 46, pp. 35-55, 2003.

\bibitem{AMS1}{E.~Aulisa, S.~Manservisi and P.~Seshaiyer};
{A computational multilevel approach for solving 2D
Navier-Stokes equations over non-matching grids,}
to appear in {Comput. Meth. Appl. Mech. Eng.}, vol. 195, pp. 4604-4616, 2006.

\bibitem{AMS2}{E.~Aulisa, S.~Manservisi and P.~Seshaiyer};
{A non-conforming computational methodology for modeling coupled problems,}
{Nonlinear Analysis}, vol. 6, pp. 1445-1454, 2005.

\bibitem{BDW}{D.~Braess, W.~Dahmen  and C.~Wieners};
A Multigrid Algorithm for the Mortar Finite Element Method,
{SIAM J. Num. Anal.}, vol. 37(1), pp. 48-69, 1999.

\bibitem{GP}{J.~Gopalakrishnan and J.E.~Pasciak};
Multigrid for the Mortar Finite Element Method,
{SIAM J. Num. Anal.}, vol. 37(3), pp. 1029-1052, 2000.

\bibitem{REDDY}{J.N.~Reddy};
{An introduction to the finite element method,} $2^{nd}$ Edition,
 McGraw-Hill, New York, 1993.

\bibitem{GR1}{V.~Girault and P.~Raviart};
{The Finite Element Method for Navier-Stokes Equations:
Theory and Algorithms}, Springer, New York, 1986.

\bibitem{TM}{R.~Temam};
{Navier-Stokes equation}, North-Holland, Amsterdam, 1979.

\bibitem{SP}{E.W.~Swim and P.~Seshaiyer};
{A nonconforming finite element method for fluid-structure interaction
problems,}
{Comput. Meth. Appl. Mech. Eng.}, vol. 195(17-18), pp. 2088-2099, 2006.

\bibitem{DGH}
J.~Donea, S.~Giuliani, J.~Halleux;
An arbitrary Lagrangian Eulerian  finite element method for transient
fluid-structure interactions,
{Comput. Meth. Appl. Mech. Eng.}, vol. 33, pp. 689-723, 1982.

\bibitem{GGM}
C.~Grandmont, V.~Guimet, Y.~Maday;
Numerical analysis of some decoupling techniques for the approximation
of the unsteady  fluid structure interaction,
{Math. Models Methods Appl. Sci.}, vol. 11, pp. 1349-1377, 2001.

\bibitem{HLZ} T.~Hughes, W.~Liu, T.~Zimmermann;
Lagrangian Eulerian  finite element formulation for incompressible
viscous  flows,
{Comput. Methods Appl. Mech. Engrg.}, vol. 29, pp. 329-349, 1981.

\bibitem{BR}{J.H.~Bramble};
{Multigrid method}, Pitman Research Notes in Math, Longman,
London, vol. 294, 1993.

\bibitem{ST}{M.~Schafer and S.~Turek};
The benchmark problem: flow around a cylinder,
in Flow simulation with high performance
computers II. (E.H. Hirschel ed.),
{Notes on Numerical Fluid Mechanics}, vol. 52, pp. 547, 1996.

\bibitem{T1}{S.~Turek};
{Efficient solvers for incompressible flow problems:
an algorithmic and computational approach},
Lecture Notes in computational science and engineering,
Springer, vol. 6, 1999.

\bibitem{GH1}{M.~Gunzburger and L.H.~Hou}; Treating inhomogeneous
essential boundary conditions in finite element methods and the
calculation of boundary stresses,
{SIAM J. Num. Anal.}, vol. 29(2), pp. 390-424, 1992.

\bibitem{GM2}{M.~Gunzburger, L.H.~Hou and T.~Svobodny};
Analysis and finite element approximations of optimal
 control problems for the stationary   Navier-Stokes
equations with Dirichlet control,
{Math. Model. Numer. Anal.}, vol. 25, pp. 711-748, 1991.

\bibitem{VA}{S.~Vanka};
Block-implicit multigrid calculation of two-dimensional recirculation
flows, {Comput. Meth Appl. Mech. Eng.},
vol. 59(1), pp. 29-48, 1986.

\bibitem{JO2} {V.~John};
A Comparison of Parallel Solvers for
 the Incompressible Navier-Stokes Equations,
{Comput. Visual. Sci.}, vol. 1(4), pp. 193-200, 1999.

\bibitem{JT} {V.~John and L.~Tobiska};
Numerical Performance of Smoothers in
Coupled Multigrid Methods for the Parallel
Solution of the Incompressible Navier-Stokes Equations,
{Int. J. Num. Meth. Fluids}, vol. 33, pp. 453-473, 2000.

\bibitem{PB}{M.F.~Paisley and N.M.~Bhatti};
Comparison of multigrid methods for neutral
and stably stratified flows over two-dimensional obstacles,
{J. Comput. Phys.}, vol. 142, pp. 581-610, 1998.

\bibitem{SM}{S.~Manservisi};
Numerical Analysis of Vanka-Type Solvers for Steady Stokes and
Navier-Stokes Stokes Flows, {SIAM J. Numer. Anal.}, vol. 44,
pp. 2025-2056, 2006.

\end{thebibliography}

\end{document}