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

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 the interface. The heat equations are discretized using an implicit Euler scheme in time, whereas a finite element method on one domain and a finite volume method with variable aspect ratio on the other one are used in space. We provide an exact formula for the spectral radius of the iteration matrix. The formula indicates 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 very good estimator in 2D and even for nonlinear thermal FSI applications.


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 segregated 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 B Azahar Monge azahar.monge@na.lu.se 1 Centre for Mathematical Sciences, Lund University, Box 118, 22100 Lund, Sweden interaction of a flexible structure with a fluid has been analyzed in [35]. 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 some FSI couplings, which is why a lot of effort goes into convergence acceleration [10]. An alternative are optimized Schwarz methods [13,15,29] and the CHAMP scheme (Conjugate Heat transfer Advanced Multidomain Partitioned) which uses a generalized Robin (mixed) condition at the interface to accelerate the iteration [24], but overlapping domains. 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 between air and steel [4].
Our prime motivation here is thermal interaction between fluids and structures, also called conjugate heat transfer. In particular, we consider 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 [19,23], gas quenching, which is an industrial heat treatment of metal workpieces [17,33] or the cooling of rocket nozzles [20,21].
A convergence analysis of the Dirichlet-Neumann iteration for an FE-FE discretization of a steady heat transfer problem can be found in [32]. Asymptotically, it is found to be the quotient of heat conductivities. In this paper, we present a convergence analysis of the unsteady heat transfer problem. For this model, a 1D stability analysis was presented by Giles [16]. There, an explicit time integration method was chosen with respect to the interface unknowns. On the other hand, Henshaw and Chand provided in [18] 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 behavior mentioned above for 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 Δ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. [30,34]. 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 it is in general a nontrivial 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 [26,28], 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 [27]. In addition, the corresponding analysis when coupling finite volumes (FVM) with FEM is described in [5,26]. 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 the temporal limit turns out to be the ratio of the product of densities and heat capacities of the materials for FEM-FEM couplings instead. Moreover, we also include numerical results where it is shown that the one dimensional formula is a very good estimator for a 2D version of the coupled heat equations and for two nonlinear FSI models, namely the cooling of a flat plate and the cooling of a flanged shaft.
An outline of the paper now follows. In Sect. 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 Sect. 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 Sect. 4 and then analyzed in 1D in Sect. 5. In Sect. 6, an extension of the analysis to 2D and different discretizations are discussed. In Sect. 7, the analytical results are compared to linear and nonlinear numerical results.

Thermal FSI methodology
We consider the basic setting where on a domain Ω 1 ⊂ R d where d corresponds to the spatial dimension, the physics is described by a fluid model, whereas on a domain Ω 2 ⊂ R d , a different model describing the structure is used. The two domains 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 Γ ⊂ ∂Ω 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 problem, these conditions are that temperature and the normal component of the heat flux are continuous across the interface.

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 ρ, momentum m = ρv and energy per unit volume ρ E as: Here, enthalpy is given by being the pressure and γ = 1.4 the adiabatic index for an ideal gas. Furthermore, q f = (q 1 , q 2 ) T represents the heat flux and S = (S i j ) i, j=1,2 the viscous shear stress tensor with . The nondimensional dynamic viscosity μ is given by the Sutherland law with the Sutherland constant Su = 110 K 273 K . 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 = (γ − 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 Γ .

Structure model
Regarding the structure model, we will consider heat conduction only. Thus, we have the nonlinear heat equation for the structure temperature Θ where denotes the heat flux vector. For alloys, the specific heat capacity c p and heat conductivity λ are temperature-dependent and highly nonlinear.
As an example, we will consider the empirical model for the steel 51CrV4 suggested in [31]. This was obtained from measurements and a least squares fit. The coefficient functions are then For the mass density, one obtains ρ = 7836 kg/m 3 . Finally, on the boundary, we have Neumann conditions q s (x, t) · n(x) = q b (x, t).

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; where T is the fluid temperature and the structure temperature and

Discretization in space
Following the partitioned coupling approach, we discretize the two models separately in space. For the fluid model (1), we use a finite volume method, leading to the following equation for all unknowns on Ω 1 , collected in the vector u ∈ R N f : where h(u, Θ Γ ) represents the nonlinear finite element spatial discretization and its dependence on the temperatures on the discrete interface to the structure, here denoted by Θ Γ . In structural mechanics, the use of finite element methods is ubiquitious. Therefore, we will also follow that approach here. Using quadratic finite element, one obtains the following nonlinear equation for all unknowns on Ω 2 : Here, M is the mass matrix, also called heat capacity matrix for this problem and A is the heat conductivity and stiffness matrix, respectively. The vector ∈ R N s consists of all discrete temperature unknowns and q Γ b (u) is the discrete heat flux vector on the coupling interface to the fluid, whereas q f b corresponds to boundary heat fluxes independent of the fluid, for example at insulated boundaries.

Time discretization
In time, we use the implicit Euler method with constant time step Δt. For the system (9)-(10) we obtain

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 Ω 1 with Dirichlet data on Γ and (12) on Ω 2 with Neumann data on Γ . Thus, one gets for the kth iteration the two decoupled equation systems with some initial condition (t = 0)| Γ = 0 Γ . The iteration is terminated according to the standard criterion where τ is a user defined tolerance.

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 rate 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. We consider 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 [4,6,7]. 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) [22]. In the fluid, the DLR TAU-Code in its 2014.2 version is employed [14], 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.

Flow over a plate
The first test case is the cooling of a flat steel plate resembling a simple work piece [7]. It is initially at a much higher temperature than the fluid and then cooled by a constant laminar air stream, see Fig. 1.
The inlet is located 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. Regarding the initial condition in the structure, a constant temperature of 900 K at t = 0 is chosen throughout. To determine the Reynolds number, a reference length ofx re f = 0.2 m is chosen.
The grid is chosen Cartesian and equidistant in the structural part. In the fluid region the thinnest cells touch the boundary and then get 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 Γ and there are 9660 cells in the fluid region and n x × n y = 120 × 9 = 1080 elements with 121 × 10 = 1210 nodes in the region of the structure.
The left plot in Fig. 4 shows the convergence behaviour of the Dirichlet-Neumann iteration against the time step Δt. One observes that the convergence rates is roughly propor-

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) [36]. Here, we have a hot flanged shaft that is cooled by cold high pressured air coming out of small tubes, see Fig. 2. We assume symmetry along the vertical 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 a straight and uniform way at a Mach number of 1.2, as well as a freestream in x-direction of Mach 0.005. The Reynolds number, based on a reference length ofx re f = 0.02 m is Re = 2500 and the Prandtl number Pr = 0.72.
The grid, see Fig. 3, consists of 279,212 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 from [4]: first freestream values are set overall in the fluid and temperatures from a thermographic camera in the structure. Then 10 −5 s of real time are computed using a time step of Δt = 10 −6 s. The right plot in Fig. 4 shows the convergence behaviour of the Dirichlet-Neumann iteration against the time step Δt. The convergence rate is again about proportional to the time step size and the iteration is again convergent even for very large time steps. If we compare the rates for the two problems, we observe that for a given Δ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 cases with strong jumps in the material coefficients, as the ones presented here. To understand this better, we perform in the next section a convergence analysis for the case of two coupled linear heat equations.

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.

Model problem
The unsteady transmission problem reads as follows, where we consider a domain Ω ⊂ R d which is cut into two sub- where t ∈ [t 0 , t f ] and n m 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. D 1 and D 2 represent the thermal diffusivities of the materials and they are defined by where ρ m represents the density and c p m the specific heat capacity of the material placed in Ω m , m = 1, 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 semidiscrete case. Henshaw and Chand applied in [18] 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 (with dual variable k) 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 β: In the one dimensional case, the transverse Fourier mode k is zero. Then, for Δt small enough, we have tanh −1/ √ D 2 Δt ≈ −1 and tanh 1/ √ D 1 Δt ≈ 1 and therefore: On the other hand, for Δt big enough, we have (20) Figure 5 shows β as a function of Δt for k = 0. It is almost constant, except for a short dynamic transition between (λ 1 /λ 2 ) √ D 2 /D 1 and λ 1 /λ 2 . Finally, one observes in (20) that the convergence rates of the Dirichlet-Neumann iteration are given by the quotient of thermal conductivities for Δt large. This suggests that strong jumps in the thermal conductivities of the materials give fast convergence.

Space discretization
We now describe a rather general space discretization of the model problem. The core property we need is that the meshes of Ω 1 and Ω 2 share the same nodes on Γ as shown in Fig. 6. Furthermore, we assume that there is a specific set of unknowns associated with the interface nodes. Otherwise, we allow at this point for arbitrary meshes on both sides. Then, letting u (m) I correspond to the unknowns on Ω m , m = 1, 2, and u Γ to the unknowns at the interface Γ , we can write a general discretization of the first equation in (16) in a compact form as: To close the system, we need an approximation of the normal derivatives on Γ . For the FVM on Ω 1 , we approximate the normal derivative of u 1 with respect to the interface using second order one-sided finite differences: Now, let φ j be a nodal FE basis function on Ω 2 for a node on Γ we observe that the normal derivative of u 2 with respect to the interface can be written as a linear functional using Green's formula [34, pp. 3]. Thus, the approximation of the normal derivative is given by Consequently, the equation is a discrete version of the fourth equation in (16) and completes the system (21). Notice that the left hand side of (24) comes from (23) and the right hand side from (22). We can now write the coupled equations (21) and (24) as an ODE for the vector of unknowns u = u wherẽ

Time discretization
Applying the implicit Euler method with time step Δt to the system (24), we get for the vector of unknowns u n+1 = (u (1),n+1 I , u where Γ I + ΔtA Γ Γ .

Dirichlet-Neumann iteration
We now employ a Dirichlet-Neumann iteration to solve the discrete system (26). This corresponds to alternately solving the discretized equations of the transmission problem (16) on Ω 1 with Dirichlet data on Γ and the discretization of (16) on Ω 2 with Neumann data on Γ . Therefore, from (26) one obtains for the kth iteration the two equation systems to be solved in succession. Here, Γ I + ΔtA with some initial condition, here u n+1,0 Γ = u n Γ . The iteration is terminated according to the standard criterion u k+1 where τ is a user defined tolerance [3]. One way to analyze this method is to write it as a splitting method for (26) and try to estimate the spectral radius of that iteration by a norm. However, the results obtained in this way are much too inaccurate. For that reason, we now rewrite (27)-(28) as an iteration for u n+1 Γ to restrict the size of the space to the dimension of u Γ which is much smaller. To this end, we isolate the term u (1),n+1,k+1 I in (27) and u (2),n+1,k+1 I in the first equation in (28) and we insert the resulting expressions into the second equation in (28). Consequently, the iteration u n+1,k+1 where 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 convergence is described by the spectral radius of the iteration matrix Σ.

One-dimensional convergence analysis
So far, the derivation was performed for a rather general discretization. In this section, we study the iteration matrix Σ 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 Ω 1 = [−1, 0], Ω 2 = [0, 1]. For the FVM discretization, we consider a primal grid, i.e, we discretize Ω 1 into N 1 + 1 equal sized grid cells of size Δx 1 = 1/(N 1 + 1), and define x i = iΔx 1 , so that x i is the center of the cell i, see Fig. 7. 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 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 Ω 2 into N 2 + 1 equal sized cells of size Δx 2 = 1/(N 2 + 1). For the coupling between a compressible fluid and a structure, there would typically be a boundary layer in the fluid, meaning that the mesh would be very fine in direction normal to the boundary, implying where the only nonzero entry is located at the jth position, the discretization matrices are given by Fig. 7 Grid cells over Ω 1 and Ω 2 for the finite volume discretization and the finite element discretization respectively Note that the iteration matrix Σ is just a real number in this case and thus its spectral radius is its modulus. One computes S (1) and S (2) by inserting the corresponding matrices specified above in (32) and (33) obtaining where α 1 i j represents the entries of the matrix (α 1 I−ΔtA 1 ) −1 and α 2 i j the entries of (M 2 + ΔtA 2 ) −1 for i, j = 1, . . . , N 1 and i, j = 1, . . . , N 2 respectively. Observe that the matrices (α 1 I − ΔtA 1 ) and (M 2 + ΔtA 2 ) are tridiagonal Toeplitz matrices but their inverses are full matrices. The computation of the exact inverses could be performed based on the recursive formula presented in [12] which runs over the entries of the matrix and consequently, it is non trivial to compute α 1 N 1 N 1 , α 1 N 1 −1N 1 and α 2 11 this way. Due to these difficulties, we propose to rewrite the matrices (α 1 I − ΔtA 1 ) −1 and (M 2 + ΔtA 2 ) −1 in terms of their eigendecomposition: where the matrix V N has the eigenvectors of any symmetric tridiagonal Toeplitz matrix of dimension N as columns. The entries of V N 1 and V N 2 are not dependent on the entries of α 1 I − ΔtA 1 or M 2 + ΔtA 2 due to their symmetry. Moreover, the matrices Λ 1 and Λ 2 are diagonal matrices having the eigenvalues of α 1 I − ΔtA 1 or M 2 + ΔtA 2 as entries respectively. These are known and given e.g. in [25, pp. 514-516]: The entries α 1 N 1 N 1 , α 1 N 1 −1N 1 and α 2 11 of the matrices (α 1 I − ΔtA 1 ) −1 and (M 2 + ΔtA 2 ) −1 , respectively, are now computed through their eigendecomposition resulting in Now, inserting (39), (40) and (41) into (34) and (35) we get for S (1) and S (2) : With this we obtain an explicit formula for the spectral radius of the iteration matrix Σ as a function of Δx 1 , Δx 2 and Δt: To simplify this, the finite sums N 1 i=1 sin 2 (iπΔx 1 ) and can be computed. We first rewrite the sum of squared sine terms into a sum of cosine 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: Inserting (48) and (49) into (47) we get after some manipulations We could not find a way of simplifying the finite sum (44) because Δx 2 depends on N 2 (i.e., Δx 2 = 1/(N 2 + 1)). However, (50) is a computable expression that gives the exact convergence rates of the Dirichlet-Neumann iteration for given Δt, Δx m , α m and λ m , m = 1, 2.
We are now interested in the asymptotics of (50) with respect to both spatial and temporal resolutions. This corresponds to the computation of two different limits: Δt → 0 for a fixed Δx 1 and Δx 1 → 0 for a fixed Δt. As an alternative, one could reformulate (50) in terms of c := Δt/Δx 2 1 and compute the limits c → 0 and c → ∞. Both choices give the same results because for a fixed Δx 1 , if Δt → 0, then c → 0 and for a fixed Δt, if Δx 1 → 0, then c → ∞.
For simplicity, we compute the asymptotics of (50) for Δt → 0 and Δx 1 → 0 with Δx 2 = r · Δx 1 where r := Δx 2 /Δx 1 is a fixed aspect ratio. This is motivated by the assumption that we have matching nodes at the interface. Thus, the resolution in the fluid in direction tangential to the wall is the same as the resolution in the structure. This means that the aspect ratio of the left subdomain cells in 2D corresponds to the ratio of grid spacings between the two subdomains in 1D. This is illustrated in Fig. 8. We obtain: (51) Fig. 8 Relation between the aspect ratio of the left subdomain cells in 2D and the ratio of grid spacings between both subdomains in 1D lim To simplify (52), it is well known that the finite sums cos 2 (iπΔx 1 ) can be computed by using Euler's formula to convert them into geometric sums. We thus obtain after some calculations: In order to compute the third sum, we rewrite the sum of squared cosine terms into a sum of sine terms using the identity cos 2 (x/2) = (1 + cos(x))/2 and then apply the same technique: Inserting (53) and (54) into (52) we get lim From the result obtained in (51) we can conclude that the convergence rate goes to zero when the time step decreases and therefore, the iteration will be fast for Δt small and can always be made to converge by decreasing Δt. This is consistent with the behavior of the cooling of a flat plat and the flanged shaft presented earlier in Fig. 4.
On the other hand, from the spatial asymptotics (55) we can observe that strong jumps in the thermal conductivities of the materials placed in Ω 1 and Ω 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.

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 for the fluid 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 Σ are given by [26,28]: When we compare these with the asymptotics obtained with FVM-FEM discretizations (51)-(55), we observe that while the spatial limit is the same, the temporal limit does not match. This arises from differences in the matrix S (1) in (30). In the FEM-FEM context, the matrices S (1) and S (2) lead to the same expression with only different material coefficients (α 1 , α 2 , λ 1 , λ 2 ). Because of this, the limits of ρ(Σ) are quotients of those coefficients. However, the situation is different in the FVM-FEM context. There, the matrix S (1) in (32) is missing several mass matrices if we compare it with S (2) in (33). This unsymmetry between S (1) and S (2) causes that the This implies that, opposed to the FVM-FEM case, where convergence can always be achieved by decreasing the time step, for an FEM-FEM coupling, a situation can occur where α 1 /α 2 > λ 1 /λ 2 and therefore, a decrease in time step can even cause divergence. This is for example the case for an air-water coupling [26].
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,26]). In 2D, the iteration matrix Σ is not easy to compute for several reasons. First of all, the matrices M 1 + ΔtA 1 and M 2 + ΔtA 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 M 1 + ΔtA 1 and M 2 + ΔtA 2 by their diagonal. Thus, we obtained an estimate of the spectral radius of the iteration matrix Σ. 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 (50) is a good estimator for convergence rates in 2D.

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 for |Σ| in (50) predicts the convergence rates in the 1D case. Secondly, we will show the validity of (50) as an estimator for the rates in the 2D case, we will also show that the theoretical asymp-totics deduced in (51) and (55) match with the numerical experiments. Finally, we illustrate the validity of (50) as an estimator for the nonlinear thermal FSI test cases introduced in Sect. 3.

Results in 1D
We first compare the semidiscrete estimator β in (18) with the discrete formula |Σ| in 1D in (50) and experimental convergences rates. The latter are obtained from implementing the Dirichlet-Neumann method (27)- (28). The results are then compared to a reference solution u re f over the whole domain Ω, obtained by choosing a tolerance of 1e−10 as a termination criterion. Figure 9 shows a comparison between β and |Σ| for r = 1, Δx = 1/20 and Δx = 1/500 and varying Δt. On the left we plot β, |Σ| and the experimental convergence rates for the FVM-FEM approach described in Sect. 5 and on the right for the FEM-FEM approach mentioned in Sect. 6. As can be seen, the experimental convergence rate matches exactly with the exact formula (50). Observe that β is almost constant and presents the same behavior as in Fig. 5. We can conclude that the formulas for the convergence rates in 1D presented in Sect. 5 match the semidiscrete one proposed in [18] when Δt/Δx 2 1. In the, less relevant case, Δt/Δx 2 1 our formula also predicts the rates accurately, while the semidiscrete estimator deviates according to (19). Finally, Fig. 9 also illustrates the differences in the temporal limit when employing different combinations of spatial discretizations as explained in previous section. In the FVM-FEM case the limit is 0 [see (51)] and in the FEM-FEM case it is α 1 /α 2 [see (56)].
The difference to the semidiscrete analysis in [18] stems from different limits taking place. The semidiscrete analysis implicitly assumes that first a limit Δx to zero has taken place for Δt fixed. Thus, a limit first Δt, then Δx to zero is not addressed by it. This can be seen in Fig. 9 in the follow-  Table 2 Temporal and spatial asymptotics of (50) for the thermal interaction of air at 273 K with steel at 900 K, water at 283 K with steel and air with water Air-water 0 0.0419 · r ing way: for fixed Δx, letting Δt become very small causes the convergence rate to move into the direction predicted by the fully discrete analysis. However, then keeping this very small Δt fixed and decreasing Δx moves that rate back in the vicinity of β.
We now want to illustrate how |Σ| in (50) gives the convergence rates and tends to the limits computed previously in (51) and (55). To this end, we present two real data examples. We consider here the thermal interaction between air at 273 K with steel at 900 K and water at 283 K with steel at 900 K. Physical properties of the materials and resulting asymptotics for these two cases are shown in Tables 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 have always fixed Δx 1 and r and vary Δt, whereas on the right we have fixed Δt and r , and vary Δx 1 . Each plot includes graphs for two different values of r . In Fig. 10 we choose r = 1 and r = 100 to illustrate the effect of a neutral or a high aspect ratio. In Fig. 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, |Σ| gives the exact convergence rates. Moreover, one observes that the rates on the left plots in Figs. 10 and 11 tend to 0 as predicted in (51) and on the right plots in Figs. 10 and 11 to δ r as predicted in (55).
Before ending this subsection, we want to illustrate the relation between the convergence rates and the aspect ratio r . To this end, Fig. 12 shows the convergence rates for the interaction between air and steel. In the left plot we have chosen Δx 1 to be coarse and on the right one to be fine. This explains why the convergence rates on the right plot are closer to the spatial limit δ r . 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.

2D FVM-FEM results
We now want to demonstrate that the 1D formula (50) 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, Ω 1 = [−1, 0] × [0, 1] and Ω 2 = [0, 1] × [0, 1]. We use a non equidistant cartesian grid with aspect ratio r on Ω 1 and an equidistant grid on Ω 2 . In order to use (50) as an estimator we decided to take the equidistant mesh width on Ω 2 as Δx 2 and the mesh width in x-direction on Ω 1 as Δx 1 .
As before, we present two real data examples described in Tables 1 and 2, namely the thermal interaction between air at 273 K with steel at 900 K and air at 273 K with water at 283 K. Figures 13 and 14 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 Δx 1 and r with variable Δt, whereas on the right we have fixed Δt and r and varying Δx 1 . As before, each plot includes two different values of r . In Fig. 13 we choose r = 1 and r = 100 as in the 1D case (see Fig. 10) and in Fig. 14 we use r = 1 and r = 1000 to illustrate the effect of a neutral or a high aspect ratio. To compute β we use the transverse Fourier mode k = iΔy, i = 0, 1, 2, . . . , N y that maximizes it [see (18)]. One observes that the convergence rates predicted by the one-dimensional formula (50) are almost exactly the ones observed in 2D. Thus, the 1D case gives a very good estimator for the 2D model problem.

Thermal FSI test cases
Finally, we want to relate the results for the two nonlinear applications (the two cooling systems introduced in Sects. 3.1 and 3.2: the cooling of a flat plate and of a flanged shaft) to our analysis. The left plot in Fig. 15 shows the convergence behaviour for the flat plate and the right one for the flanged shaft. We plot the experimental convergence rates, the onedimensional formula (50), the semidiscrete estimator (18) for the maximizing Fourier mode and the spatial limit δ r specified in (55).
In order to apply the 1D formula (50) here, some assumptions need to be made, since we partly have unstructured meshes and nonuniform temperatures. Thus, we assume air at 273 K on the first subdomain with steel at 900 K on the second subdomain for the cooling of a flat plate and air at 273 K with steel at 1145 K 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 900 and 1145 K are obtained from the nonlinear coefficient functions (3) and (4) by inserting Θ = 900 K or Θ = 1145 K respectively. This gives λ = 39.82 and c p = 1.3684e3 for steel at 900 K and λ = 39.8 and c p = 572.75 for steel at 1145 K.
Furthermore, for the cooling of a flat plate, we take Δx 1 = 9.3736e−5 which is the width of the fluid cells touching the interface in the y-direction and Δx 2 = 1.6667 which is the width of the structure cells in both directions. Thus, we have an aspect ratio of r = 1.7780e4. On the other hand, choosing Δx 1 and Δ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 Δ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 Δ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 the left plot in Fig. 15 one observes with these choices that (50) predicts the rates accurately for the cooling of a flat plate. Note that the semidiscrete estimator β does not show any change with Δt. Remember that β is almost always constant, except for a short dynamic transition between (λ 1 /λ 2 ) √ D 2 /D 1 and λ 1 /λ 2 as shown in Fig. 5. Here, we would have to choose a Δt larger than 1e6 to see the transition.
Finally, on the right plot in Fig. 15 one can see that (50) 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 1145 K. Again, β is almost constant.

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 λ 1 /λ 2 ) or time (0). This was confirmed by numerical results, where we also demonstrated that the 1D case gives excellent estimates for the 2D case. In addition, numer- All in all, strong jumps in the coefficients of the coupled PDEs imply fast convergence. On the other hand, the coupling iteration will be slow when the material coefficients are continuous over all the subdomains, i.e, λ 1 = λ 2 , and therefore δ 1 ∼ 1. For coupling of structures and compressible flows, the aspect ratio in the fluid 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 Δ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.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.