A MATHEMATICAL MODEL AND NUMERICAL SOLUTION OF INTERFACE PROBLEMS FOR STEADY STATE HEAT CONDUCTION

We study interface (or transmission) problems arising in the steady state heat conduction for layered medium. These problems are related to the elliptic equation of the form Au :=−∇(k(x)∇u(x))= F(x), x ∈Ω⊂ R2, with discontinuous coefficient k = k(x). We analyse two types of jump (or contact) conditions across the interfaces Γδ = Ω1 ∩Ωδ and Γδ = Ωδ ∩Ω2 of the layered medium Ω := Ω1 ∪Ωδ ∪Ω2. An asymptotic analysis of the interface problem is derived for the case when the thickness (2δ > 0) of the layer (isolation) Ωδ tends to zero. For each case, the local truncation errors of the used conservative finite difference scheme are estimated on the nonuniform grid. A fast direct solver has been applied for the interface problems with piecewise constant but discontinuous coefficient k = k(x). The presented numerical results illustrate high accuracy and show applicability of the given approach.

The second-order homogeneous finite difference schemes for one-dimensional problems with discontinuous coefficients on uniform meshes using the jump conditions have been given first in [12].Later, various immersed methods have been developed for numerical solution of such interface problems (see [5,[9][10][11]14] and references therein).
In one of the previous studies (see [4]), an immersed finite element space is used to solve the elliptic interface problems by a finite volume element method.Special nodal basis functions are introduced in a triangle whose interior intersects with the interface so that the jump conditions are satisfied.Further, finite difference methods for elliptic equations of the form ∇(β(x)∇u) + κ(x)u(x) = f (x) in a region Ω in 1 or 2 space dimensions are developed in [9].Here Ω is assumed to be a simple region (e.g., a rectangle) and a uniform rectangular grid is used.Across the irregular surface Γ of codimension 1 contained in Ω the functions β, κ, and f are assumed to be discontinuous, and along Γ, the source f may have a delta function.An immersed interface elliptic problem with discontinuities and singularities in a circular domain has been studied in [10].On the contact surface, the flux and the potential (or temperature) here are taken discontinuous.
We consider more general problem where the domain consists of the following three parts: conductor-isolator-conductor (see, Figure 2.1(b)).Evidently, the interface problem is a limit case of this problem.We provide an asymptotic analysis of the problem conductor-isolator-conductor and show which interface problem is a limit case.Then we consider an interface problem with ideal contact conditions.Since thickness of the isolation is small enough with respect to the dimension of the domain isolation, our approach is based on conservative finite difference schemes not on uniform, but on a nonuniform, mesh with the same order of accuracy.
In this paper, we consider the following three types of interface problems: (P1) the interface problem for two-layered nonhomegeneous medium with ideal contact conditions (1.3), (P2) the interface problem for nonhomegeneous three-layered medium (conductorisolation-conductor) with ideal contact conditions (1.3), (P3) the interface problem for two-layered nonhomegeneous medium with continuous flux and discontinuous temperature across the interface.In the first part of the paper, we prove that the problem (P3) is a limit case of the problem (P2), when the thickness of the isolation tends to zero.In the second part, we derive the conservative finite difference schemes of a nonuniform mesh, which has a truncation error of order O(h 2 ).Then we obtain the finite difference schemes of orders O(h) and O(h 2 ) on the interface for Neumann transmission conditions.In the final part, we illustrate numerical results related to considered problems.

The mathematical model of layered conductivities without isolation: ideal contact
Let us consider the steady state heat conduction problem [u] x1=ξ = 0, in the domain Ω = Ω 1 ∪ Γ ξ ∪ Ω 2 , occupied by the homegeneous bodies Ω 1 and Ω 2 , with different heat conductivities along the Ox 1 -axis and k 2 (x) along the Ox 2 -axis.The piecewise constant coefficient k 1 (x) is assumed to be discontinuous (k (1)  1 = k (2) 1 ) and the coefficient k 2 (x) > 0 is assumed to be continuously differentiable function.The source function F(x) is assumed to be continuous.The contact conditions between the bodies Ω 1 and Ω 2 are assumed to be ideal ones via the interface (2.5) Hence we assume that the layered body occupies the domain Ω with the boundary ∂Ω := Then we can prove that the solution u ∈ C 2 (Ω 1 ∪ Ω 2 ) ∩ C(Ω) of the interface problem is also the solution of the variational problem (1.2) given by (2.4) for discontinuous coefficient k 1 (x).(1.2), for the discontinuous coefficient k 1 (x) given by (2.4).

Z. M. Seyidmamedov and E. Ozbilge 5
Proof.Let us multiply the both sides of (2.1) by v = v(x), integrate them on Ω 1 and Ω 2 seperately, then apply integration by parts, (2.6) Taking into account the homogeneous Dirichlet condition (2.3) and summing up the identities, we get for all v = v(x).Since u = u(x) satisfies the Dirichlet transmission condition [u] x1=ξ = 0, it is natural to require that arbitrary function v = v(x) also satisfies this condition.Using this condition in the first integral on the right-hand side, we obtain This integral idenditity with the Neumann transmission condition completes the proof.
This theorem shows the equivalence of problem (1.1), with the discontinuous coefficient k 1 (x) given by (2.4), and the transmission problem (2.1)-(2.3),although problem (1.1) does not contain any transmission condition.This suggests a possibility of construction of such a finite difference analogue of the interface problem, which has the similar structure.Specifically, the problem is to construct the homogeneous finite difference scheme, which has the same form for all mesh points, including ones on the interface Γ ξ .

The interface problem for layered conductivities with isolation: asymptotic analysis
Consider now the steady state heat conduction problem formed by the homogeneous conductive bodies occupying the domains with conductivities k (1)  1 (x) and consider k (2)  1 (x), respectively, and the isolation between Ω 1 and Ω 2 occupying the domain (Figure 2.1(b)), It is assumed that the conductivity of the isolation, with the thickness 2δ > 0, is small enough, γ = const 1.
We assume that heat transfers between conductivities and the isolation across the interfaces according to ideal contact conditions (3.3)-(3.4).
Introducing the piecewise continuous coefficient we can show, as in the proof of Theorem 2.1, that the solution u(x ) is also the solution of the variational problem (1.2), for given coefficient (3.9).
Z. M. Seyidmamedov and E. Ozbilge 7 In the case of the finite thickness isolation Ω δ , problem (3.1)-(3.5)represents an interface problem with two ideal contact conditions across the interfaces Γ − δ and Γ + δ .In practice, the isolator can be given as a thin boundary layer with small enough thickness δ > 0. If the value of the thickness is less than the mesh size h 1 along the direction Ox 1 , that is, 2δ < h 1 , then the interface conditions (3.3)-(3.4)cannot be approximated on this mesh.To derive a finite difference approximation of the interface problem (3.1)-(3.5)for this case, one needs to derive an asymptotic analysis of this problem, when δ → 0.
Proposition 3.1.The limit case, when δ → 0, of the transmission problem (3.1)-(3.5)with ideal contact conditions, is the following transmission problem with nonideal contact conditions: , where ξ ∈ (−δ,δ) and δ > 0 is an arbitrary small parameter and summing up the obtained integral identities imply We apply here the mean value theorem, dividing first the both sides of the above integral identity by l 2 = 0. Then we have where x 2 ∈ (0,l 2 ).We use here the second interface condition (3.3) for the flux on Γ − δ : F(x)dx = 0. (3.13) Going to the limit as −l 1 → −δ − 0 (the first integral drops), we obtain Since the parameter ξ ∈ (−δ,δ) is an arbitrary one, integrating the both sides of the above identity with respect to this parameter on (−δ,δ), we obtain Let us divide now the both sides by 2δ = 0. Then going to the limit as δ,γ → 0 and requiring σ := γ/(2δ) = const, we get where ξ → 0 as δ → 0. This condition can be rewritten in the following form: , respectively, by the same way, we can obtain the second limit condition across the interface Γ 0 = {(0, x 2 ) ∈ R 2 : x 2 ∈ (0,l 2 )} remains continuous, when the thickness tends to zero.However, as these conditions show, the temperature u(x) becomes discontinuous across the interface Γ 0 .The jump of the function u(x) across the interface is expressed via the flux ϕ 0 (x 2 ) on Γ 0 by formula (3.19).
The above asymptotic analysis shows that the limit case δ → 0 of the interface problem (3.1)-(3.5)with ideal contact conditions (3.3)-(3.4) is the interface problem (3.10) whose solution has discontinuity across the interface.In practice, this analysis is necessary, especially for the class of problems, in which the conductivity γ > 0 and the thickness 2δ > 0 are of the same order of small parameters.

Conservative finite difference schemes for interface problems on nonuniform mesh
Our goal here is a finite difference approximation of transmission problems (2.1)-( 2.3) and (3.10).For simplicity, we assume ξ = 0, so the interface Γ ξ of the domain Ω, shown in Figure 2.1(a), lies on Ox 2 -axis.The most existing numerical approaches use a uniform mesh and the interface points (ξ,x 2 ), x 2 ∈ (0,l 2 ) are assumed to be between mesh points x i0 − j < ξ < x i0 + j (see, e.g., [9]).We will develop here finite difference equations on a nonuniform mesh, assuming that the interface lies on mesh points.
Let us denote by a nonuniform rectangular mesh with mesh steps h (ip) j > 0, j = 1,2, assuming that x (0) Then, according to the above assumption, the interface mesh points are defined as follows: (0,x (i2) 2 ) : ( Below, we will use the following notations without index at the mesh point (x (i1) 1 ,x (i2) 2 ): For the first left and right finite differences, we will use the standard notations In a similar way, u x2 and u x2 can also be defined.The first finite difference corresponding to the mesh step p = 0.5(h − p + h + p ) will be denoted as follows: We define the mesh point (x (i1) 1 ,x (i2) 2 ) ∈ w h to be a regular mesh point, if the interface Γ ξ does not come between any points in the standard five-point stencil (Figure 4.1).In view of these definitions, the standard five-point conservative finite difference scheme approximating the elliptic operator (1.1) at the regular mesh point (x (i1) 1 ,x (i2) 2 ), with the coefficient k 1 (x), x ∈ Ω, is defined as follows [9,11]: Here the coefficients a p are defined, for example, by the formulas This scheme has the local truncation error of O(h 2 ),

Approximation of the Neumann transmission condition (2.2) on nonuniform mesh. We begin by the consideration of the finite difference approximation of the transmission problem (2.1)-(2.3) with ideal contact conditions.
To approximate the Neumann interface condition (2.2), we use first the simplest finite difference equation where Γ h ξ is the set of interface mesh points.Taking into account formulas (2.4) and (4.5), we get the following: Z. M. Seyidmamedov and E. Ozbilge 11 Let us now calculate the truncation error of the finite difference approximation (4.6).For this aim, we use Taylor's formula on the right-hand side of (4.8): (4.9) Taking into account the Neumann interface condition (2.2), we can eliminate the first and the second terms on the right-hand side.Hence and we get an O(h) truncation error along the interface Γ ξ .
Consider now the spectral case on an orthotropic material, when the heat conduction equation (2.1) has the form where k 2 (x) ∈ C(Ω) and k 1 (x) has the same form (2.4), that is, the heat conduction coefficient k (i) 1 (x) is continuous, but the coefficient is discontinuous.In this case, by using (4.11) on the right-hand side of (4.10), we get Taking into account the right-hand side of this expression, we consider the following approximation of the Neumann transmission condition: instead of (4.6).Here a 2 (x) = k 2 (x 1 ,x 2 − h 2 /2), according to (4.5).This scheme has the local truncation error O(h 2 1 ), as formula (4.12) shows.To analyse scheme (4.13), we divide both sides by h 1 = 0 and use the above notations for finite derivatives.Then the finite difference equation has the following canonical form: This is identical to the scheme (4.4).Thus, due to the conservativeness of the finite difference schemes, even in the case of discontinuous coefficient k 1 (x), the scheme (4.14) has the same form (4.4) at the interface mesh points.

Approximation of the transmission conditions (3.10).
Due to the discontinuity of the solution u(x) across the interface, we introduce the following: y − := y − − 0,x 2 , y + := y + + 0,x 2 y − = y + (4.15) on the mesh points along the interface Γ ξ .We approximate the interface conditions (3.10) by the following finite difference equations: where the coefficients are defined by (4.5).To estimate the truncation error, we use the above Taylor formula: By the first Neumann transmission condition (3.10), the local truncation error on the mesh points of the interface is of O(h 1 ).The same error has the second scheme (4.16).
In the case of the elliptic equation (4.11), we can use the above technique to construct the following schemes: (4.18) Both schemes (4.18) have a local truncation error of O(h 2 ) on the mesh points of the interface Γ ξ .

Numerical examples
The series of computational experiments are done to confirm the mathematical model of layered conductivities without isolation and with thin isolation with ideal contact conditions, as well as the expected accuracy of the presented finite difference approximations.We illustrate some results by the following two examples in which k 2 (x) = 1 is taken.Z. M. Seyidmamedov and E. Ozbilge 13 Example 5.1.In this example, we study the verification of the finite difference scheme (4.4) approximating the elliptic problem (2.1)-(2.3)with discontinuous coefficient k 1 (x).For this problem, we use the following exact solution: with appropriate Dirichlet condition, and the piecewise constant coefficient k 1 (x) is given by where Ω :={(x 1 ,x 2 ) ∈ R 2 : x 1 ∈ (0,1), x 2 ∈ (0,1)} and ξ = 0.5.For the function u(x 1 ,x 2 ), u(x 1 ,0. with appropriate Dirichlet condition, and the discontinuous coefficient (3.9) is given by where The domains Ω 1 , Ω δ , and Ω 2 are taken as and       To have an overview of p, the order of convergence (experimentally) of each scheme, we chose two different step sizes h 1 and h 2 .Then we calculated the corresponding errors, say E 1 and E 2 , respectively.We assumed that E 1 = ch (5.6) The above formula was used to estimate the order of convergence for each finite difference method.The results were obtained via various values of h.We observed that the order of convergence in the domain Ω δ and Ω 1 ∪ Ω 2 were 1 and 2, respectively for Example 5.2.

Figure 2 . 1 .
Figure 2.1.The geometry of layered conductivities without isolation and with thin isolation is given in part (a) and part (b), respectively.

Figure 5 . 1 .
Figure 5.1.For the first example, the obtained solution and the relative error are given in part (a) and part (b), respectively.
x2)∈Γ − δ = 0.008686, and [k 1 ∂u(x 1 ,x 2 )/∂x 1 ] (x1,x2)∈Γ + δ = 0.010345 are obtained.The results of the numerical solution on different meshes are given in Table 5.1.The fourth and fifth columns of the table show the values of the fluxes on the interfaces Γ − δ and Γ + δ , respectively, for δ = 0.025.These results are in agreement with Proposition 3.1.Relative error of the obtained numerical solution is given in Figure 5.2(b).The figures of the cross sections of the approximate solution and relative error with the plane x 2 = 0.5 are given in Figures 5.

3
(a) and5.3(b),respectively.The geometry of the behaviour of the relative error at the cross section x 2 = 0.5 for different thickness δ = h;2h;3h of the layer (isolation) is illustrated in Figure5.4.The minimum relative error, in all cases, corresponds to the central point x 1 = 0.5 of the isolation, is as expected.

Figure 5 . 2 .
Figure 5.2.For the second example, the exact solution and the relative error of the obtained numerical solution are given in part (a) and part (b), respectively.

(N+1)/ 2 ) 3 NFigure 5 . 3 .
Figure 5.3.For the second example, the cross section of the approximate solution and the relative error with the plane x 2 = 0.5 are given in part (a) and part (b), respectively.

Table 5 .
1. Results of the numerical solution on different meshes.