Numerical Methods For Unsteady Thermal Fluid Structure Interaction

XXXX, 1–39 © De Gruyter YYYY
Numerical Methods for Unsteady Thermal

Fluid Structure Interaction
Philipp Birken and Azahar Monge
Abstract. We discuss thermal fluid structure interaction processes, where a simulation of the
time dependent temperature field is of interest. Thereby, we consider partitioned coupling
schemes with a Dirichlet-Neumann method. We present an analysis of the method on a model
problem of discretized coupled linear heat equations. This shows that for large quotients in
the heat conductivities, the convergence rate will be very small. The time dependency makes
the use of time adaptive implicit methods imperative. This gives rise to the question, how
accurate the appearing nonlinear systems should be solved, which is discussed in detail for
both the nonlinear and linear case. The efficiency of the resulting method is demonstrated
using realistic test cases.
Keywords. Thermal Fluid Structure Interaction, Conjugate Heat Transfer, Dirichlet-Neumann

method, Time adaptivity, Termination criteria.
AMS classification. 65F10, 65L04, 65M22, 74F04, 74F10.
1 Introduction
Fluid structure interaction occurs when a deformable or moving structure interacts
with a surrounding or internal fluid flow. Our specific field of interest is thermal in-
teraction between fluids and structures, also called conjugate heat transfer. Examples
for thermal fluid structure interaction are cooling of gas-turbine blades, thermal anti-
icing systems of airplanes [12], supersonic reentry of vehicles from space [31, 25], gas
quenching, which is an industrial heat treatment of metal workpieces [23, 40] or the
cooling of rocket nozzles [27, 26]. These problems are usually too complex to solve
them analytically, and therefore, numerical simulations of conjugate heat transfer are
essential in many applications.
The efficient numerical simulation of fluid structure interaction (FSI) models is one
of the important current challenges in scientific computing as stated in [11]:
“The issue of coupling models of different events at different scales and
governed by different physical laws is largely wide open and represents an
enormously challenging area for future research.”
In this article, we focus on the coupling between air and an alloy. When being
cooled or heated, the alloy experiences thermomechanical effects that change its in-
ternal structure. A first example are steel forging processes. One possibility is to use
2 P. Birken and A. Monge
(a) Inductive heating. (b) Thermo-mechanical (c) Local air-cooling.

forming.
Figure 1. Gas quenching. Left and center picture: Institute of Mechanics and Dynamics,
University of Kassel. Right picture: Institute of Metal Forming Technology, University
of Kassel.
cold high pressured air as a cooling medium (see Figure 1). Knowledge about the
time dependent temperature field is imperative to predict where martensite will be in a
finished steel part. This allows to predict material properties generated by the forging
process.
Another example is the cooling of rocket thrust chambers. The temperature achieved
in the turbines has increased over the years due to progress in the building materials
and therefore, advanced cooling methods are needed. This is both to avoid critical
damage to the rocket nozzle when in use and to develop reusable rocket stages. An
important case is the Ariane 5, see Figure 2a, which is used to deliver payloads into the
geostationary transfer orbit (GTO) or the low Earth orbit (LEO). The first stage rocket
engine for the Ariane 5 is the Vulcain 2, see figure 2b.
The first stage rocket engine is only used during the launch of the rocket. Therefore,
recovering and reusing the first stage will reduce the cost of space access and the
environment impact. Related to this, the company SpaceX is developing a set of new
technologies for an orbital launch system that may be reused many times in a manner
similar to the reusability of aircraft. The first controlled vertical splashdown of an
orbital rocket stage on the ocean surface was achieved in April 2014. The next two
flights in January and April 2015 attempted to land the returning first stage on a floating
platform. Both of them were guided accurately to the target, but they did not succeed
in landing vertically on the floating platform and were destroyed [3]. Finally, the first
vertical landing was achieved on December 21, 2015, when the first stage of Falcon 9
Flight 20 successfully landed on solid ground [1]. On April 8, 2016, Flight 23 achieved
the first soft landing on a drone ship in the Atlantic Ocean [2], see Figure 3.
In view of these achievements, it is of interest to simulate several cycles of the
combustion-cooling process of the thrust chamber. This will tell how many times the
engine can be used without damaging the structure. More details about this can be
found in [27, 26].
Numerical Methods for Unsteady Thermal Fluid Structure Interaction 3
(a) Ariane 5 on the launch (b) The Vulcain engine (c) Sketch of the rocket
pad. in a museum. thurst chamber.
Figure 2. Cooling of rocket thrust chambers. Left picture: DLR German Aerospace
Center, CC BY 2.0. Center picture: Pline, CC BY-SA 3.0
(a) Unsuccessul vertical (b) First stage landing ver- (c) First stage landed on au-
landing attempt. tically on solid ground in tonomous droneship in April 2016.
December 2015.
Figure 3. SpaceX reusable launch system development program. SpaceX Photos, CC0
1.0.
Figure 2c shows a combustion chamber with the nozzle, through which hot gas can
escape. The nozzle is delimited by a structure that can be damaged due to the high
temperature of the gas flowing inside. In order to avoid this, a cooling fluid flows
through small channels contained inside the structure. This results in a system with
two thermal interactions between fluids and structures. On one hand, between the hot
gas coming out from the combustion chamber and the structure recovering the nozzle.
On the other hand, between the cooling fluid and the structure. A sketch of the coupling
surfaces can be consulted in Figure 4 where Γs,cf corresponds to the interface between
the structure and the cooling fluid and Γs,hg to the interface between the structure and
the hot gas.
Figure 4. Sketch of the coupling surfaces. Figure taken from [27].
With regards to the space discretization, the use of the finite element method (FEM)
is ubiquitious for the structure. For the fluid, typically finite volume methods (FVM)
are used. However, there are some approaches to use finite elements for both problems,
in particular when using a monolithic method [4].
1.1 Partitioned Coupling Methods

To simulate FSI problems there exist two main methods: monolithic and partitioned
ones. In the monolithic method, a new code is tailored for the coupled equations,
whereas the partitioned approach allows to reuse existing software for each sub-problem.
The coupling is done by a master program which calls interface functions of the other
codes [16]. If the data transfer between the subsolvers is done only once per time step,
we are using a loosely coupled scheme [17]. However, for stability reasons, often a
strongly coupled scheme needs to be used [28]. In this case the data exchange at ev-
ery time step is repeated until a convergence criterion is satisfied. Here, we focus on
partitioned methods.
At the boundary, one imposes that the temperature and the heat flux have to be con-
tinuous accross the interface. To obtain such a solution, a Dirichlet-Neumann iteration
can be employed. This consists of solving the fluid problem with Dirichlet boundary
conditions at the interface and then the structure problem with Neumann boundary
conditions at the interface, resulting in a fixed point iteration. Figure 5 illustrates this.
This iteration is a basic method in both domain decomposition and fluid structure in-
teraction.
In the domain decomposition context, it has two main problems, namely slow con-
vergence and the need for an implementation using a red-black colouring. The slow
Figure 5. Illustration of the FSI solver.
convergence can be improved using a relaxation procedure [34]. In fluid structure

interaction, there are only two domains, coupled along an interface, making the appli-
cation straight forward. However, the convergence rate is not great for the coupling
between a compressible fluid and a structure [14], which is why a lot of effort goes
into convergence acceleration. On the other hand, the Dirichlet-Neumann iteration
was reported to be a very fast solver for thermal fluid structure interaction [8].
As mentioned before, we study transient processes and therefore a numerical method
needs to be chosen for the time discretization. In [29], the implicit midpoint rule is used
in a monolithic scheme to analyze energy conservation of an aeroelasticity problem.
Already in [5, 13], it is suggested to use an explicit high order Runge-Kutta scheme
for both subproblems with data exchange at each stage. However, the resulting scheme
has limited time steps due to the explicit nature of the method. The order of coupling
schemes on moving meshes is analyzed in [22], but only first order convergence is
proved for p-th order schemes. Moreover, higher order implicit Runge-Kutta schemes
on moving meshes are analyzed in [43] (in 1D) and in [44] (in 3D). There, so called
explicit first stage, singly diagonally implicit Runge-Kutta schemes (ESDIRK) are em-
ployed and higher order in time is proved by numerical results. The master program
of the FSI procedure can be extended to SDIRK methods and furthermore, time adap-
tivity can be added into this framework as explained in [9, 10]. Another interesting
alternative is subcycling, where in the subsolvers, different time step sizes are chosen
[21, 37]. So far, this has not been brought together with time adaptive methods.
Summarizing, a partitioned method for unsteady thermal FSI involves two sub-
solvers (fluid and structure) and two embedded iterations. There exists an outer loop
corresponding to the time integration. Then, at each timestep, a subiteration (the
Dirichlet-Neumann method) couples the subsolvers, which have another nonlinear it-
eration so provide solutions to the subproblems.
An outline of the article now follows. In section 2, we describe the model and
discretization, as well as the coupling conditions and the partitioned solution method,
the Dirichlet-Neumann iteration. A model problem, namely two coupled discretized
heat equations, is presented in section 3. Our analysis to determine the convergence
rate of the Dirichlet-Neumann method for the model problem can be found in section
4. Numerical results are included to illustrate the theoretical analysis. In section 5,

we describe how a fast solver is obtained using time adaptivity, followed by an in
depth analysis of how tolerances in the various nested iterations should be chosen to
avoid oversolving while giving the desired accuracy in section 6. Numerical results to
demonstrate the efficiency of the obtained method are presented in section 7.
2 Treating Thermal FSI Problems

The basic setting we are in is that on a domain Ω1 ⊂ Rd where d corresponds to the
spatial dimension, the physics is described by a fluid model, whereas on a domain
Ω2 ⊂ Rd , a different model describing the structure is used. The two domains are
almost disjoint in that they are connected via an interface. The part of the interface
where the fluid and the structure are supposed to interact is called the coupling inter-
face Γ ⊂ ∂Ω1 ∪ ∂Ω2 . Note that Γ might be a true subset of the intersection, because
the structure could be insulated. At the interface Γ, coupling conditions are prescribed
that model the interaction between fluid and structure. For the thermal coupling prob-
lem, these conditions are that temperature and the normal component of the heat flux
are continuous across the interface.
2.1 Fluid Model

We model the fluid using the time dependent Reynolds Averaged Navier-Stokes equa-
tions (URANS), which are a second order system of conservation laws (mass, momen-
tum, energy) modeling turbulent compressible flow. We consider the two dimensional
case, written in conservative variables density ρ, momentum m = ρv and energy per
unit volume ρE, where a ˜ denotes the Favre average and the overbar the ensemble
average [6]:
∂t ρ + ∇ · ρṽ = 0,
2 2
X 1 X R

∂t ρṽ + ∂xj (ρṽi ṽj ) = −∂xj pδij + ∂xj S̃ij + Sij , i = 1, 2 (2.1)
Re
j=1 j=1
2
X 1 R 00 ] 00 00 Wj
∂t ρẼ + ∇ · (ρH̃ ṽj ) = ∂xj ( Sij − Sij )vi − ρvj + Sij vi − ρvj k +
f g .
Re ReP r
j=1
The Reynolds stresses

R 00 00
Sij = −ρv]
i vj
and the turbulent energy
d
1X 0 0
k= vj vj
2
j=1
are modelled using the Spallart-Allmaras model [39]. Furthermore, qf = (q1 , q2 )T

represents the heat flux and S = (Sij )i,j=1,2 the viscous shear stress tensor. As the
equations are dimensionless, the Reynolds number Re and the Prandtl number P r
appear. The system is closed by the equation of state for the pressure p = (γ − 1)ρe,
the Sutherland law representing the correlation between temperature and viscosity,
as well as the Stokes hypothesis. Additionally, we prescribe appropriate boundary
conditions at the boundary of Ω1 except for Γ, where we have the coupling conditions.
In the Dirichlet-Neumann coupling, a temperature value is enforced at Γ.
2.2 Structure Model

Regarding the structure model, we will consider heat conduction only. Thus, we have
the nonlinear heat equation for the structure temperature Θ
d
ρ(x)cp (Θ) Θ(x, t) = −∇ · q(x, t), (2.2)
dt
where
qs (x, t) = −λ(Θ)∇Θ(x, t)
denotes the heat flux vector. For alloys, the specific heat capacity cp and heat conduc-
tivity λ are temperature-dependent and highly nonlinear. How to model thermome-
chanical effects is subject of a lot of current research with important questions being
how to relate the microstructure to macroscopical models.
As an example, an empirical model for the steel 51CrV4 was suggested in [35]. This
was obtained simply by doing measurements and then a least squares fit to a chosen
curve. The coefficient functions are then
λ(Θ) = 40.1 + 0.05Θ − 0.0001Θ2 + 4.9 · 10−8 Θ3 (2.3)
and
!
e−cp1 (Θ)/10 + e−cp2 (Θ)/10
cp (Θ) = −10 ln (2.4)
2
with
cp1 (Θ) = 34.2e0.0026Θ + 421.15 (2.5)
and
cp2 (Θ) = 956.5e−0.012(Θ−900) + 0.45Θ. (2.6)
For the mass density one has ρ = 7836 kg/m3 .

Finally, on the boundary, we have Neumann conditions qs (x, t) · n(x) = qb (x, t).
2.3 Coupling Conditions

As mentioned before at the beginning of this section, the coupling conditions are that
temperature and the normal component of the heat flux are continuous across the in-
terface, i.e;
T (x, t) = Θ(x, t), x ∈ Γ, (2.7)
where T is the fluid temperature and Θ the structure temperature and
qf (x, t) · n(x) = qs (x, t) · n(x), x ∈ Γ. (2.8)
2.4 Discretization in Space

Following the partitioned coupling approach, we discretize the two models separately
in space. For the fluid, we use a finite volume method, leading to
d
u + h(u, ΘΓ ) = 0, (2.9)
dt
where h(u, ΘΓ ) represents the spatial discretization and its dependence on the tem-
peratures on the discrete interface to the structure, here denoted by ΘΓ .
Regarding structural mechanics, the use of finite element methods is ubiquitious.
Therefore, we will also follow that approach here, leading to the nonlinear equation
for all unknowns on Ω2 :
d
M(Θ) Θ + A(Θ)Θ = qfb + qΓb (u). (2.10)
dt
Here, M is the mass matrix, also called heat capacity matrix for this problem and
A is the heat conductivity matrix. The vector Θ consists of all discrete temperature
unknowns and qΓb (u) is the discrete heat flux vector on the coupling interface to the
fluid, whereas qfb corresponds to boundary heat fluxes independent of the fluid, for
example at insulated boundaries.
2.5 Time Discretization

For sake of completeness, we now write down a discretization using the implicit Euler
method. More advanced time integration schemes are discussed in Section 5. For the
coupled system (2.9)-(2.10) we obtain
un+1 − un + ∆tn h(un+1 , Θn+1

Γ ) = 0, (2.11)
M(Θn+1 )(Θn+1 − Θn ) + ∆tn A(Θn+1 )Θn+1 = ∆tn (qfb + qΓb (uk+1 )). (2.12)
2.6 The Dirichlet-Neumann Method

The Dirichlet-Neumann method is a basic iterative substructuring method in domain
decomposition and it is a common choice for treating FSI problems. Therefore, we
now employ it to solve the system (2.11)-(2.12). This corresponds to alternately solv-
ing equation (2.11) on Ω1 with Dirichlet data on Γ and (2.12) on Ω2 with Neumann
data on Γ.
Thus, one gets for the k-th iteration the two decoupled equation systems
uk+1 − un + ∆tn h(uk+1 , ΘkΓ ) = 0, (2.13)
M(Θk+1 )(Θk+1 − Θn ) + ∆tn A(Θk+1 )Θk+1 = ∆tn (qfb + qΓb (uk+1 )), (2.14)
with some initial condition Θ0Γ . The iteration is terminated according to the standard
criterion
kΘk+1
Γ − ΘkΓ k ≤ τ (2.15)
where τ is a user defined tolerance.

The convergence rate of the Dirichlet-Neumann iteration is not great for the cou-
pling between a compressible fluid and a structure [14], which is why a lot of effort
goes into convergence acceleration. On the other hand, the Dirichlet-Neumann itera-
tion was reported to be very fast solver for thermal fluid structure interaction. More
specifically, in [8] the iteration is extremely efficient and achieves a very accurate so-
lution with at most two iterations per timestep. To analyze the convergence behaviour
of the Dirichlet-Neumann iteration will explain why this coupling method is very effi-
cient in some FSI models and very inefficient in some others.
In principle, the convergence rate of the Dirichlet-Neumann method is analyzed
in any standard book on domain decomposition methods, e.g. [34, 41]. There, the
iteration matrix of the discretized equations is derived with respect to the interface
unknowns and the convergence rate is the spectral radius of that. However, due to the
nonlinearity of the thermal fluid structure interaction model explained previously, it is
not possible to compute the spectral radius of the iteration matrix in this case. For this
reason, we perform in the next section a convergence analysis of the coupling of two
linear heat equations. We chose this model because it is a basic building block in fluid
structure interaction.
3 A Model Problem: Coupled Heat Equations

In FSI models, the solid is typically discretized using finite elements. On the other
hand, although the finite element method (FEM) is applicable to computational fluid
dynamics, the finite volume method (FVM) is generally a better choice for the dis-
cretization of the fluid. This method guarantees the conservation of fluxes through a
particular control volume. Therefore, we present here a convergence analysis of the
unsteady transmission problem with mixed discretizations (FVM-FEM).
For this model problem, Henshaw and Chand provided in [24] a method to analyze
stability and convergence speed of the Dirichlet-Neumann iteration in 2D based on ap-
plying the continuous Fourier transform to the semi-discretized equations. Their result
depends on ratios of thermal conductivities and diffusivities of the materials. How-
ever, in the fully discrete case we have observed that the iteration behaves differently
in some cases. Therefore, we propose a complementary analysis for the fully discrete
case here.
3.1 Model Problem

The unsteady transmission problem is as follows, where we consider a domain Ω ⊂ Rd
which is cut into two subdomains Ω1 ∪ Ω2 = Ω with transmission conditions at the
interface Γ = Ω1 ∩ Ω2 :
∂um (x, t)
αm − ∇ · (λm ∇um (x, t)) = 0, t ∈ [t0 , tf ], x ∈ Ωm ⊂ R2 , m = 1, 2,
∂t
um (x, t) = 0, t ∈ [t0 , tf ], x ∈ ∂Ωm \Γ,
u1 (x, t) = u2 (x, t), x ∈ Γ,
∂u2 (x, t) ∂u1 (x, t)
λ2 = −λ1 , x ∈ Γ,
∂n2 ∂n1
um (x, 0) = u0m (x), x ∈ Ωm ,
(3.1)
where nm is the outward normal to Ωm for m = 1, 2.

The constants λ1 and λ2 describe the thermal conductivities of the materials on Ω1
and Ω2 respectively. D1 and D2 represent the thermal diffusivities of the materials and
they are defined by
λm
Dm = , with αm = ρm cpm (3.2)
αm
where ρm represents the density and cpm the specific heat capacity of the material
placed in Ωm , m = 1, 2.
We discretize this problem with a constant mesh width with respect to both spatial
components (∆y := ∆x = 1/(N + 1)) resulting in N 2 interior space discretization
points in both Ω1 and Ω2 . We use the implicit Euler method for the time discretization.
3.2 Semidiscrete Analysis

Before we present in the next section an analysis for the fully discrete equations,
we want to describe previous results which analyze the behaviour of the Dirichlet-
Neumann iteration for the transmission problem in the semi discrete case.
On one hand, a one dimensional stability analysis was presented by Giles [20].
There, while using the implicit Euler method in the subsolvers, an explicit time inte-
gration method was chosen with respect to the interface unknowns.
On the other hand, Henshaw and Chand provided in [24] a method to analyze con-
vergence speed of the Dirichlet-Neumann iteration. There, one applies the implicit
Euler method for the time discretization on both equations in (3.1) but keeps the space
continuous. Then, they applied the Fourier transform in space in order to transform the
second order derivatives into algebraic expressions. This converts the partial differen-
tial equations into a system of purely algebraic equations. Once we have a coupled
system of algebraic equations, we can insert one into the other and obtain the conver-
gence rate, called amplification factor β in [24]. They then derive the formula

1
r
λ1 D2 tanh − √
D2 ∆t
β = −
. (3.3)
λ2 D1 tanh √ 1
D1 ∆t
√ √ √
For ∆t big enough, we have tanh −1/ D2 ∆t ≈ −1/ D2 ∆t and tanh 1/ D1 ∆t ≈
√
1/ D1 ∆t and therefore:
r √
λ1 D2 D1 ∆t λ1
β≈ √ = . (3.4)
λ2 D1 D2 ∆t λ2
One observes in (3.4) that the rates of the iteration behave as the quotient of thermal
conductivities when ∆t → 0. This suggests that strong jumps in the thermal conduc-
tivities of the materials cause fast convergence.
3.3 Space Discretization

We now describe a rather general space discretization of the model problem. The core
property that we need is that the meshes of Ω1 and Ω2 are compatible on Γ (they share
the same nodes on Γ) as shown in Figure 6. Furthermore, we need that there is a
specific set of unknowns associated with the interface nodes.
Figure 6. Splitting of Ω between finite volumes and finite elements.
(1)
Then, letting uI correspond to the unknowns on Ω1 and uΓ to the unknowns at
the interface Γ, we can write a general discretization of the first equation in (3.1) in a
compact form as:
(1) (1) (1) (1)

M1 u̇I + MIΓ u̇Γ + A1 uI + AIΓ uΓ = 0. (3.5)
On the other hand, a general discretization of the first equation in (3.1) on Ω2 can
be written as:
(2) (2) (2) (2)

M2 u̇I + MIΓ u̇Γ + A2 uI + AIΓ uΓ = 0. (3.6)
(2)
where uI correspond to the unknowns on Ω2 .
However, the system (3.5)-(3.6) is not enough to describe (3.1). Additionally, we
need an approximation of the normal derivatives at Γ. If an FVM is used over Ω1 ,
we approximate the normal derivative with respect to u1 using second order one-sided
finite differences:
∂u1 λ1
−λ1 ≈ (4u1,N (t) − u1,N −1 (t) − 3uΓ ). (3.7)
∂n1 2∆x
On the other hand, if an FEM is used over Ω1 and φj is a nodal basis function for a
node on Γ we observe that the normal derivative with respect to u2 can be written as
linear functionals using Green’s formula [41, pp. 3]. Thus, the approximation of the
normal derivative is given by
Z Z
∂u2
λ2 φj dS = λ2 (∆u2 φj + ∇u2 ∇φj )dx
Γ ∂n2 Ω2
Z Z (3.8)
d
= α2 u2 φj + λ2 ∇u2 ∇φj dx.
Ω2 dt Ω2
Consequently, the equation
(2) (2) (2) (2) (2) (2) (1) (1) (1) (1) (1) (1)
MΓΓ u̇Γ + MΓI u̇I + AΓΓ uΓ + AΓI uI = −MΓΓ u̇Γ − MΓI u̇I − AΓΓ uΓ − AΓI uI ,
(3.9)
is a discrete version of the fourth equation in (3.1) and completes the system (3.5)-
(3.6). We now reformulate the coupled equations (3.5), (3.6) and (3.9) into an ODE
T
(1) (2)
for the vector of unknowns u = uI , uI , uΓ
M̃u̇ + Ãu = 0, (3.10)
where
 (1)   (1) 
M1 0 MIΓ A1 0 AIΓ
(2) (2)
M̃ =  0 M2 MIΓ  , Ã =  0 A2 AIΓ .
   
(1) (2) (1) (2) (1) (2) (1) (2)
MΓI MΓI MΓΓ + MΓΓ AΓI AΓI AΓΓ + AΓΓ
3.4 Time Discretization

Applying the implicit Euler method with time step ∆t to the system (3.9), we get for
(1),n+1 (2),n+1 n+1 T
the vector of unknowns un+1 = (uI , uI , uΓ )
Aun+1 = M̃un , (3.11)
where
 (1) (1) 
M1 + ∆tA1 0 MIΓ + ∆tAIΓ
(2) (2) 
A = M̃ + ∆tÃ =  0 M2 + ∆tA2 MIΓ + ∆tAIΓ  ,

(1) (1) (2) (2)
MΓI + ∆tAΓI MΓI + ∆tAΓI MΓΓ + ∆tAΓΓ
(1) (2) (1) (2)

with MΓΓ = MΓΓ + MΓΓ and AΓΓ = AΓΓ + AΓΓ .
3.5 Fixed Point Iteration

We now employ a Dirichlet-Neumann iteration to solve the discrete system (3.11).
This corresponds to alternately solving the discretized equations of the transmission
problem (3.1) on Ω1 with Dirichlet data on Γ and the discretization of (3.1) on Ω2 with
Neumann data on Γ.
Therefore, from (3.11) one gets for the k-th iteration the two equation systems
(1),n+1,k+1 (1) (1) (1),n (1)

(M1 + ∆tA1 )uI = −(MIΓ + ∆tAIΓ )un+1,k
Γ + M1 uI + MIΓ unΓ , (3.12)
Âûk+1 = M̂un − bk , (3.13)

to be solved in succession. Here,
(2) (2)
! (2)
!
M2 + ∆tA2 MIΓ + ∆tAIΓ 0 M2 MIΓ
Â = (2) (2) (2) (2) , M̂ = (1) (2) ,
MΓI + ∆tAΓI MΓΓ + ∆tAΓΓ MΓI MΓI MΓΓ
and
!
k 0
b = (1) (1) (1),n+1,k+1 (1) (1) , (3.14)
(MΓI + ∆tAΓI )uI + (MΓΓ + ∆tAΓΓ )un+1,k
Γ
!
(2),n+1,k+1
k+1 uI
û = ,
uΓn+1,k+1
with some initial condition, here un+1,0

Γ = unΓ . The iteration is terminated according
k+1 k
to the standard criterion kuΓ − uΓ k ≤ τ where τ is a user defined tolerance [7].
One way to analyze this method is to write it as a splitting method for (3.11) and
try to estimate the spectral radius of that iteration. However, the results obtained in
this way are much too inaccurate. For that reason, we now rewrite (3.12)-(3.13) as
an iteration for un+1
Γ to restrict the size of the space to the dimension of uΓ . To this
(1),n+1,k+1 (2),n+1,k+1
end, we isolate the term uI in (3.12) and uI in the first equation
in (3.13) and we insert the resulting expressions into the second equation in (3.13).
Consequently, the iteration un+1,k+1
Γ = Σun+1,k
Γ + ψ n is obtained with iteration matrix
−1 (1)
Σ = −S(2) S , (3.15)
where
(m) (m) (m) (m) (m) (m)

S(m) = (MΓΓ + ∆tAΓΓ ) − (MΓI + ∆tAΓI )(Mm + ∆tAm )−1 (MIΓ + ∆tAIΓ ),
(3.16)
for m = 1, 2 and ψ n contains terms that depend only on the solutions at the previous
time step. Notice that Σ is a discrete version of the Steklov-Poincaré operator.
Thus, the Dirichlet-Neumann iteration is a linear iteration and the rate of conver-
gence is described by the spectral radius of the iteration matrix Σ.
4 Convergence Analysis
The derivation so far was for a rather general discretization. It is now possible to look
at specific discretizations. In this section, we study the iteration matrix Σ for a specific
FVM-FEM discretization.
The subdomains are here Ω1 = [−1, 0] × [0, 1], Ω2 = [0, 1] × [0, 1]. An equidis-
tant grid is chosen i.e, ∆x = ∆y = 1/(N + 1). At each discrete point of the fi-
nite volume discretization xi,j , we integrate over the cell Ii,j = [xi−1/2,j , xi+1/2,j ] ×
[xi,j−1/2 , xi,j+1/2 ] and use the flux function
λ1
F (uL , uR ) = − (uR − uL ), (4.1)
∆x
to approximate the flux, which results in a second order scheme. For the FEM dis-
cretization, we use triangular elements distributed as sketched in Figure 7 and the
following pyramidal test functions
 x+y
∆x − 1, if x = (x, y) ∈ Region 1,


y

∆x , if x ∈ Region 2,





∆x−x




φk (x, y) = 1 − x+y
∆x , if x ∈ Region 4, (4.2)
 ∆x−y




 x





0, otherwise.

5
5
6 4
6 4
3
1 2 5 1 3
6 4 2
5 1 3
6 4 26 5 4
1 3 1 3
2 2
Figure 7. Sketch of the regions for the pyramidal test functions defined in (4.2).
The discretization matrices are given in this case by

   
−B I 0 B −I 0
 .   .. 
λ1  I −B . .  λ2  −I B . 
A1 =  , A2 = ,
  
∆x2
 .. .. 2
∆x 
 .. ..

 . . I 
  . . −I 

0 I −B 0 −I B
where
 
4 −1 0
 .. 
 −1 4 .
B= ,
 
.. ..

 . . −1 

0 −1 4
and I ∈ RN ×N is an identity matrix. Note that each block of the matrices Am ∈

2 2
RN ×N has size N × N .
 
N Ñ 0
 T .. 
 Ñ N . 
M2 = α2  ,
 
.. ..

 . . Ñ 

T
0 Ñ N
where
   
5/6 −1/12 0 −1/12 1/4 0
 ..   . 
 −1/12 5/6 .   0 −1/12 . . 
N= , Ñ = .
   
.. ..   .. ..

 . . −1/12 


 . . 1/4 

0 −1/12 5/6 0 0 −1/12
2 2
N ×N has size N × N as well. We consider here
block of the matrix M2 ∈ RT
Each
2
Ej = 0 · · · 0 I 0 · · · 0 ∈ RN ×N where the only nonzero block is the
j-th block of size N × N . Thus,
(1) λ1 (2) λ2 (2)

AIΓ = 2
EN , AIΓ = − 2 E1 , MIΓ = α2 E1 Ñ,
∆x ∆x
 
5/12 −1/24 0
 .. 
(2)
 −1/24 5/12 . 
MΓΓ = α2  ,
 
.. ..

 . . −1/24 

0 −1/24 5/12
 
2 −1/2 0
 .. 
(1) 3λ1 (2) λ2  −1/2 2 . 
AΓΓ = I, AΓΓ = ,
 
2∆x2 ∆x2
 .. ..

 . . −1/2 

0 −1/2 2
(2) (m)
where MΓΓ and AΓΓ ∈ RN ×N for m = 1, 2.
(1) λ1 (2) λ2 (2)

AΓI = (4ETN − ETN −1 ), AΓI = − 2 ET1 , MΓI = α2 ET1 Ñ.
2∆x2 ∆x
(1) (1) (1)
In this specific case, M1 = α1 I, MIΓ = MΓΓ = MΓI = 0. In particular,
(1) (1) (1)

S(1) = ∆tAΓΓ − ∆t2 AΓI (α1 I − ∆tA1 )−1 AIΓ , (4.3)
(2) (2) (2) (2) (2) (2)

S(2) = (MΓΓ + ∆tAΓΓ ) − (MΓI + ∆tAΓI )(M2 + ∆tA2 )−1 (MIΓ + ∆tAIΓ ). (4.4)
One computes S(1) and S(2) by inserting the corresponding matrices specified above
in (4.3) and (4.4) obtaining
3λ1 ∆t λ2 ∆t2
S(1) = 2
I − 1 4 (4ETN − ETN −1 )(α1 I − ∆tA1 )−1 EN , (4.5)
2∆x 2∆x

(2) 1 5 1 λ2 ∆t 1 1
S = α2 tridiag − , , − + tridiag − , 2, −
24 12 24 ∆x2 2 2
(4.6)
λ2 ∆t λ2 ∆t
− α2 N2 − 2
I ET1 (M2 + ∆tA2 )−1 E1 α2 N2 − I .
∆x ∆x2
In the two-dimensional case, the iteration matrix Σ is a matrix of size N × N .
This makes the iteration matrix Σ difficult to compute for several reasons. First of all,
the matrices α1 I − ∆tA1 and M2 + ∆tA2 are sparse block tridiagonal matrices, and
consequently, their inverses are not straight forward to compute. A block-by-block
algorithm for inverting a block tridiagonal matrix is explained in [36]. However, the
algorithm is based on the iterative application of the Schur complement [46], and it
results in a sequence of block matrices and inverses of block matrices that we did not
find possible to compute exactly. Moreover, the diagonal blocks of α1 I − ∆tA1 and
M2 + ∆tA2 are tridiagonal but their inverses are full matrices [18].
Due to these difficulties, we propose here to approximate Σ. One can observe that
α1 I − ∆tA1 and M2 + ∆tA2 are strictly diagonally dominant matrices, and therefore,
we propose to approximate them by their block diagonal. Thus,
−1 !
2λ21 ∆t2 λ1 ∆t α1 ∆x2 + 4λ1 ∆t λ1 ∆t

(1) 3λ1 ∆t
S ≈ I − tridiag − 2 , ,− 2 , (4.7)
2∆x2 ∆x4 ∆x ∆x2 ∆x

(2) 1 5 1 λ2 ∆t 1 1
S ≈ α2 tridiag − , , − + tridiag − , 2, −
24 12 24 ∆x2 2 2
2 (4.8)
2

α2 ∆x + 12λ2 ∆t
− tridiag (b, a, b)−1 ,
12∆x2
with
5α2 ∆x2 + 24λ2 ∆t α2 ∆x2 + 12λ2 ∆t

a= , b=− .
6∆x2 12∆x2
Now, we compute the eigenvalues of the proposed approximations of S(1) and S(2) .
The eigenvalues of a tridiagonal Toeplitz matrix are known and given e.g. in [32, pp.
514-516]:
3λ1 ∆t(α1 ∆x2 + 2λ1 ∆t(2 − cos(jπ∆x))) − 4λ21 ∆t2

µ1j = , (4.9)
2∆x2 (α1 ∆x2 + 2λ1 ∆t(2 − cos(jπ∆x)))
2(α2 ∆x2 (5 − cos(jπ∆x)) + 12λ2 ∆t(2 − cos(jπ∆x)))2 − (α2 ∆x2 + 12λ2 ∆t)2
µ2j = ,
24∆x2 (α2 ∆x2 (5 − cos(jπ∆x)) + 12λ2 ∆t(2 − cos(jπ∆x)))
(4.10)
for j = 1, .., N . Here µ1j are the eigenvalues of the approximation of S(1) and µ2j the
eigenvalues of the approximation of S(2) . Note that the eigenvectors are common for
symmetric tridiagonal Toeplitz matrices, i.e, the eigenvectors do not depend on the
specific entries.
Thus, we obtain an estimate of the spectral radius of the iteration matrix Σ:
−1
−1
µ1
ρ(Σ) = ρ S(2) S(1) = µmax S(2) S(1) ≈ 21
µ1
2

12(α2 ∆x (5 − cos(π∆x)) + 12λ2 ∆t(2 − cos(π∆x)))
=
α1 ∆x2 + 2λ1 ∆t(2 − cos(π∆x))
3λ1 ∆t(α1 ∆x2 + 2λ1 ∆t(2 − cos(π∆x))) − 4λ21 ∆t2

· := σ.
2(α2 ∆x2 (5 − cos(π∆x)) + 12λ2 ∆t(2 − cos(π∆x)))2 − (α2 ∆x2 + 12λ2 ∆t)2
(4.11)
Furthermore, computing the limits of (4.11) when ∆t → 0 and ∆x → 0 we get
(12α2 ∆x2 (5 − cos(π∆x))) · 0

lim ρ(Σ) ≈ lim σ = = 0, (4.12)
∆t→0 ∆t→0 α1 ∆x2 (2α22 ∆x4 (5 − cos(π∆x))2 − α22 ∆x4 )
122 λ2 ∆t(6λ21 ∆t2 − 4λ21 ∆t2 ) λ1

lim ρ(Σ) ≈ lim σ = = =: δ. (4.13)
∆x→0 ∆x→0 2λ1 ∆t(2 · 122 λ22 ∆t2 − 122 λ22 ∆t2 ) λ2
Therefore, from the result obtained in (4.13) we expect that strong jumps in the
physical properties of the materials placed in Ω1 and Ω2 will imply fast convergence.
This is the case when modelling thermal fluid structure interaction, where often a fluid
with low thermal conductivity and density is coupled with a structure having higher
thermal conductivity and density.
When c → ∞, ∆t ∆x2 and (3.4) matches with the asymptotic computed in
(4.13). However, when c < 1, the semicontinuous analysis fails and the discrete
analysis just presented fills the gap.
4.1 Numerical Results

In this section we present a set of numerical experiments designed to show how the va-
lidity of the approximation of ρ(Σ) as an estimator for the rates of the coupled problem
formulated above. We also show that the theoretical asymptotics deduced in (4.12) and
(4.13) match with the numerical experiments.
We first compare the semidiscrete estimator β with the discrete estimator σ and the
experimental convergences rates. The experimental convergence rates CR are com-
puted with respect to a reference solution uref over the whole domain Ω using the
formula
ku3 − uref k2
CR = ,
ku2 − uref k2
where u2 and u3 are the second and third iterates of the Dirichlet-Neumann iteration.
Figure 8 shows a comparison between β and σ. On the left we plot β, σ and the
experimental convergence rates with c 1 and on the right we plot the same but
with c 1. We can conclude that the estimator for the convergence rates presented
in the previous section is minimally better than the one proposed by the semidiscrete
analysis in [24] when c 1. Moreover, when c 1 our estimator also predicts the
rates accurately and the semidiscrete estimator deviates.
0.3 0.3
<
Conv. Rate
0.25 - 0.25
0.2 0.2
0.15 0.15
0.1 0.1
<
0.05 0.05 Conv. Rate
-
0 0
0 0.002 0.004 0.006 0.008 0.01 0 0.1 0.2 0.3 0.4 0.5
"t "t
(a) c 1 (b) c 1
Figure 8. Here, D1 = 1, D2 = 0.5, λ1 = 0.3 and λ2 = 1. The circles correspond to σ, the

crosses to the experimental convergence rates and the continuous line to β. ∆x = 1/20
and on the left the curves are restricted to the discrete values ∆t = 1e − 2/50, 2 · 1e −
2/50, ..., 50 · 1e−2/50 and on the right to the values ∆t = 1e−2, 2 · 1e−2, ..., 50 · 1e−2.
We now want to illustrate how the formula (4.11) predicts the convergence rates
and tends to the limits computed previously. To this end, we present two real data
examples. We consider here the thermal interaction between air at 273K with steel at
900K and water at 283K with steel at 900K. Physical properties of the materials and
resulting asymptotics for these two cases are shown in table 1 and 2 respectively.
Table 1. Physical properties of the materials. λ is the thermal conductivity, ρ the density,
C the specific heat capacity and α = ρC.
Material λ (W/mK) ρ (kg/m3 ) C (J/kgK) α (J/K m3 )

Air 0.0243 1.293 1005 1299.5
Water 0.58 999.7 4192.1 4.1908e6
Steel 48.9 7836 443 3471348
Figures 9 and 10 show the convergence rates for the interactions between air-steel
and water-steel respectively. On the left we always plot the rates with respect to the
Table 2. The limits of the convergence rates when ∆t → 0 and ∆x → 0.
Case ∆t → 0 ∆x → 0
Air-Steel 0 4.9693e-4
Water-Steel 0 0.0119
variation of ∆t and for a fixed ∆x. On the right we plot the behaviour of the rates for a
fixed ∆t and varying ∆x.
We observe from figures 9 and 10 that the approximation σ predicts the convergence
rates quite well because the difference with respect to the experimental rates is really
small. Moreover, the rates in 9a and 10a tend to 0 as predicted in (4.12) and the rates
in 9b and 10b tend to δ as predicted in (4.13).
10 -3 10 -3
<
Conv. Rate
<
Conv. Rate
-4 -4 /
10 10
log
log
10 -5 10 -5
10 -6 10 -6
10 -1 10 0 10 1 10 -2 10 -1 10 0
log( "t) log( "x)
(a) The curves are restricted to the discrete values (b) The curves are restricted to the discrete values
∆t = 10/40, 2 · 10/40, ..., 40 · 10/40 and ∆x = ∆x = 1/3, 1/4, ..., 1/60 and ∆t = 10.
1/20.
Figure 9. Air-Steel thermal interaction with respect ∆t on the left and ∆x on the right.
From figure 9 we can observe that the convergence rates are really fast (factor of
∼ 1e−4) when there exist strong jumps in the coefficient of the materials. For instance,
when performing the thermal coupling between air and steel the Dirichlet-Neumann
iteration only needs two iterations to achieve a tolerance of 1e − 10.
5 Time Adaptive Methods

The standard in FSI is to use fixed time step sizes. We find that it is much more efficient
to use time adaptive methods. Thus the question arises how such a method can be
implemented in a partitioned way. The method described here has been suggested in
10 -2 10 -1
< <
Conv. Rate Conv. Rate
/
10 -2
10 -3
log
log
10 -3
10 -4
10 -4
10 -5 10 -5
10 -1 10 0 10 1 10 -2 10 -1 10 0
log( "t) log( "x)
(a) The curves are restricted to the discrete values (b) The curves are restricted to the discrete values
∆t = 10/40, 2 · 10/40, ..., 40 · 10/40 and ∆x = ∆x = 1/3, 1/4, ..., 1/60 and ∆t = 10.
1/20.
Figure 10. Water-Steel thermal interaction with respect ∆t on the left and ∆x on the right.
α α 0
√
1 1−α α α = 1 − 2/2
√
bi 1−α α α̂ = 2 − 54 2
√
b̂i 1 − α̂ α̂ α − α̂ = −1 + 34 2
bi − b̂i α̂ − α α − α̂
Table 3. Butcher array for SDIRK2.
[10, 9]. Basis is a master program that does the global time stepping and calls functions
available in the subsolvers.
Here, an A-stable singly diagonally implicit Runge-Kutta method (SDIRK) is used.
Consider an autonomous initial value problem
u̇(t) = f (u(t)), u(0) = u0 .
An SDIRK method is then defined as

i
X
i n
U = u + ∆tn aij f (Uj ), i = 1, ..., s (5.1)
j=1
s
X
un+1 = un + ∆tn bi f (Uk )
i=1
with given coefficients aij , bj , cj . Here, we use the two stage method SDIRK2, which
is defined by the coefficients in the Butcher array in Table 3. For this method, asi = bi ,
which implies that the last line is superfluous (first-is-last-property). Besides saving
computational effort, this guarantees L-stability. The vectors
ki = f (Ui )
are called stage derivates and s is the number of stages. Since the starting vector
i−1
X
si = un + ∆tn aij kj , i = 1, ..., s − 1
j=1
is known, (5.1) is just a sequence of implicit Euler steps.

In the following it is assumed that at time tn , the step size ∆tn is prescribed globally.
Applying a DIRK method to equations (2.9)-(2.10) results in the coupled system of
equations to be solved at Runge-Kutta stage i:
ui − sui − ∆tn aii h(ui , ΘΓi ) = 0, (5.2)
f
[M + ∆tn aii A(Θ)]Θi − MsΘ Γ
i − qb − qb (ui ) = 0. (5.3)
Here, sui and sΘi are the given starting vectors in the subproblems.
The coupled equations (5.2)-(5.3) are then solved using a Dirichlet-Neumann cou-
pling, as would be done for the implicit Euler method. Following the analysis just
presented, temperature is prescribed for the equation with smaller heat conductivity,
here the fluid, and heat flux is given on Γ for the structure. Choosing these conditions
the other way around leads to an unstable scheme. The canonical starting guess for the
Dirichlet-Neumann iteration at stage i is the starting vector si .
Regarding implementation, it is a safe assumption that both the fluid and the solid
solver are able to carry out time steps of implicit Euler type. Then the master program
of the FSI procedure can be extended to SDIRK methods very easily. The master
program just has to call at stage i in iteration k the backward Euler routines with time
step size aii ∆tn and starting vectors si and with the boundary data appropriate for
iteration k.
To obtain time adaptivity, the technique of embedded Runge-Kutta methods is used.
Thereby, the local error l̂ is estimated by the solvers separately:
2
(i)
X
(i) (i)
l̂ ≈l = (bj − b̂j )kj , i = 1, 2, (5.4)
j=1
(i)
where kj are the stage derivatives on the subdomains for each Runge-Kutta stage.
The computation of (5.4) will typically not be implemented in the subsolvers and has
to be added to both codes, as well as that all stage derivatives have to be stored by the
subsolvers. These error estimates are then reported back to the master program and
aggregated in the form of a weighted scaled norm [9]:
v
u1 n
u X r
li 1
klkW SN = t = (n1 kl(1) k2W SN + n2 kl(2) k2W SN ).
n T OLuj + T OL n
j=1
(5.5)
Thus, the master solver either needs to know the number of unknowns n1 and n2 in the
subsolvers, or instead the quantities ni kl(i) k2W SN have to be provided to it. Based on
this, the new time step is chosen to comply with a user defined error tolerance T OL
for the time integration:
−1/k
∆tn+1 = ∆tn · klkW SN (5.6)
Finally, if the possibility of rejected time steps is taken into account, the current
solution pair (u, Θ) has to be stored as well.
5.1 Extrapolation from Time Integration

To find better starting values for iterative processes in implicit time integration schemes,
it is common to use extrapolation based on knowledge about the trajectory of the so-
lution of the initial value problem [15, 33]. In the spirit of partitioned solvers, it was
suggested in [8] to use extrapolation of the interface temperatures only. On top, this
strategy could be used as well within the subsolvers, which we will not consider here
and leave that to the discretion of the developers of those codes. We now present
extrapolation methods for SDIRK2 from [8].
At the first stage, we have the old time step size ∆tn−1 with value Θn−1 and the
current time step size ∆tn with value Θn . We are looking for the value Θ1n at the next
stage time tn + c1 ∆tn . Linear extrapolation results in

c1 ∆tn c1 ∆tn
Θ1n ≈ Θn + c1 ∆tn (Θn − Θn−1 )/∆tn−1 = 1 + Θn − Θn−1 . (5.7)
∆tn−1 ∆tn−1
At the second stage, we linearly extrapolate Θn at tn and Θ1n at tn + c1 ∆t to obtain

1 1 1
Θn+1 ≈ Θn + ∆tn (Θn − Θn )/(c1 ∆tn ) = 1 − Θn + Θ1n . (5.8)
c1 c1
6 Choosing Tolerances
In the time adaptive setting, a core input parameter is a tolerance TOL. This, based on
the error estimate, steers the time step size according to (5.5)-(5.6). For this to work,
it is imperative that the iteration error is smaller than the time integration error, so as
not to destroy the error estimate of the latter. The Dirichlet-Neumann iteration is for-
mulated as a fixed point iteration in the interface unknowns with termination criterion
(2.15). The tolerance there thus has to be chosen that the iteration error not only on the
interface, but on the whole domain is smaller than the time integration error. At the
same time, oversolving should be avoided to prevent unnecessary computational cost.
A good rule of thumb is to divide the tolerance for the time integration procedure by
five [38], ending up with
τ = T OL/5. (6.1)
On the assumption of exact solves in the subsolvers, errors in the interface will
nevertheless be translated into errors in the subproblems. The boundary condition from
the interface enters the right hand side of the nonlinear equations on the subproblem.
Thus, by way of the inverse operator, errors in the interface values are translated into
errors in the subdomains. For parabolic problems as in heat transfer, these are bounded
operators and thus the error transfer can be controlled. Right now, we assume that it is
of order one and do not adjust (6.1) further. However, a more careful analysis of this
point to get better estimates is required, for example using discrete trace inequalities.
Actually, it is common to use iterative solvers for the subproblems as well, typically
Newton or Multigrid methods (compare [6] for an overview on the state of the art). So
the question needs to be raised again, what type of termination criteria to choose and
what kind of tolerance. For a nonlinear subproblem of the form
F(u) = 0,
the standard would be to use a relative termination criterion
kF(uk )k ≤ τr kF(u0 )k
with relative tolerance τr . The errors introduced from these systems into the fixed
point iteration will be on the order of the residual. Thus, it makes sense to use (6.1)
for the tolerances in the subsolvers as well.
In the linear case, the same reasoning applies when we use the relative termination
criterion
kA(xk )xk+1 − bk ≤ τr kbk. (6.2)
However, it turns out that this has some pitfalls and that in fact the use of the nonstan-
dard relative criterion
kA(xk )xk+1 − bk ≤ τr kA(xk )xk − bk (6.3)
can speed up the iteration significantly. In the following analysis, taken from [7], an
absolute criterion
kA(xk )xk+1 − bk ≤ τa , (6.4)
with an absolute tolerance τa is possible as well. The first criterion, in case of conver-
gence of the iteration, leads to an automatic tightening of the tolerance over the course
of the iteration.
To determine good tolerances, we will now analyze the errors arising from these
more careful.
6.1 Perturbed Nested Fixed Point Iteration

Consider two functions F : Ω1 → Ω2 and S : Ω2 → Ω1 with Ω1 , Ω2 ⊂ Rn closed
and the fixed point equation
x = S(F(x)) (6.5)
which we assume to have a unique solution x∗ . Let then both the evaluation of F and
of S be perturbed, namely S by δk and F by k , and consider the iteration:
xk+1 = S(F(xk ) + k ) + δk . (6.6)
We assume that this iteration is well defined and that this sequence has the limit x .
Then, we have the following theorem [7].
Theorem 6.1. Let F and S be Lipschitz continuous with Lipschitz constants LF and
LS , respectively. Assume that LF LS < 1. Then we have, if k = δk = for all k, that
1 + LS
kx − x∗ k ≤ . (6.7)
1 − LS LF
In the case k = and δk = δ, we obtain
LS + δ
kx − x∗ k ≤ . (6.8)
1 − LS LF
Finally, x = x∗ if and only if both δk and k converge to zero.
Proof: We have due to the Lipschitz continuity
kxk+1 − x∗ k = kS(F(xk ) + k ) + δk − x∗ k = kS(F(xk ) + k ) + δk − S(F(x∗ ))k

≤ LS kF(xk ) − F(x∗ ) + k k + δk ≤ LS LF kxk − x∗ k + LS k + δk
≤ (LS LF )2 kxk−1 − x∗ k + L2S LF k−1 + LS LF δk−1 + LS k + δk
   
k k
Lj+1 j
LjS LjF δk−j 
X X
≤ (LS LF )k+1 kx0 − x∗ k +  S LF k−j
+
j=0 j=0
and thus in the limit xk+1 → x ,

k
X k
X
kx − x∗ k ≤ LS lim (LS LF )j k−j + lim (LS LF )j δk−j (6.9)
k→∞ k→∞
j=0 j=0
For a constant perturbation overall, e.g. k = δk = for all k, we obtain in the limit
k
X 1 + LS
kx − x∗ k ≤ (1 + LS ) lim (LS LF )j = ,
k→∞ 1 − LS LF
j=0
which proves the inequality (6.7). If we have constant but separate perturbations and
δ of S and F, we obtain (6.8) from
k k
X X LS + δ
kx − x∗ k ≤ LS lim (LS LF )j + δ lim (LS LF )j = .
k→∞ k→∞ 1 − LS LF
j=0 j=0
In the general case, due to positivity, the right hand side of (6.9) is zero if and only
if both k and δk are such that for φk = k or φk = δk ,
k
X
lim (LS LF )j φk−j = 0.
k→∞
j=0
This is the case if and only if both k and δk converge to zero.
6.2 Application: Dirichlet-Neumann Coupling for Transmission

Problem
As an application of the theory just presented, we consider the steady transmission
problem, where the Laplace equation with right hand side f (x, y) on a domain Ω is
cut into two domains Ω = Ω1 ∪ Ω2 using transmission conditions at the interface
Γ = Ω1 ∩ Ω2 :
∆um (x, y) = f (x, y), (x, y) ∈ Ωm ⊂ R2 , m = 1, 2,

um (x, y) = 0, (x, y) ∈ ∂Ωm Γ, (6.10)
u1 (x, y) = u2 (x, y), (x, y) ∈ Γ,
∂u2 (x, y) ∂u1 (x, y)
=− , (x, y) ∈ Γ.
n2 n1
We now employ a standard Dirichlet-Neumann iteration to solve it. This corre-
sponds to alternately solving the problems
A1 uk+1
1 = b1 (uk2 ) (6.11)
and
A2 uk+1
2 = b2 (uk+1
1 ), (6.12)
where problem (6.11) originates from a linear discretization of the transmission prob-
lem (6.10) on Ω1 only with Dirichlet data on Γ given by uk2 on the coupling interface
and problem (6.12) corresponds to a linear discretization of (6.10) on Ω2 only with
Neumann data on Γ given by the discrete normal derivative of u1 on Γ.
By considering (6.11)-(6.12) as a nested iteration, we obtain a fixed point formula-
tion
uΓ = S(F(uΓ ))
where uΓ = u2 |Γ , F = DnΓ A−1 −1

1 b1 (uΓ ) and S = PΓ A2 b2 (u1 ). Hereby DnΓ is the
matrix that computes the discrete normal derivatives in Ω1 on Γ and PΓ is the discrete
trace operator with respect to Γ.
In practice, the linear equation systems are solved iteratively, typically using the
conjugate gradient method (CG) up to a relative tolerance of τ . Thus, we obtain a
perturbed nested fixed point iteration of the form (6.6) and the question is now again
if we can quantify this perturbation. We have
−1
uk+1 k
1 = A1 b1 (uΓ ) + k (6.13)
and
−1
uk+1 k+1
2 = A2 b2 (u1 ) + δk . (6.14)
For the iteration (6.13) we obtain
−1 −1
kk k = kuk+1 k k+1 k
1 − A1 b1 (uΓ )k ≤ kA1 kkA1 u1 − b1 (uΓ )k.
Again, the second factor is what is tested in the termination criterion. In the case of
the relative criterion (6.3), this results in
kA1 uk+1 k k k
1 − b1 (uΓ )k ≤ τr kA1 u1 − b1 (uΓ )k.
We thus have managed to shift the index in A1 uk1 , but not in b1 (ukΓ ). For this we add
zero and use the triangle inequality:
τr kA1 uk1 − b1 (ukΓ )k ≤ τr (kA1 uk1 − b1 (uk−1 k−1 k

Γ )k + kb1 (uΓ ) − b1 (uΓ )k).
We can now repeat this argument on kA1 uk1 − b1 (ukΓ )k, thereby multiplying again
with τr and adding a difference of right hand sides:
 
k
τrj kb1 (uk−j k−j+1 
X
kk k ≤ kA−1  k 1 0
1 k τr kA1 u1 − b1 (uΓ )k + Γ ) − b1 (uΓ )k .
j=1
The right hand sides are equal except at the boundary and we have
kb1 (uk−j k−j+1

Γ ) − b1 (uΓ )k = 1/∆x2 kuΓk−j − uk−j+1
Γ k.
Furthermore, the last system is tested against the unperturbed initial data leading to
 
k
τrj /∆x2 kuk−j − uk−j+1
X
kk k ≤ kA−1  k+1 kA1 u01 − b1 (u0Γ )k +
1 k τr Γ Γ k .
j=1
For k to infinity, the first term on the right converges to zero if and only if τr < 1 and
from Lemma 1 we know that the second term on the right converges to zero if and
only if kuk−j
Γ − uk−j+1
Γ k converges to zero, which is the case if τr < 1 and if ukΓ is
a convergent sequence. Note that these are the perturbed iterates of the twin iteration.
Thus, we have to assume that that iteration converges.
If we choose the absolute termination criterion (6.4) or the relative one based on the
right hand side (6.2), we obtain a bound of the form
kk k ≤ kA−1 k
1 kτr kb(uΓ )k,
respectively
kk k ≤ kA−1
1 kτa .
Here, we cannot make a statement on the limit of k .

In the second case, meaning the iteration with Neumann data (6.14), we obtain
−1 −1
kδk k = kuk+1 k+1 k+1 k+1
2 − A2 b2 (u1 )k ≤ kA2 kkA2 u2 − b2 (u1 )k
and analogous arguments produce the same results for δk , meaning that we obtain
convergence to zero under the assumptions that τr < 1 and that the perturbations k
are convergent. By theorem 6.1, we then have that when using the relative criterion we
obtain convergence to the exact solution for any τr < 1.
6.3 Testcase: Transmission Problem

We now consider the transmission problem (6.10). Specifically, we use Ω1 = [0, 1] ×
[0, 1], Ω2 = [1, 2] × [0, 1] and
π π
f (x, y) = sin πy 2 (π cos x2 − π 2 x2 sin x2 ) (6.15)
2 2
π 2
+ sin x (2π cos πy 2 − 4π 2 y 2 sin πy 2 ).
2
This was chosen such that the solution is
π
u(x, y) = sin πy 2 sin x2 , (6.16)
2
which satisfies the boundary conditions.
We discretize this problem using central differences with a constant mesh width of
∆x = ∆y. As initial guess for the Dirichlet-Neumann procedure, we employ a vector
of all zeros. All linear systems are solved using a multigrid method. Specifically, we
use a V cycle with exact solves for systems of dimension smaller than 16. Restriction
and prolongation are standard full weighting, respectively bilinear interpolation [42],
appropriately extended for the Neumann problem. The smoother is a damped Jacobi
method with ω = 2/3, which is applied thrice both as a pre- and a postsmoother. The
exact solution and the discrete solution with ∆x = 1/64 can be seen in Figure 11.
We now look at the convergence properties of the fixed point schemes for different
termination criteria. Thus we choose T OL = 1e − 10 and ∆x = 1/64 and then run the
Figure 11. Exact and discrete solution with ∆x = 1/64
iteration for different tolerances in multigrid. In the case of the termination criterion
(6.3), the iteration recovers the reference solution for any tolerance τr , as predicted
by the theory. For the relative termination criterion (6.2) and the absolute termination
criterion (6.4) the situation is more complicated. We observe that for a given tolerance,
kuk+1
Γ − ukΓ k does not converge to zero, but instead oscillates around a value that is
proportional to the tolerance in multigrid. The error kuΓ − u∗Γ k2 is then about half that
value. This means that in practice, these schemes behave as if there were a constant
perturbation and that the error can be controlled by choosing an appropriately small
tolerance in the linear solver. Note that (6.2) is implemented in the MATLAB version
of CG and that thus, a native implementation of the Dirichlet-Neumann iteration with
CG in MATLAB will produce questionable results if not used with care.
4 5
x 10 x 10
7 2
Relative to r Relative to r
Relative to b 1.8 Relative to b
6 Absolute Absolute
1.6
5 1.4
1.2
MG Iterations
MG Iterations
4
1
3
0.8
2 0.6
0.4
1
0.2
0 0
−1 −2 −3 −4 −5 −6 −7 −8 −1 −2 −3 −4 −5 −6 −7
10 10 10 10 10 10 10 10 10 10 10 10 10 10 10
Tolerance Tolerance
Figure 12. Total multigrid iterations over tolerance for different termination criteria for
the transmission problem with rhs (6.15) and ∆x = 1/128 (left) and ∆x = 1/256 (right).
We now compare the different termination criteria with regards to efficiency for
values of TOL more relevant in practice. Hereby, we assume that the user wants to
have a solution that is TOL close to the exact one. Based on the theory discussed here,
there are three choices: Using (6.3) with τr = 1e − 1 independent of TOL, using
(6.4) with τa = T OL or (6.2) with τr = T OL, meaning that for the latter ones, we
have to solve the inner iteration more accurately the more accurate we want the outer
one. Hereby, we choose ∆x = 1/128 and ∆x = 1/256. As turns out, the error is
in all cases about 0.4 · T OL. Since costs in the fixed point iterations not arising from
solving the linear systems is negligible and the cost of one multigrid iteration is fixed,
we compare the total number of multigrid iterations needed until termination. The
results are depicted in Figure 12. As we can see, the scheme corresponding to (6.2) is
the slowest. This is more pronounced on the finer grid, where it is about a factor of
two slower than the other schemes. The only scheme where the computational effort
has a linear behavior with regards to tolerance is that corresponding to (6.3). Thus,
while slower for large tolerances, it beats the scheme with (6.4) at some point.
4 5
x 10 x 10
7 2
Relative to r Relative to r
Relative to b Relative to b
1.8
6 Absolute Absolute
1.6
5
1.4
MG Iterations
MG Iterations
4 1.2
1
3
0.8
2
0.6
1 −1 −2 −3 −4 −5 −6 −7
0.4 −1 −2 −3 −4 −5 −6
10 10 10 10 10 10 10 10 10 10 10 10 10
Tolerance Tolerance
Figure 13. Total multigrid iterations over tolerance for different termination criteria for
the transmission problem with rhs (6.15) multiplied by 100 and ∆x = 1/128 (left) and
∆x = 1/256 (right).
An absolute termination criterion can lead to oversolving. This is illustrated in the

next test, where we just multiply the right hand side (6.15) with 100. The results are
shown in Figure 13. Again, criterion (6.2) leads the slowest scheme. However, the
scheme corresponding to (6.3) is now faster than that with criterion (6.4) for T OL <
1e − 1.
7 Numerical Results
The following results illustrate the performance of the complete solver described so
far for realistic test cases of a steel-air coupling. Thus, the models from Section 2 are
employed and discretized using FEM in the structure and FVM in the fluid. In partic-
ular, the DLR TAU-Code is employed [19], which is a cell-vertex-type finite volume
method with AUSMDV as flux function and a linear reconstruction to increase the or-
der of accuracy. The finite element code uses quadratic finite elements [47] and is the
inhouse code Native of the Institute for Static and Dynamic at the University of Kas-
sel. The coupling between the solvers is done using the Component Template Library
(CTL) of the University of Braunschweig [30]. The following numerical results are
taken from [8].
7.1 Flow over a Plate

As a first test case, the cooling of a flat plate resembling a simple work piece is consid-
ered. The work piece is initially at a much higher temperature than the fluid and then
cooled by a constant air stream, see Figure 14.
Figure 14. Test case for the coupling method
The inlet is given on the left, where air enters the domain with an initial velocity of
Ma∞ = 0.8 in horizontal direction and a temperature of 273 K. Then, there are two
succeeding regularization regions of 50 mm to obtain an unperturbed boundary layer.
In the first region, 0 ≤ x ≤ 50, symmetry boundary conditions, vy = 0, q = 0, are
applied. In the second region, 50 ≤ x ≤ 100, a constant wall temperature of 300 K is
specified. Within this region the velocity boundary layer fully develops. The third part
is the solid (work piece) of length 200 mm, which exchanges heat with the fluid, but is
assumed insulated otherwise, thus qb = 0. Therefore, Neumann boundary conditions
are applied throughout. Finally, the fluid domain is closed by a second regularization
region of 100 mm with symmetry boundary conditions and the outlet.
Regarding the initial conditions in the structure, a constant temperature of 900 K
at t = 0 s is chosen throughout. To specify reasonable initial conditions within the
fluid, a steady state solution of the fluid with a constant wall temperature Θ = 900 K
is computed.
The grid is chosen cartesian and equidistant in the structural part, where in the fluid
region the thinnest cells are on the boundary and then become coarser in y-direction
(see Figure 15). To avoid additional difficulties from interpolation, the points of the
primary fluid grid, where the heat flux is located in the fluid solver, and the nodes of
the structural grid are chosen to match on the interface Γ.
We now compare the different schemes for a whole simulation of 100 seconds real
time. If not mentioned otherwise, the initial time step size is ∆t = 0.5s. To first give
an impression on the effect of the time adaptive method, we look at fixed time step
(a) Entire mesh (b) Mesh zoom
Figure 15. Full grid (left) and zoom into coupling region (right)
T OL Fixed time step size Time adapt. Time adapt. with lin. extr.
10−2 ∆t = 5s 64 31 19
10−3 ∆t = 5s 82 39 31
10−4 ∆t = 0.5s 802 106 73
10−5 ∆t = 0.5s 1014 857 415
Table 4. Total number of iterations for 100 secs of real time. Fixed time step sizes versus
adaptive steering. For time adaptive calculations, the initial time step is ∆t0 = 0.5s.
versus adaptive computations in Table 4. Thus, the different tolerances for the time
adaptive case lead to different time step sizes and tolerances for the nonlinear system
over the course of the algorithm, whereas in the fixed time step size, they steer only
how accurate the nonlinear systems are solved. For the fixed time step case, we chose
∆t = 0.5s and ∆t = 5s, which roughly corresponds to an error of 10−2 and 10−3 ,
respectively 10−4 . Thus, computations in one line of table 4 correspond to similar
errors. As it can be seen, the time adaptive method is in the worst case a factor two
faster and in the best case a factor of eight. Thus the time adaptive computation serves
from now on as the base method for the construction of a fast solver.
Finally, we consider extrapolation based on the time integration scheme. As can
be seen in Table 4, linear extrapolation speeds up the computations between 20% and
50%. Overall, we are thus able to simulate 100 seconds of real time for this problem
for an engineering tolerance using only 19 calls to fluid and the structure solver each.
7.2 Cooling of a Flanged Shaft

To illustrate a realistic example of gas quenching, the cooling of a flanged shaft by
cold high pressured air is used as a second test case. The complete process consists of
the inductive heating of a steel rod, the forming of the hot rod into a flanged shaft, a
transport to a cooling unit and the cooling process. Here, we consider only the cooling,
meaning that we have a hot flanged shaft that is cooled by cold high pressured air com-
ing out of small tubes [45]. We perform a two dimensional cut through the domain and
assume symmetry along the vertical axis, resulting in one half of the flanged shaft and
two tubes blowing air at it, see figure 16. Since the air nozzles are evenly distributed
around the flanged shaft, we use an axisymmetric model in the structure. The heat flux
from the two-dimensional simulation of the fluid at the boundary of the flanged shaft
is impressed axially symmetrical on the structure.
Figure 16. Sketch of the flanged shaft
We assume that the air leaves the tube in a straight and uniform way at a Mach
number of 1.2. Furthermore, we assume a freestream in x-direction of Mach 0.005.
This is mainly to avoid numerical difficulties at Mach 0, but could model a draft in the
workshop. The Reynolds number is Re = 2500 and the Prandtl number P r = 0.72.
The grid consists of 279212 cells in the fluid, which is the dual grid of an unstruc-
tured grid of quadrilaterals in the boundary layer and triangles in the rest of the domain,
and 1997 quadrilateral elements in the structure. It is illustrated in Figure 17.
To obtain initial conditions for the subsequent tests, we use the following procedure:
We define a first set of initial conditions by setting the flow velocity to zero throughout
and choose the structure temperatures at the boundary points to be equal to tempera-
tures that have been measured by a thermographic camera. Then, setting the y-axis on
the axis of revolution of the flange, we set the temperature at each horizontal slice to
the temperature at the corresponding boundary point. Finally, to determine the actual
initial conditions, we compute 10−5 seconds of real time using the coupling solver
with a fixed time step size of ∆t = 10−6 s. This means, that the high pressured air
is coming out of the tubes and the first front has already hit the flanged shaft. This
solution is illustrated in Figure 18 (left).
Now, we compute 1 second of real time using the time adaptive algorithm with
different tolerances and an initial time step size of ∆t = 10−6 s. This small initial step
(a) Entire mesh (b) Mesh zoom
Figure 17. Full grid (left) and zoom into shaft region (right)
size is necessary to prevent instabilities in the fluid solver. During the course of the
computation, the time step size is increased until it is on the order of ∆t = 0.1s, which
demonstrates the advantages of the time adaptive algorithm and reaffirms that it is this
algorithm that we need to compare to. In total, the time adaptive method needs 22,
41, 130 and 850 time steps to reach t = 1s for the different tolerances, compared to
the 106 steps the fixed time step method would need. The solution at the final time is
depicted in Figure 18 (right). As can be seen, the stream of cold air is deflected by the
shaft.
Finally, we consider extrapolation based on the time integration scheme. In Table 5,
the total number of iterations for 1 second of real time is shown. As before, the extrap-
olation causes a noticable decrease in the total number of fixed point iterations. The
speedup from linear extrapolation is between 18% and 34%, compared to the results
obtained without extrapolation.
T OL none lin.
10−2 51 42
10−3 126 97
10−4 414 309
10−5 2768 1805
Table 5. Total number of iterations for 1 sec of real time for different extrapolation
methods in time.
Figure 18. Temperature distribution in fluid and structure at t = 0s (left) and t = 1s

(right)
8 Summary and Conclusions

We considered numerical methods for unsteady thermal fluid structure iteraction. As
a basic paradigm, we employed a partitioned method with the use of the Dirichlet-
Neumann iteration. In this framework, a time adaptive SDIRK2 method with linear
extrapolation of initial values is a fast solver. We cite three reasons for this: the time
adaptivity, the very fast convergence rate of the Dirichlet-Neumann iteration for this
problem and a smart choice of tolerances in the various nested solvers.
The fast convergence rate can be explained by estimating the spectral radius of the
iteration matrix based on a model problem. This shows that the quotients of material
properties play a crucial role. An analysis of the tolerances to be chosen shows that
in the linear case, a nonstandard relative termination criterion based on the current
residual leads to better performance compared to a standard termination criterion and
to avoid pitfalls with possible wrong limits.
Bibliography
[1] Background on Tonight’s Launch, http://www.spacex.com/news/2015/12/21/
background-tonights-launch, Accessed: 2016-08-23.
[2] CRS-8 Launch and Landing, http://www.spacex.com/news/2016/04/09/
crs-8-launch-and-landing, Accessed: 2016-08-23.
[3] The Why and How of Landing Rockets, http://www.spacex.com/news/2015/06/
24/why-and-how-landing-rockets, Accessed: 2016-08-23.
[4] Y. Bazilevs, K. Takizawa and T. E. Tezduyar, Computational Fluid-Structure Interaction:
Methods and Applications, Wiley, 2013.
[5] O.O. Bendiksen, A new approach to computational aeroelasticity, AIAA Paper 91-09
(1991), 1712–1727.
[6] P. Birken, Numerical Methods for the Unsteady Compressible Navier-Stokes Equations,
Habilitation Thesis, University of Kassel, 2012.
[7] , Termination criteria for inexact fixed point schemes, Numer. Linear Algebra
Appl. (2015).
[8] P. Birken, T. Gleim, D. Kuhl and A. Meister, Fast Solvers for Unsteady Thermal Fluid
Structure Interaction, Int. J. Numer. Meth. Fluids 79(1) (2015), 16–29.
[9] P. Birken, K.J. Quint, S. Hartmann and A. Meister, Choosing norms in adaptive FSI
calculations, PAMM 10 (2010), 555–556.
[10] , A time-adaptive fluid-structure interaction method for thermal coupling, Comp.
Vis. in Science 13(7) (2011), 331–340.
[11] D.L. Brown et al., Applied Mathematics at the U.S. Department of Energy: Past, Present
and a View to the Future, tech. rep., Office of Science, U.S. Department of Energy, 2008.
[12] J.M. Buchlin, Convective Heat Transfer and Infrared Thermography, J. Appl. Fluid Mech.
3 (2010), 55–62.
[13] G.A. Davis and O.O. Bendiksen, Transonic Panel Flutter, AIAA Paper 93-1476 (1993).
[14] S. Deparis, M.A. Fernández and L. Formaggia, Acceleration of a fixed point algorithm for
fluid-structure interaction using transpiration conditions, M2AN 37(4) (2003), 601–616.
[15] P. Erbts and A. Düster, Accelerated staggered coupling schemes for problems of ther-
moelasticity at finite strains, Comp. & Math. with Appl. 64 (2012), 2408–2430.
[16] C. Farhat, CFD-based Nonlinear Computational Aeroelasticity, in Encyclopedia of Com-
putational Mechanics (2004), 459–480, ch. 13.
[17] C. Farhat, K.G. van der Zee and P. Geuzaine, Provably second-order time-accurate
loosely-coupled solution algorithms for transient nonlinear computational aeroelasticity,
Comput. Methods Appl. Mech. Engrg. 195 (2006), 1973–2001.
[18] C.M. Fonseca and J. Petronilho, Explicit inverses of some tridiagonal matrices, Linear
Algebra Appl. 325(1-3) (2001), 7–21.
[19] T. Gerhold, O. Friedrich, J. Evans and M. Galle, Calculation of Complex Three-
Dimensional Configurations Employing the DLR-TAU-Code, AIAA Paper 97-0167
(1997).
[20] M.B. Giles, Stability Analysis of Numerical Interface Conditions in Fluid-Structure Ther-
mal Analysis, Int. J. Numer. Meth. Fluids 25 (1997), 421–436.
[21] M. Guenther, A. Kvaerno and P. Rentrop, Multirate partitioned Runge-Kutta methods,
BIT 41 (2001), 504–514.
[22] H. Guillard and C. Farhat, On the significance of the geometric conservation law for
flow computations on moving meshes, Comput. Methods Appl. Mech. Engrg. 190 (2000),
1467–1482.
[23] U. Heck, U. Fritsching and Bauckhage K., Fluid flow and heat transfer in gas jet quench-
ing of a cylinder, Int. J. Numer. Methods Heat Fluid Flow 11 (2001), 36–49.
[24] W.D. Henshaw and K.K. Chand, A composite grid solver for conjugate heat transfer in
fluid-structure systems, J. Comput. Phys. 228 (2009), 2708–3741.
[25] M. Hinderks and R. Radespiel, Investigation of Hypersonic Gap Flow of a Reentry Nose-
cap with Consideration of Fluid Structure Interaction, AIAA Paper 6 (2006), 2708–3741.
[26] D. Kowollik, V. Tini, S. Reese and M. Haupt, 3D fluid-structure interaction analysis of
a typical liquid rocket engine cycle based on a novel viscoplastic damage model, Int. J.
Numer. Methods Engrg. 94 (2013), 1165–1190.
[27] D.S.C. Kowollik, P. Horst and M.C. Haupt, Fluid-structure interaction analysis applied to
thermal barrier coated cooled rocket thrust chambers with subsequent local investigation
of delamination phenomena, Progress in Propulsion Physics 4 (2013), 617–636.
[28] P. Le Tallec and J. Mouro, Fluid structure interaction with large structural displacements,
[29] R. Massjung, Discrete conservation and coupling strategies in nonlinear aeroelasticity,
[30] H. G. Matthies, R. Niekamp and J. Steindorf, Algorithms for strong coupling procedures,
[31] R.C. Mehta, Numerical Computation of Heat Transfer on Reentry Capsules at Mach 5,
AIAA-Paper 178 (2005).
[32] C.D. Meyer, Matrix Analysis and Applied Linear Algebra, 2000.
[33] H. Olsson and G. Söderlind, Stage value predictors and efficient Newton iterations in
implicit Runge-Kutta methods, SIAM J. Sci. Comput. 20 (1998), 185–202.
[34] A. Quarteroni and A. Valli, Domain Decomposition Methods for Partial Differential
Equations, Oxford Science Publications, 1999.
[35] K. J. Quint, S. Hartmann, S. Rothe, N. Saba and K. Steinhoff, Experimental validation
of high-order time integration for non-linear heat transfer problems, Comput. Mech. 48
(2011), 81–96.
[36] M.G. Reuter and J.C. Hill, An efficient, block-by-block algorithm for inverting a block
tridiagonal, nearly block Toeplitz matrix, Comp. Sci. Disc. 5 (2012).
[37] V. Savcenco, W. Hundsdorfer and J. G. Verwer, A multirate time stepping strategy for
stiff ordinary differential equations, BIT Numerical Mathematics 47 (2007), 137–155.
[38] G. Söderlind and L. Wang, Adaptive time-stepping and computational stability, J. Comp.
Appl. Math. 185 (2006), 225 – 243.
[39] P. R. Spalart and S. R. Allmaras, A One-Equation Turbulence Model for Aerodynamic
Flows, AIAA 30th Aerospace Science Meeting 92–0439 (1992).
[40] P. Stratton, I. Shedletsky and M. Lee, Gas Quenching with Helium, Solid State Phenom-
ena 118 (2006), 221–226.
[41] A. Toselli and O. Widlund, Domain Decomposition Methods - Algorithms and Theory,
Springer, 2004.
[42] U. Trottenberg, C. W. Oosterlee and S. Schüller, Multigrid, Elsevier Academic Press,
2001.
[43] A. van Zuijlen and H. Bijl, Implicit and explicit higher order time integration schemes
for structural dynamics and fluid-structure interaction computations, Comp. & Struct. 83
(2005), 93–105.
[44] A. van Zuijlen, A. de Boer and H. Bijl, Higher-order time integration through smooth
mesh deformation for 3D fluid-structure interaction simulations, J. Comp. Phys. 224
(2007), 414–430.
[45] U. Weidig, N. Saba and K. Steinhoff, Massivumformprodukte mit funktional gradierten
Eigenschaften durch eine differenzielle thermo-mechanische Prozessführung, WT-Online
(2007), 745–752.
[46] F. Zhang, The Schur complement and its applications, Springer, 2005.
[47] O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method, Butterworth Heinemann,
2000.
Author information
Philipp Birken, Lund University, Centre for the Mathematical Sciences, Numerical Analysis,
Box 118, 22100 Lund, Sweden.
E-mail: philipp.birken@na.lu.se
Azahar Monge, Lund University, Centre for the Mathematical Sciences, Numerical Analysis,
Box 118, 22100 Lund, Sweden.
E-mail: azahar.monge@na.lu.se

Numerical Methods For Unsteady Thermal Fluid Structure Interaction

Uploaded by

Copyright:

Available Formats

Numerical Methods For Unsteady Thermal Fluid Structure Interaction

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Numerical Methods For Unsteady Thermal Fluid Structure Interaction

Uploaded by

Copyright:

Available Formats

XXXX, 1–39 © De Gruyter YYYY

Numerical Methods for Unsteady Thermal

Keywords. Thermal Fluid Structure Interaction, Conjugate Heat Transfer, Dirichlet-Neumann

AMS classification. 65F10, 65L04, 65M22, 74F04, 74F10.

(a) Inductive heating. (b) Thermo-mechanical (c) Local air-cooling.

Figure 4. Sketch of the coupling surfaces. Figure taken from [27].

1.1 Partitioned Coupling Methods

Figure 5. Illustration of the FSI solver.

convergence can be improved using a relaxation procedure [34]. In fluid structure

4. Numerical results are included to illustrate the theoretical analysis. In section 5,

2 Treating Thermal FSI Problems

2.1 Fluid Model

The Reynolds stresses

are modelled using the Spallart-Allmaras model [39]. Furthermore, qf = (q1 , q2 )T

2.2 Structure Model

λ(Θ) = 40.1 + 0.05Θ − 0.0001Θ2 + 4.9 · 10−8 Θ3 (2.3)

cp1 (Θ) = 34.2e0.0026Θ + 421.15 (2.5)

cp2 (Θ) = 956.5e−0.012(Θ−900) + 0.45Θ. (2.6)

For the mass density one has ρ = 7836 kg/m3 .

2.3 Coupling Conditions

T (x, t) = Θ(x, t), x ∈ Γ, (2.7)

where T is the fluid temperature and Θ the structure temperature and

qf (x, t) · n(x) = qs (x, t) · n(x), x ∈ Γ. (2.8)

2.4 Discretization in Space

2.5 Time Discretization

un+1 − un + ∆tn h(un+1 , Θn+1

2.6 The Dirichlet-Neumann Method

uk+1 − un + ∆tn h(uk+1 , ΘkΓ ) = 0, (2.13)

where τ is a user defined tolerance.

3 A Model Problem: Coupled Heat Equations

3.1 Model Problem

where nm is the outward normal to Ωm for m = 1, 2.

3.2 Semidiscrete Analysis

3.3 Space Discretization

Figure 6. Splitting of Ω between finite volumes and finite elements.

(1) (1) (1) (1)

(2) (2) (2) (2)

Consequently, the equation

M̃u̇ + Ãu = 0, (3.10)

3.4 Time Discretization

Aun+1 = M̃un , (3.11)

(1) (2) (1) (2)

3.5 Fixed Point Iteration

(1),n+1,k+1 (1) (1) (1),n (1)

Âûk+1 = M̂un − bk , (3.13)

with some initial condition, here un+1,0

(m) (m) (m) (m) (m) (m)

The discretization matrices are given in this case by

and I ∈ RN ×N is an identity matrix. Note that each block of the matrices Am ∈

(1) λ1 (2) λ2 (2)

(1) λ1 (2) λ2 (2)

(1) (1) (1)

(2) (2) (2) (2) (2) (2)

5α2 ∆x2 + 24λ2 ∆t α2 ∆x2 + 12λ2 ∆t

3λ1 ∆t(α1 ∆x2 + 2λ1 ∆t(2 − cos(jπ∆x))) − 4λ21 ∆t2

Thus, we obtain an estimate of the spectral radius of the iteration matrix Σ:

(12α2 ∆x2 (5 − cos(π∆x))) · 0

122 λ2 ∆t(6λ21 ∆t2 − 4λ21 ∆t2 ) λ1

4.1 Numerical Results

xk+1 = S(F(xk ) + k ) + δk . (6.6)

kxk+1 − x∗ k = kS(F(xk ) + k ) + δk − x∗ k = kS(F(xk ) + k ) + δk − S(F(x∗ ))k

and thus in the limit xk+1 → x ,

This is the case if and only if both k and δk converge to zero.

τr kA1 uk1 − b1 (ukΓ )k ≤ τr (kA1 uk1 − b1 (uk−1 k−1 k

Here, we cannot make a statement on the limit of k .