On the convergence rate of the Dirichlet-Neumann iteration for unsteady thermal fluid structure interaction

TL;DR

This paper analyzes the convergence rate of the Dirichlet-Neumann iteration in thermal fluid-structure interaction, providing an exact spectral radius formula and confirming its accuracy through numerical experiments.

Contribution

It derives an exact formula for the spectral radius of the iteration matrix in discretized thermal FSI, linking convergence rate to time step and material properties.

Findings

Convergence rate depends on aspect ratio and material conductivity.

Decreasing time step improves convergence.

The 1D formula estimates 2D and nonlinear cases effectively.

Abstract

We consider the Dirichlet-Neumann iteration for partitioned simulation of thermal fluid-structure interaction, also called conjugate heat transfer. We analyze its convergence rate for two coupled fully discretized 1D linear heat equations with jumps in the material coefficients across these. These are discretized using implicit Euler in time, a finite element method on one domain, a finite volume method on the other one and variable aspect ratio. We provide an exact formula for the spectral radius of the iteration matrix. This shows that for large time steps, the convergence rate is the aspect ratio times the quotient of heat conductivities and that decreasing the time step will improve the convergence rate. Numerical results confirm the analysis and show that the 1D formula is a good estimator in 2D and even for nonlinear thermal FSI applications.

Figures

Tables

Table 1. Table 1: Physical properties of the materials. λ 𝜆 \lambda is the thermal conductivity, ρ 𝜌 \rho the density, c p subscript 𝑐 𝑝 c_{p} the specific heat capacity and α = ρ ​ c p 𝛼 𝜌 subscript 𝑐 𝑝 \alpha=\rho c_{p} .

Material	$\lambda$ (W/mK)	$\rho$ (kg/${\text{m}}^3$)	$c_p$ (J/kgK)	$\alpha$ (J/K ${\text{m}}^3$)
Air	0.0243	1.293	1005	1299.5
Water	0.58	999.7	4192.1	4.1908e6
Steel	48.9	7836	443	3471348

Table 2. Table 2: Temporal and spatial asymptotics of ( 66 ) for the thermal interaction of air at 273 ​ K 273 𝐾 273K with steel at 900 ​ K 900 𝐾 900K , water at 283 ​ K 283 𝐾 283K with steel and air with water.

Case	$\Delta t\to 0$	$\Delta x\to 0$
Air-Steel	0	4.9693e-4 ⋅r
Water-Steel	0	0.0119 ⋅r
Air-Water	0	0.0419 ⋅r

Equations

$${\partial }_t\rho +\nabla \cdot \rho \mathbf{v}$$

$${\partial }_t\rho v_i+\sum _{j=1}^2{\partial }_{x_j}(\rho v_i{v}_j+p{\delta }_{ij})$$

$${\partial }_t\rho E+\nabla \cdot (\rho H{v}_j)$$

$$\rho (\mathbf{x})c_p(\Theta )\frac{d}{dt}\Theta (\mathbf{x},t)=-\nabla \cdot \mathbf{q}(\mathbf{x},t),$$

$${\mathbf{q}}_s(\mathbf{x},t)=-\lambda (\Theta )\nabla \Theta (\mathbf{x},t)$$

$$\lambda (\Theta )=40.1+0.05\Theta -0.0001{\Theta }^2+4.9\cdot 10^{-8}{\Theta }^3$$

$$c_p(\Theta )=-10\ln \left(\frac{e^{-c_{p1}(\Theta )/10}+e^{-c_{p2}(\Theta )/10}}{2}\right)$$

$$c_{p1}(\Theta )=34.2e^{0.0026\Theta }+421.15$$

$$c_{p2}(\Theta )=956.5e^{-0.012(\Theta -900)}+0.45\Theta .$$

$$T(\mathbf{x},t)=\mathbf{\Theta }(\mathbf{x},t),\mathbf{x}\in \Gamma ,$$

$${\mathbf{q}}_f(\mathbf{x},t)\cdot \mathbf{n}(x)={\mathbf{q}}_s(\mathbf{x},t)\cdot \mathbf{n}(x),\mathbf{x}\in \Gamma .$$

$$\frac{d}{dt}\mathbf{u}+\mathbf{h}(\mathbf{u},{\mathbf{\Theta }}_{\Gamma })=0,$$

$$\mathbf{M}(\mathbf{\Theta })\frac{d}{dt}\mathbf{\Theta }+\mathbf{A}(\mathbf{\Theta })\mathbf{\Theta }={\mathbf{q}}_b^f+{\mathbf{q}}_b^{\Gamma }(\mathbf{u}).$$

$${\mathbf{u}}^{n+1}-{\mathbf{u}}^n+\Delta t\mathbf{h}({\mathbf{u}}^{n+1},{\mathbf{\Theta }}_{\Gamma }^{n+1})=0,$$

$$\mathbf{M}({\mathbf{\Theta }}^{n+1})({\mathbf{\Theta }}^{n+1}-{\mathbf{\Theta }}^n)+\Delta t\mathbf{A}({\mathbf{\Theta }}^{n+1}){\mathbf{\Theta }}^{n+1}=\Delta t({\mathbf{q}}_b^f+{\mathbf{q}}_b^{\Gamma }({\mathbf{u}}^{n+1})).$$

$${\mathbf{u}}^{n+1,k+1}-{\mathbf{u}}^n+\Delta t\mathbf{h}({\mathbf{u}}^{n+1,k+1},{\mathbf{\Theta }}_{\Gamma }^{n+1,k})=0,$$

$$\mathbf{M}({\mathbf{\Theta }}^{n+1,k+1})({\mathbf{\Theta }}^{n+1,k+1}-{\mathbf{\Theta }}^n)+\Delta t\mathbf{A}({\mathbf{\Theta }}^{n+1,k+1}){\mathbf{\Theta }}^{n+1,k+1}=\Delta t({\mathbf{q}}_b^f+{\mathbf{q}}_b^{\Gamma }({\mathbf{u}}^{n+1,k+1})),$$

$$\Vert {\mathbf{\Theta }}_{\Gamma }^{k+1}-{\mathbf{\Theta }}_{\Gamma }^k\Vert \le \tau$$

$$\begin{array}{rl}{\alpha }_m\frac{\partial u_m(\text{x},t)}{\partial t}-\nabla \cdot ({\lambda }_m\nabla u_m(\text{x},t))& =0,t\in [t_0,t_f],\text{x}\in {\Omega }_m\subset {\mathbb{R}}^d,m=1,2,\\ u_m(\text{x},t)& =0,t\in [t_0,t_f],\text{x}\in \partial {\Omega }_m\backslash \Gamma ,\\ u_1(\text{x},t)& =u_2(\text{x},t),\text{x}\in \Gamma ,\\ {\lambda }_2\frac{\partial u_2(\text{x},t)}{\partial {\text{n}}_2}& =-{\lambda }_1\frac{\partial u_1(\text{x},t)}{\partial {\text{n}}_1},\text{x}\in \Gamma ,\\ u_m(\text{x},0)& =u_m^0(\text{x}),\text{x}\in {\Omega }_m,\end{array}$$

$$D_m=\frac{{\lambda }_m}{{\alpha }_m},\text{mbox}with{\alpha }_m={\rho }_m{c}_{p_m}$$

$$\beta =\left|-\frac{{\lambda }_1}{{\lambda }_2}\sqrt{\frac{D_2}{D_1}}\frac{\tanh \left(-\frac{1}{\sqrt{D_2\Delta t}}\right)}{\tanh \left(\frac{1}{\sqrt{D_1\Delta t}}\right)}\right|.$$

$$\beta \approx \frac{{\lambda }_1}{{\lambda }_2}\sqrt{\frac{D_2}{D_1}}.$$

$$\beta \approx \frac{{\lambda }_1}{{\lambda }_2}\sqrt{\frac{D_2}{D_1}}\frac{\sqrt{D_1\Delta t}}{\sqrt{D_2\Delta t}}=\frac{{\lambda }_1}{{\lambda }_2}.$$

$${\text{M}}_1{\dot{\text{u}}}_I^{(1)}+{\text{M}}_{I\Gamma }^{(1)}{\dot{\text{u}}}_{\Gamma }+{\text{A}}_1{\text{u}}_I^{(1)}+{\text{A}}_{I\Gamma }^{(1)}{\text{u}}_{\Gamma }=0.$$

$${\text{M}}_2{\dot{\text{u}}}_I^{(2)}+{\text{M}}_{I\Gamma }^{(2)}{\dot{\text{u}}}_{\Gamma }+{\text{A}}_2{\text{u}}_I^{(2)}+{\text{A}}_{I\Gamma }^{(2)}{\text{u}}_{\Gamma }=0.$$

$$-{\lambda }_1\frac{\partial u_1}{\partial {\text{n}}_1}\approx \frac{{\lambda }_1}{2\Delta x}(4u_{1,N}(t)-u_{1,N-1}(t)-3u_{\Gamma }).$$

$$\begin{array}{rl}& {\lambda }_2{\int }_{\Gamma }\frac{\partial u_2}{\partial {\text{n}}_2}{\varphi }_{j}dS={\lambda }_2{\int }_{{\Omega }_2}(\Delta u_2{\varphi }_j+\nabla u_2\nabla {\varphi }_j)d\text{x}\\ & ={\alpha }_2{\int }_{{\Omega }_2}\frac{d}{dt}{u}_2{\varphi }_j+{\lambda }_2{\int }_{{\Omega }_2}\nabla u_2\nabla {\varphi }_{j}d\text{x}.\end{array}$$

$${\text{M}}_{\Gamma \Gamma }^{(2)}{\dot{\text{u}}}_{\Gamma }+{\text{M}}_{\Gamma I}^{(2)}{\dot{\text{u}}}_I^{(2)}+{\text{A}}_{\Gamma \Gamma }^{(2)}{\text{u}}_{\Gamma }+{\text{A}}_{\Gamma I}^{(2)}{\text{u}}_I^{(2)}=-{\text{M}}_{\Gamma \Gamma }^{(1)}{\dot{\text{u}}}_{\Gamma }-{\text{M}}_{\Gamma I}^{(1)}{\dot{\text{u}}}_I^{(1)}-{\text{A}}_{\Gamma \Gamma }^{(1)}{\text{u}}_{\Gamma }-{\text{A}}_{\Gamma I}^{(1)}{\text{u}}_I^{(1)},$$

$$\tilde{\text{M}}\dot{\text{u}}+\tilde{\text{A}}\text{u}=0,$$

$$\displaystyle\begin{split}\tilde{\textbf{M}}=\left(\begin{array}[]{ccc}\textbf{M}_{1}&\mathbf{0}&\textbf{M}_{I\Gamma}^{(1)}\\ \mathbf{0}&\textbf{M}_{2}&\textbf{M}_{I\Gamma}^{(2)}\\ \textbf{M}_{\Gamma I}^{(1)}&\textbf{M}_{\Gamma I}^{(2)}&\textbf{M}_{\Gamma\Gamma}^{(1)}+\textbf{M}_{\Gamma\Gamma}^{(2)}\end{array}\right),\ \ \tilde{\textbf{A}}=\left(\begin{array}[]{ccc}\textbf{A}_{1}&\mathbf{0}&\textbf{A}_{I\Gamma}^{(1)}\\ \mathbf{0}&\textbf{A}_{2}&\textbf{A}_{I\Gamma}^{(2)}\\ \textbf{A}_{\Gamma I}^{(1)}&\textbf{A}_{\Gamma I}^{(2)}&\textbf{A}_{\Gamma\Gamma}^{(1)}+\textbf{A}_{\Gamma\Gamma}^{(2)}\end{array}\right).\end{split}$$

$$\text{A}{\text{u}}^{n+1}=\tilde{\text{M}}{\text{u}}^n,$$

Full text

Keywords: Thermal Fluid Structure Interaction, Coupled Problems, Transmission Problem, Fixed Point Iteration, Dirichlet-Neumann Iteration

1 Introduction

The Dirichlet-Neumann iteration is a basic method in both domain decomposition and fluid structure interaction (FSI). In the latter case, the iteration arises in a partitioned approach [11], where different codes for the sub-problems are reused and the coupling is done by a master program which calls interface functions of the other codes. This allows to reuse existing software for each sub-problem, in contrast to a monolithic approach, where a new code is tailored for the coupled equations. To satisfy coupling conditions at the interface, the subsolvers are iterated by providing Dirichlet- and Neumann data for the other solver in a sequential manner, giving rise to its name.

In the domain decomposition context, the iteration has two main problems, namely slow convergence and the need for an implementation using a red-black colouring. The slow convergence can be slightly improved using a relaxation procedure. In fluid structure interaction, there are typically only two domains, coupled along an interface, making the application straight forward. The convergence rate for the interaction of a flexible structure with a fluid has been analyzed in [30]. There, the added mass effect is proven to be dependent on the step size for compressible flows and independent for incompressible flows. However, the convergence rate is not great for the coupling between a compressible fluid and a structure [10], which is why a lot of effort goes into convergence acceleration. Furthermore, for incompressible fluids it is known that the ratio of densities of the materials plays an important role [1, 9]. Finally, the Dirichlet-Neumann iteration was reported to be a very fast solver for thermal fluid structure interaction [4].

Our prime motivation here is thermal interaction between fluids and structures, also called conjugate heat transfer. There are two domains with jumps in the material coefficients across the connecting interface. Conjugate heat transfer plays an important role in many applications and its simulation has proved essential [2]. Examples for thermal fluid structure interaction are cooling of gas-turbine blades, thermal anti-icing systems of airplanes [8], supersonic reentry of vehicles from space [21, 17], gas quenching, which is an industrial heat treatment of metal workpieces [15, 28] or the cooling of rocket nozzles [18, 19].

For the case of coupled heat equations, a 1D stability analysis was presented by Giles [14]. There, an explicit time integration method was chosen with respect to the interface unknows. On the other hand, Henshaw and Chand provided in [16] a method to analyze stability and convergence speed of the Dirichlet-Neumann iteration in 2D based on applying the continuous Fourier transform to the semi-discretized equations. They show that the ratios of thermal conductivities and diffusivities of the materials play an important role. This is similar to the result names above in classical FSI with incompressible fluids where the performance is affected by the ratio of densities of the materials [1, 9].

However, in the fully discrete case we observe in some cases that the iteration behaves differently, because some aspects of the problem are not taken into account by the semidiscrete analysis: The effect of $\Delta t$ is not accurately represented and neither are possibly different mesh widths in the two problems. This matters particularly for compressible fluids where a high aspect ratio grid is needed to accurately represent the boundary layer. This leads to geometric stiffness that significantly influences the convergence rate, as we will show.

For the fully discrete case, the convergence rate is in principle analyzed in any standard book on domain decomposition methods, e.g. [26, 29]. There, the iteration matrix is derived in terms of the stiffness and mass matrices of finite element discretizations and the convergence rate is the spectral radius of that. However, this does not provide a quantitative answer, since the spectral radius is unknown. Computing the spectral radius is in general a non trivial task. In our context, the material properties are discontinuous across the interface and as a consequence, computing the spectral radius of the iteration matrix is even more difficult.

In [23, 25], a convergence analysis of the Dirichlet-Neumann iteration for the unsteady transmission problem using finite element methods (FEM) on both subdomains is presented. A similar analysis using finite differences (FDM) on one domain and FEM on the other one can be found in [24]. In addition, the corresponding analysis when coupling finite volumes (FVM) with FEM is described in [5, 23]. All these results assume equal mesh sizes on both subdomains, i.e, the aspect ratio is equal to one.

Thus, we present here a complete discretization of the coupled problem using FVM in space on one domain and FEM on the other one with variable aspect ratio $r$. We consider this to be a relevant case, because these are the standard discretizations for the subproblems. The implicit Euler method is used for the temporal discretization. Then, we derive the spectral radius of the iteration matrix exactly in terms of the eigendecomposition of the resulting matrices for the one dimensional case. The asymptotic convergence rates when approaching the continuous case in either time or space are also determined. In the spatial limit, the convergence rate turns out to be proportional to the aspect ratio $r$, whereas in the temporal limit, we obtain 0. Note that for FEM-FEM couplings, this is not the case. Moreover, we also include numerical results where it is shown that the one dimensional formula is a good estimator for a 2D version of the coupled heat equations and for two non linear FSI models, namely the cooling of a flat plate and the cooling of a flanged shaft.

An outline of the paper now follows. In section 2, we describe the model and discretization, as well as the coupling conditions and the Dirichlet-Neumann iteration. Two thermal FSI test cases are introduced in section 3: the cooling of a flat plate and of a flanged shaft. For these, we present numerical convergence rates, motivating further analysis. A model problem, consisting of two coupled discretized heat equations, is presented in section 4 and then analyzed in 1D in section 5. In section 6, extension of the analysis to 2D and different discretizations are discussed. In section 7, the analytical results are compared to linear and nonlinear numerical results.

2 Thermal FSI Methodology

The basic setting we are in is that on a domain $\Omega_{1}\subset\mathbb{R}^{d}$ where $d$ corresponds to the spatial dimension, the physics is described by a fluid model, whereas on a domain $\Omega_{2}\subset\mathbb{R}^{d}$, 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 interface $\Gamma\subset\partial\Omega_{1}\cup\partial\Omega_{2}$. Note that $\Gamma$ might be a true subset of the intersection, because the structure could be insulated. At the interface $\Gamma$, coupling conditions are prescribed that model the interaction between fluid and structure. For the thermal coupling problem, 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 compressible Navier-Stokes equations, which are a second order system of conservation laws (mass, momentum, energy) modeling compressible flow. We consider the two dimensional case, written in conservative variables density $\rho$, momentum ${\bf m}=\rho{\bf v}$ and energy per unit volume $\rho E$ as:

[TABLE]

Here, enthalpy is given by $H=E+p/\rho$ with $p=(\gamma-1)\rho(E-1/2|v|^{2})$ being the pressure and $\gamma=1.4$ the adiabatic index for an ideal gas. Furthermore, ${\bf q}_{f}=(q_{1},q_{2})^{T}$ represents the heat flux and $\textbf{S}=(S_{ij})_{i,j=1,2}$ the viscous shear stress tensor. As the equations are dimensionless, the Reynolds number $Re$ and the Prandtl number $Pr$ appear. The system is closed by the equation of state for the pressure $p=(\gamma-1)\rho 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 $\Omega_{1}$ except for $\Gamma$, where we have the coupling conditions. In the Dirichlet-Neumann coupling, a temperature value is enforced at $\Gamma$.

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 $\Theta$

[TABLE]

where

[TABLE]

denotes the heat flux vector. For alloys, the specific heat capacity $c_{p}$ and heat conductivity $\lambda$ are temperature-dependent and highly nonlinear.

As an example, an empirical model for the steel 51CrV4 was suggested in [27]. This was obtained measurements and a least squares fit to a chosen curve. The coefficient functions are then

[TABLE]

and

[TABLE]

with

[TABLE]

and

[TABLE]

For the mass density one has $\rho=$7836\text{,}\mathrm{kg}\text{,}{\mathrm{m}}^{-3}$$.

Finally, on the boundary, we have Neumann conditions $\mathbf{q}_{s}(\mathbf{x},t)\cdot\mathbf{n}(\mathbf{x})=q_{b}(\mathbf{x},t)$.

2.3 Coupling Conditions

As mentioned before, the coupling conditions are that temperature and the normal component of the heat flux are continuous across the interface, i.e;

[TABLE]

where $T$ is the fluid temperature and $\mathbf{\Theta}$ the structure temperature and

[TABLE]

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 the following equation for all unknowns on $\Omega_{1}$:

[TABLE]

where ${\bf h}({\bf u},\boldsymbol{\Theta}_{\Gamma})$ represents the spatial discretization and its dependence on the temperatures on the discrete interface to the structure, here denoted by $\boldsymbol{\Theta}_{\Gamma}$.

Regarding structural mechanics, the use of finite element methods is ubiquitious. Therefore, we will also follow that approach here, using quadratic finite element one gets the following nonlinear equation for all unknowns on $\Omega_{2}$:

[TABLE]

Here, $\mathbf{M}$ is the mass matrix, also called heat capacity matrix for this problem and $\mathbf{A}$ is the heat conductivity and stiffness matrix. The vector ${\bf\Theta}$ consists of all discrete temperature unknowns and ${\bf q}_{b}^{\Gamma}({\bf{\bf u}})$ is the discrete heat flux vector on the coupling interface to the fluid, whereas ${\bf q}_{b}^{f}$ corresponds to boundary heat fluxes independent of the fluid, for example at insulated boundaries.

2.5 Time Discretization

In time, we use the implicit Euler method with constant time step $\Delta t$. For the system (9)-(10) we obtain

[TABLE]

[TABLE]

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 (11)-(12). This corresponds to alternately solving equation (11) on $\Omega_{1}$ with Dirichlet data on $\Gamma$ and (12) on $\Omega_{2}$ with Neumann data on $\Gamma$.

Thus, one gets for the $k$-th iteration the two decoupled equation systems

[TABLE]

[TABLE]

with some initial condition $\mathbf{\Theta}_{\Gamma}^{0}$. The iteration is terminated according to the standard criterion

[TABLE]

where $\tau$ is a user defined tolerance.

3 Thermal FSI Test Cases

In this section we present two thermal FSI test cases that are solved using the methodology explained in the previous section. The aim of this paper is to estimate the convergence rates of the Dirichlet-Neumann iteration used as a solver for thermal FSI problems. Therefore, we first want to illustrate the behavior for two examples before proceeding to the convergence analysis in the next section. Two different test cases are discussed: the cooling of a flat plate and the cooling of a flanged shaft. For the first problem, structured grids are used and for the second, unstructured grids.

For the coupling, the Dirichlet-Neumann method as presented in (13)-(14) is used. A fixed tolerance of $1e-8$ is chosen for all involved equation solvers. The coupling code used has been developed in a series of papers [7, 6, 4]. It’s main feature is time adaptivity, which is not employed here. The coupling between the solvers is done using the Component Template Library (CTL) [20]. In the fluid, the DLR TAU-Code in its 2014.2 version is employed [13], which is a cell-vertex-type finite volume method with AUSMDV as flux function and a linear reconstruction to increase the order of accuracy. The finite element code uses quadratic finite elements and is the inhouse code Native of the Institute for Static and Dynamic at the University of Kassel.

3.1 Flow over a plate

The first test case is the cooling of a flat steel plate resembling a simple work piece [7]. The work piece is initially at a much higher temperature than the fluid and then cooled by a constant laminar air stream, see figure 1.

The inlet is given on left, where air enters the domain with an initial velocity of $\mbox{Ma}_{\infty}=0.8$ in horizontal direction and a temperature of 273K. Regarding the initial condition in the structure, a constant temperature of 900K at $t=0$ is chosen throughout.

The grid, see figure 2, is chosen cartesian and equidistant in the structural part. In the fluid region the thinest cells are on the boundary and then become coarser in $y$-direction with a maximal aspect ratio of $r=1.7780e5$. The points of the primary fluid grid and the nodes of the structural grid match on the interface $\Gamma$ and there are 9660 cells in the fluid region and $n_{x}\times n_{y}=120\times 9=1080$ elements with $121\times 10=1210$ nodes in the region of the structure.

Figure 5a shows the convergence behaviour of the Dirichlet-Neumann iteration against the time step $\Delta t$. One observes how the convergence rates is roughly proportional to the time step $\Delta t$. Furthermore, even for $\Delta t=1$ a reduction of the error by a factor of ten per iteration is achieved.

3.2 Cooling of a flanged shaft

The second test case is the cooling of a flanged steel shaft by cold high pressured air (this process is also known as gas quenching) [31]. Here, we have a hot flanged shaft that is cooled by cold high pressured air coming out of small tubes, see figure 3. We assume symmetry along the horizontal axis in order to consider one half of the flanged shaft and two tubes blowing air at it. We also assume that the air leaves the tube in straight and uniform way at a Mach number of 1.2. Moreover, we assume a freestream in $x$-direction of Mach 0.005. The Reynolds number is $Re=2500$ and the Prandtl number $Pr=0.72$.

The grid, see figure 4, consists of 279212 cells in the fluid, which is the dual grid of an unstructured grid of quadrilaterals in the boundary layer and triangles in the rest of the domain, and 1997 quadrilateral elements in the structure. Regarding the initial conditions, we use the procedure explained in [4].

Figure 5b shows the convergence behaviour of the Dirichlet-Neumann iteration against the time step $\Delta t$. The convergence rate is again about proportional to the time step size and again convergent even for very large time steps. If we compare the rates for the two problems, we observe that for a given $\Delta t$, the iteration is about a factor ten faster for the plate.

Summarizing, the Dirichlet-Neumann iteration is a very fast solver for thermal FSI. To understand this better, we perform in the next section a convergence analysis for the case of two coupled linear heat equations.

4 A Model Problem: Coupled Heat Equations

We present here a convergence analysis of the unsteady transmission problem with mixed discretizations. In particular, we choose a finite volume method (FVM) on the first subdomain and a finite element method (FEM) on the second subdomain.

4.1 Model Problem

The unsteady transmission problem is as follows, where we consider a domain $\Omega\subset\mathbb{R}^{d}$ which is cut into two subdomains $\Omega_{1}\cup\Omega_{2}=\Omega$ with transmission conditions at the interface $\Gamma=\Omega_{1}\cap\Omega_{2}$:

[TABLE]

where $\textbf{n}_{m}$ is the outward normal to $\Omega_{m}$ for $m=1,2$.

The constants $\lambda_{1}$ and $\lambda_{2}$ describe the thermal conductivities of the materials on $\Omega_{1}$ and $\Omega_{2}$ respectively. $D_{1}$ and $D_{2}$ represent the thermal diffusivities of the materials and they are defined by

[TABLE]

where $\rho_{m}$ represents the density and $c_{p_{m}}$ the specific heat capacity of the material placed in $\Omega_{m}$, $m=1,2$.

We always use the implicit Euler method for time discretization. With regards to the spatial discretization, we use FVM on $\Omega_{1}$ and FEM on $\Omega_{2}$.

4.2 Semidiscrete Analysis

Before we present in the next section an analysis for the fully discrete equations, we want to describe previous results about the behaviour of the Dirichlet-Neumann iteration for the transmission problem in the semi discrete case.

Henshaw and Chand applied in [16] the implicit Euler method for the time discretization on both equations in (16) but kept the space continuous. Then, they applied the Fourier transform in space in order to transform the second order derivatives into algebraic expressions. Once they have a coupled system of algebraic equations, they insert one into the other and obtain the Dirichlet-Neumann convergence rate $\beta$:

[TABLE]

For $\Delta t$ small enough, we have $\tanh\left(-1/\sqrt{D_{2}\Delta t}\right)\approx-1$ and $\tanh\left(1/\sqrt{D_{1}\Delta t}\right)\approx 1$ and therefore:

[TABLE]

On the other hand, for $\Delta t$ big enough, we have $\tanh\left(-1/\sqrt{D_{2}\Delta t}\right)\approx-1/\sqrt{D_{2}\Delta t}$ and $\tanh\left(1/\sqrt{D_{1}\Delta t}\right)\approx 1/\sqrt{D_{1}\Delta t}$ and therefore:

[TABLE]

Figure 6 shows $\beta$ as a function of $\Delta t$. It is almost constant, except for a short dynamic transition between $(\lambda_{1}/\lambda_{2})\sqrt{D_{2}/D_{1}}$ and $\lambda_{1}/\lambda_{2}$.

Finally, one observes in (20) that the convergence rates of the Dirichlet-Neumann iteration are given by the quotient of thermal conductivities for $\Delta t$ large. This suggests that strong jumps in the thermal conductivities of the materials give fast convergence.

4.3 Space Discretization

We now describe a rather general space discretization of the model problem. The core property we need is that the meshes of $\Omega_{1}$ and $\Omega_{2}$ share the same nodes on $\Gamma$ as shown in Figure 7. Furthermore, we need that there is a specific set of unknowns associated with the interface nodes. Otherwise, we allow at this point for arbitraty meshes on both sides.

Then, letting $\textbf{u}_{I}^{(1)}$ correspond to the unknowns on $\Omega_{1}$ and $\textbf{u}_{\Gamma}$ to the unknowns at the interface $\Gamma$, we can write a general discretization of the first equation in (16) in a compact form as:

[TABLE]

On the other hand, a general discretization of the first equation in (16) on $\Omega_{2}$ can be written as:

[TABLE]

where $\textbf{u}_{I}^{(2)}$ correspond to the unknowns on $\Omega_{2}$.

To close the system, we need an approximation of the normal derivatives at $\Gamma$. For the FVM on $\Omega_{1}$, we approximate the normal derivative with respect to $u_{1}$ using second order one-sided finite differences:

[TABLE]

On the other hand, let $\phi_{j}$ be a nodal FE basis function on $\Omega_{2}$ for a node on $\Gamma$ we observe that the normal derivative with respect to $u_{2}$ can be written as a linear functional using Green’s formula [29, pp. 3]. Thus, the approximation of the normal derivative is given by

[TABLE]

Consequently, the equation

[TABLE]

is a discrete version of the fourth equation in (16) and completes the system (21)-(22). Notice that the left hand side of (25) comes from (LABEL:approxnormalderfem) and the right hand side from (23). We can now write the coupled equations (21), (22) and (25) as an ODE for the vector of unknowns $\textbf{u}=\left(\textbf{u}_{I}^{(1)},\textbf{u}_{I}^{(2)},\textbf{u}_{\Gamma}\right)^{T}$

[TABLE]

where

[TABLE]

4.4 Time Discretization

Applying the implicit Euler method with time step $\Delta t$ to the system (25), we get for the vector of unknowns $\textbf{u}^{n+1}=(\textbf{u}_{I}^{(1),n+1},\textbf{u}_{I}^{(2),n+1},\textbf{u}_{\Gamma}^{n+1})^{T}$

[TABLE]

where

[TABLE]

with $\textbf{M}_{\Gamma\Gamma}=\textbf{M}_{\Gamma\Gamma}^{(1)}+\textbf{M}_{\Gamma\Gamma}^{(2)}$ and $\textbf{A}_{\Gamma\Gamma}=\textbf{A}_{\Gamma\Gamma}^{(1)}+\textbf{A}_{\Gamma\Gamma}^{(2)}$.

4.5 Dirichlet-Neumann Iteration

We now employ a Dirichlet-Neumann iteration to solve the discrete system (27). This corresponds to alternately solving the discretized equations of the transmission problem (16) on $\Omega_{1}$ with Dirichlet data on $\Gamma$ and the discretization of (16) on $\Omega_{2}$ with Neumann data on $\Gamma$.

Therefore, from (27) one gets for the $k$-th iteration the two equation systems

[TABLE]

[TABLE]

to be solved in succession. Here,

[TABLE]

and

[TABLE]

[TABLE]

with some initial condition, here $\textbf{u}_{\Gamma}^{n+1,0}=\textbf{u}_{\Gamma}^{n}$. The iteration is terminated according to the standard criterion $\|\textbf{u}_{\Gamma}^{k+1}-\textbf{u}_{\Gamma}^{k}\|\leq\tau$ where $\tau$ is a user defined tolerance [3].

One way to analyze this method is to write it as a splitting method for (27) 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 (32)-(33) as an iteration for $\textbf{u}_{\Gamma}^{n+1}$ to restrict the size of the space to the dimension of $\textbf{u}_{\Gamma}$ which is much smaller. To this end, we isolate the term $\textbf{u}_{I}^{(1),n+1,k+1}$ in (32) and $\textbf{u}_{I}^{(2),n+1,k+1}$ in the first equation in (33) and we insert the resulting expressions into the second equation in (33). Consequently, the iteration $\textbf{u}_{\Gamma}^{n+1,k+1}=\Sigma\textbf{u}_{\Gamma}^{n+1,k}+\psi^{n}$ is obtained with iteration matrix

[TABLE]

where

[TABLE]

for $m=1,2$ and $\psi^{n}$ contains terms that depend only on the solutions at the previous time step. Notice that $\Sigma$ is a discrete version of the Steklov-Poincaré operator.

Thus, the Dirichlet-Neumann iteration is a linear iteration and the rate of convergence is described by the spectral radius of the iteration matrix $\Sigma$.

5 One-Dimensional Convergence Analysis

The derivation so far was for a rather general discretization. In this section, we study the iteration matrix $\Sigma$ for a specific FVM-FEM discretization in 1D. We will give an exact formula for the convergence rates. The behaviour of the rates when approaching both the continuous case in time and space is also given.

Specifically, we use $\Omega_{1}=[-1,0]$, $\Omega_{2}=[0,1]$. For the FVM discretization, we consider a primal grid, i.e, we discretize $\Omega_{1}$ into $N_{1}$ equal sized grid cells of size $\Delta x_{1}=1/(N_{1}+1)$, and define $x_{i}=i\Delta x_{1}$, so that $x_{i}$ is the center of the cell $i$, see figure 8. The edges of cell $i$ are then $x_{i-1/2}$ and $x_{i+1/2}$ and they form the corresponding dual grid. Moreover, we use the flux function

[TABLE]

to approximate the flux, which results in a second order scheme. For the FEM discretization, we use the standard piecewise-linear polynomials as test functions. Here we discretize $\Omega_{2}$ into $N_{2}$ equal sized cells of size $\Delta x_{2}=1/(N_{2}+1)$.

For the coupling between a compressible fluid and a structure, there would be a boundary layer in the fluid, meaning that the mesh would be very fine in direction normal to the boundary, implying $\Delta x_{1}\ll\Delta x_{2}$.

With \textbf{e}_{m,j}=\left(\begin{array}[]{ccccccc}0&\cdots&0&1&0&\cdots&0\end{array}\right)^{T}\in\mathbb{R}^{N_{m}} where the only nonzero entry is located at the $j$-th position, the discretization matrices are given by

[TABLE]

[TABLE]

[TABLE]

where $\textbf{A}_{m}$, $\textbf{M}_{m}\in\mathbb{R}^{N_{m}\times N_{m}}$, $\textbf{A}_{I\Gamma}^{(m)}$, $\textbf{M}_{I\Gamma}^{(2)}\in\mathbb{R}^{N_{m}\times 1}$ and $\textbf{A}_{\Gamma I}^{(m)}$, $\textbf{M}_{\Gamma I}^{(2)}\in\mathbb{R}^{1\times N_{m}}$ for $m=1,2$.

In this case, $\textbf{M}_{1}=\alpha_{1}\textbf{I}$, $\textbf{M}_{I\Gamma}^{(1)}=\textbf{M}_{\Gamma\Gamma}^{(1)}=\textbf{M}_{\Gamma I}^{(1)}=\mathbf{0}$. Thus,

[TABLE]

[TABLE]

Note that the iteration matrix $\Sigma$ is just a real number in this case and thus its spectral radius is its modulus. One computes $\textbf{S}^{(1)}$ and $\textbf{S}^{(2)}$ by inserting the corresponding matrices specified above in (48) and (49) obtaining

[TABLE]

[TABLE]

where $\alpha_{ij}^{1}$ represents the entries of the matrix $(\alpha_{1}\textbf{I}-\Delta t\textbf{A}_{1})^{-1}$ and $\alpha_{ij}^{2}$ the entries of $(\textbf{M}_{2}+\Delta t\textbf{A}_{2})^{-1}$ for $i,j=1,...,N_{1}$ and $i,j=1,...,N_{2}$ respectively. Observe that the matrices $(\alpha_{1}\textbf{I}-\Delta t\textbf{A}_{1})^{-1}$ and $(\textbf{M}_{2}+\Delta t\textbf{A}_{2})$ are tridiagonal Toeplitz matrices but their inverses are full matrices. The computation of the exact inverses is based on a recursive formula which runs over the entries [12] and consequently, it is non trivial to compute $\alpha_{N_{1}N_{1}}^{1}$, $\alpha_{N_{1}-1N_{1}}^{1}$ and $\alpha_{11}^{2}$ this way.

Due to these difficulties, we propose to rewrite the matrices $(\alpha_{1}\textbf{I}-\Delta t\textbf{A}_{1})^{-1}$ and $(\textbf{M}_{2}+\Delta t\textbf{A}_{2})^{-1}$ in terms of their eigendecomposition:

[TABLE]

[TABLE]

where the matrix $\textbf{V}_{N}$ has the eigenvectors of any symmetric tridiagonal Toeplitz matrix of dimension $N$ as columns. The entries of $\textbf{V}_{N_{1}}$ and $\textbf{V}_{N_{2}}$ are not dependent on the entries of $\alpha_{1}\textbf{I}-\Delta t\textbf{A}_{1}$ or $\textbf{M}_{2}+\Delta t\textbf{A}_{2}$ due to their symmetry. Moreover, the matrices $\Lambda_{1}$ and $\Lambda_{2}$ are diagonal matrices having the eigenvalues of $\alpha_{1}\textbf{I}-\Delta t\textbf{A}_{1}$ or $\textbf{M}_{2}+\Delta t\textbf{A}_{2}$ as entries respectively. These are known and given e.g. in [22, pp. 514-516]:

[TABLE]

The entries $\alpha_{N_{1}N_{1}}^{1}$, $\alpha_{N_{1}-1N_{1}}^{1}$ and $\alpha_{11}^{2}$ of the matrices $(\alpha_{1}\textbf{I}-\Delta t\textbf{A}_{1})^{-1}$ and $(\textbf{M}_{2}+\Delta t\textbf{A}_{2})^{-1}$, respectively, are now computed through their eigendecomposition resulting in

[TABLE]

[TABLE]

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

Now, inserting (55), (56) and (57) into (50) and (51) we get for $\textbf{S}^{(1)}$ and $\textbf{S}^{(2)}$:

[TABLE]

[TABLE]

With this we obtain an explicit formula for the spectral radius of the iteration matrix $\Sigma$ as a function of $\Delta x_{1}$, $\Delta x_{2}$ and $\Delta t$:

[TABLE]

To simplify this, the finite sums $\sum_{i=1}^{N_{1}}\sin^{2}(i\pi\Delta x_{1})$ and $\sum_{i=1}^{N_{2}}\sin^{2}(i\pi\Delta x_{2})$ can be computed. We first rewrite the sum of squared sinus terms into a sum of cosinus terms using the identity $\sin^{2}(x/2)=(1-\cos(x))/2$. Then, the resulting sum can be converted into a geometric sum using Euler’s formula. We thus obtain after some calculations:

[TABLE]

[TABLE]

Inserting (64) and (65) into (63) we get after some manipulations

[TABLE]

We could not find a way of simplifying the finite sum (60) because $\Delta x_{2}$ depends on $N_{2}$ (i.e., $\Delta x_{2}=1/(N_{2}+1)$). However, (66) is a computable formula that gives exactly the convergence rates of the Dirichlet-Neumann iteration for given $\Delta t$, $\Delta x_{m}$, $\alpha_{m}$ and $\lambda_{m}$, $m=1,2$.

We are now interested in the asymptotics of (66) for $\Delta t\rightarrow 0$ and $\Delta x_{1}\rightarrow 0$ with $\Delta x_{2}=r\cdot\Delta x_{1}$ where $r:=\Delta x_{2}/\Delta x_{1}$ is a fixed aspect ratio. This is motivated by the assumption that the resolution in the fluid in direction tangential to the wall would be similar to the resolution in the structure. We obtain:

[TABLE]

[TABLE]

To simplify (68), it is well known that the finite sums $\sum_{i=1}^{N_{1}}\cos(i\pi\Delta x_{1})$, $\sum_{i=1}^{N_{2}}\cos(i\pi r\Delta x_{1})$ and $\sum_{i=1}^{N_{1}}\cos^{2}(i\pi\Delta x_{1})$ can be computed by using Euler’s formula to convert them into geometric sums. We thus obtain after some calculations:

[TABLE]

In order to compute the third sum, we rewrite the sum of squared cosinus terms into a sum of sinus terms using the identity $\cos^{2}(x/2)=(1+\cos(x))/2$ and the apply the same technique:

[TABLE]

Inserting (69) and (70) into (68) we get

[TABLE]

From the result obtained in (67) we can conclude that the convergence rate goes to zero when the time step decreases and therefore, the iteration will be fast for $\Delta t$ small and can always be made to converge by decreasing $\Delta t$. This is consistent with the behavior of the cooling of a flat plat and the flanged shaft presented earlier in figures 5a and 5b.

On the other hand, from the spatial asymptotics (LABEL:finallimitdeltax1D) we can observe that strong jumps in the thermal conductivities of the materials placed in $\Omega_{1}$ and $\Omega_{2}$ will imply fast convergence. This is often the case when modelling thermal fluid structure interaction, since fluids typically have lower thermal conductivities than structures.

Finally, the aspect ratio $r$ also influences the behavior of the fixed point iteration, i.e, the rates will become smaller the higher the aspect ratio, e.g. the higher the Reynolds number in the fluid. This phenomenon is not unknown for PDE discretizations and is referred to as geometric stiffness. As is the case here, refining the mesh to reduce the aspect ratio would lead to faster convergence of the iterative method.

Before presenting numerical results we want to show the results obtained for different space discretization combinations with the same constant mesh width on both subdomains.

6 Extension of the Analysis

In this section we want to extend the results presented in the previous section by reviewing similar analysis for other choices of space discretizations. In particular, FEM-FEM coupling and 2D FVM-FEM with $r=1$.

Firstly, when one uses a linear FEM discretization in 1D and the same mesh width on both subdomains (i.e, $r=1$) and applies the same analysis as in the previous section, the corresponding limits for the spectral radius of the iteration matrix $\Sigma$ are given by [23, 25]:

[TABLE]

[TABLE]

When we compare these with the asymptotics obtained with FVM-FEM discretizations (67)-(LABEL:finallimitdeltax1D), we observe that while the spatial limit is the same, the temporal limit does not match. This arises from differences in the matrix $\textbf{S}^{(1)}$ in (46). In the FEM-FEM context, the matrices $\textbf{S}^{(1)}$ and $\textbf{S}^{(2)}$ lead to the same expression with only different material coefficients ($\alpha_{1}$, $\alpha_{2}$, $\lambda_{1}$, $\lambda_{2}$). Because of this, the limits of $\rho(\Sigma)$ are quotients of those coefficients. However, the situation is different in the FVM-FEM context. There, the matrix $\textbf{S}^{(1)}$ in (48) is missing several mass matrices if we compare it with $\textbf{S}^{(2)}$ in (49). This unsymmetry between $\textbf{S}^{(1)}$ and $\textbf{S}^{(2)}$ causes that the limit of $\rho(\Sigma)$ when $\Delta t\rightarrow 0$ is not balanced between the numerator and the denominator, resulting in 0.

This implies that opposed to the FVM-FEM case, where convergence can always be achieved by decreasing the time step, that for an FEM-FEM coupling, a situation can occur where $\alpha_{1}/\alpha_{2}>\lambda_{1}/\lambda_{2}$ and therefore, a decrease in time step can cause divergence. This is for example the case for an air-water coupling [23].

Secondly, for an aspect ratio of $r=1$, we were able to extend the 1D results for both FVM-FEM and FEM-FEM to 2D in the following sense (see [5, 23]). In 2D, the iteration matrix $\Sigma$ is not easy to compute for several reasons. First of all, the matrices $\textbf{M}_{1}+\Delta t\textbf{A}_{1}$ and $\textbf{M}_{2}+\Delta t\textbf{A}_{2}$ are sparse block tridiagonal matrices, and consequently their inverses are not straight forward to compute. Moreover, the diagonal blocks of the same matrices are tridiagonal but their inverses are full matrices.

Due to these difficulties, we approximated the strictly diagonally dominant matrices $\textbf{M}_{1}+\Delta t\textbf{A}_{1}$ and $\textbf{M}_{2}+\Delta t\textbf{A}_{2}$ by their diagonal. Thus, we obtained an estimate of the spectral radius of the iteration matrix $\Sigma$. This estimator tends to the exact same limits as for the 1D case for both combination of discretizations.

We did not find a way to further extend these results to the high aspect ratio case. However, we will show now by numerical experiments that already the 1D formula (66) is a good estimator for convergence rates in 2D.

7 Numerical Results

We now present numerical experiments designed to illustrate the validity of the theoretical results of the previous sections. Firstly, we will confirm that the theoretical formula $|\Sigma|$ in (66) predicts the convergence rates in the 1D case. Secondly, we will show the validity of (66) as an estimator for the rates in the 2D case, we will also show that the theoretical asymptotics deduced in (67) and (LABEL:finallimitdeltax1D) match with the numerical experiments. Finally, we illustrate the validity of (66) as an estimator for the non linear thermal FSI test cases introduced in section 3.

7.1 1D FVM-FEM Results

We first compare the semidiscrete estimator $\beta$ in (18) with the discrete formula $|\Sigma|$ in 1D in (66) and experimental convergences rates. The latter are computed with respect to a reference solution $u_{ref}$ over the whole domain $\Omega$.

Figure 9 shows a comparison between $\beta$ and $|\Sigma|$ for $r=1$, $\Delta x=1/20$ and variable $\Delta t$. On the left we plot $\beta$, $|\Sigma|$ and the experimental convergence rates with $\Delta t/\Delta x^{2}\ll 1$ and on the right we plot the same but with $\Delta t/\Delta x^{2}\gg 1$. As can be seen, the experimental convergence rate matches exactly with the exact formula (66). Observe that $\beta$ is almost constant and represents the lower branch in figure 6. To arrive at the jump we would have to choose $\Delta t\gg 1$. We can conclude that the formulas for the convergence rates in 1D presented in the previous section are minimally better than the semidiscrete one proposed in [16] when $\Delta t/\Delta x^{2}\gg 1$. In the, less relevant case, $\Delta t/\Delta x^{2}\ll 1$ our formula also predicts the rates accurately, while the semidiscrete estimator deviates.

We now want to illustrate how $|\Sigma|$ in (66) gives the convergence rates and tends to the limits computed previously in (67) and (LABEL:finallimitdeltax1D). 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.

Figures 10 and 11 show the convergence rates for the interactions between air and steel and between water and steel, respectively. On the left we always have fixed $\Delta x_{1}$ and $r$, but variable $\Delta t$, whereas on the right we have fixed $\Delta t$ and $r$, but varying $\Delta x_{1}$. Each plot includes graphs for two different values of $r$. In figure 10 we choose $r=1$ and $r=100$ to illustrate the effect of a neutral or a high aspect ratio. In figure 11 we use $r=0.01$ and $r=1$ to illustrate how the rates are affected by a small or a neutral aspect ratio.

Again, $|\Sigma|$ gives the exact convergence rates. Moreover, one observes that the rates in 10a and 11a tend to 0 as predicted in (67) and in 10b and 11b to $\delta_{r}$ as predicted in (LABEL:finallimitdeltax1D). Furthermore, there is a roughly proportional relation between the convergence rate and the aspect ratio. For coupling with compressible flows, we typically have a high aspect ratio and therefore, the Dirichlet-Neumann iteration will be slowed down. Furthermore, this shows that it is very important to take the aspect ratio into account to make a reasonable prediction of the convergence rate at all.

7.2 2D FVM-FEM Results

We now want to demonstrate that the 1D formula (66) is a good estimator for the convergence rates in 2D. Thus, we now consider a 2D version of (16) consisting of two coupled linear heat equations on two identical unit squares, e.g, $\Omega_{1}=[-1,0]\times[0,1]$ and $\Omega_{2}=[0,1]\times[0,1]$. We use a non equidistant cartesian grid with aspect ratio $r$ on $\Omega_{1}$ and an equidistant grid on $\Omega_{2}$. In order to use (66) as an estimator we decided to take the equidistant mesh width on $\Omega_{2}$ as $\Delta x_{2}$ and the mesh width in $x$-direction on $\Omega_{1}$ as $\Delta x_{1}$.

As before, we present two real data examples described in table 1 and 2, namely the thermal interaction between air at $273K$ with steel at $900K$ and air at $273K$ with water at $283K$.

Figures 12 and 13 show the convergence rates for the interactions between air and steel and between air and water in 2D respectively. On the left we always plot the rates for fixed $\Delta x_{1}$ and $r$ with variable $\Delta t$. On the right we plot the behaviour of the rates for fixed $\Delta t$ and $r$ and varying $\Delta x_{1}$. As before, notice that each plot includes two cases for two different $r$ values. In figure 12 we choose $r=1$ and $r=100$ as in the 1D case (see figure 10) and in figure 13 we use $r=1$ and $r=1000$ to illustrate the effect of a neutral or a high aspect ratio. One observes that the convergence rates predicted by the one-dimensional formula (66) are almost exactly the ones observed in 2D. Thus, the 1D model problem case gives a very good estimator for the 2D model problem.

7.3 Thermal FSI Test Cases

Finally, we want to relate the results for the two nonlinear applications (the two cooling systems introduced in sections 3.1 and 3.2: the cooling of a flat plate and of a flanged shaft) to our analysis. Figure 14a shows the convergence behaviour for the flat plate and 14b for the flanged shaft. We plot the experimental convergence rates, the one-dimensional formula (66), the semidiscrete estimator (20) and the spatial limit $\delta_{r}$ specified in (LABEL:finallimitdeltax1D).

In order to apply the 1D formula (66) here, some assumptions need to be taken, since we partly have unstructured meshes and nonuniform temperatures. Thus, we assume air at $273K$ on the first subdomain with steel at $900K$ on the second subdomain for the cooling of a flat plate and air at $273K$ with steel at $1145K$ for the cooling of a flanged shaft. The density, heat capacity and heat conductivity of air and the density of steel are given in table 1. In addition, the heat conductivities and heat capacities of steel at $900K$ and $1145K$ are obtained from the nonlinear coefficient functions (3) and (4) by inserting $\Theta=900K$ or $\Theta=1145K$ respectively. This gives $\lambda=39.82$ and $c_{p}=1.3684e3$ for steel at $900K$ and $\lambda=39.8$ and $c_{p}=572.75$ for steel at $1145K$.

Furthermore, for the cooling of a flat plate, we take $\Delta x_{1}=9.3736e-5$ which is the width of the fluid cells touching the interface in the $y$-direction and $\Delta x_{2}=1.6667$ which is width of the structure cells in both directions. Thus, we have an aspect ratio of $r=1.7780e4$. On the other hand, choosing $\Delta x_{1}$ and $\Delta x_{2}$ for the cooling of a flanged shaft is more difficult due to the unstructured grids. In order to get an upper bound for the aspect ratio $r$, we choose $\Delta x_{1}=1.6538e-4$ which is the minimum width of all the fluid cells touching the interface in direction normal to the wall and $\Delta x_{2}=1.1364$ which is the maximum width of all the structure cells touching the interface tangential to the wall. This gives $r=6.8713e3$.

From figure 14a one observes with these choices, (66) predicts the rates accurately for the cooling of a flat plate. Note that the semidiscrete estimator $\beta$ does not show any change with $\Delta t$. Remember that $\beta$ is almost always constant, except for a short dynamic transition between $(\lambda_{1}/\lambda_{2})\sqrt{D_{2}/D_{1}}$ and $\lambda_{1}/\lambda_{2}$ as shown in figure 6. Here, we would have to choose a $\Delta t$ larger than $10E6$ to see the transition.

Finally, in figure 14b one can see that (66) predicts the convergence rates for the cooling of a flanged shaft to be only slightly smaller compared to the actual performance. This could be due to either the unstructured grids used or to the nonconstant temperature in the structure, which varies from room temperature to $1145K$. Again, $\beta$ is almost constant.

8 Summary and Conclusions

We considered the Dirichlet-Neumann iteration for thermal FSI and studied the convergence rates. To this end, we considered the coupling of two heat equations on two identical domains. We assumed structured grids on both subdomains, but allowed for high aspect ratio grids in one domain. An exact formula for the convergence rates was derived for the 1D case. Furthermore, we determined the limits of the convergence rates when approaching the continuous case either in space ($r\lambda_{1}/\lambda_{2}$) or time ([math]). This was confirmed by numerical results, where we also demonstrated that the 1D case gives excellent estimates for the 2D case. In addition, numerical experiments show that the linear analysis is relevant for nonlinear thermal FSI problems.

All in all, strong jumps in the coefficients of the coupled PDEs will imply fast convergence. In the domain decomposition context, the coupling will be slow because the material coefficients are continuous over all the subdomains, i.e, $\lambda_{1}=\lambda_{2}$, and therefore $\delta_{1}\sim 1$. For coupling of structures and compressible flows, the aspect ratio in the structure has to be taken into account, since the convergence rate is proportional to it. For the nonlinear cooling problems considered here, the convergence rate was still around 0.1 for large $\Delta t$. When encountering divergence anyhow, this can be solved by reducing the time step. Note that in a time adaptive setting, it is standard to allow for a feedback loop between the nonlinear solver and the time stepper.