Solving stationary inverse heat conduction in a thin plate

We consider a steady state heat conduction problem in a thin plate. In the application, it is used to connect two cylindrical containers and fix their relative positions. At the same time it serves to measure the temperature on the inner cylinder. We derive a two dimensional mathematical model, and use it to approximate the heat conduction in the thin plate. Since the plate has sharp edges on the sides the resulting problem is described by a degenerate elliptic equation. To find the temperature in the interior part from the exterior measurements, we formulate the problem as a Cauchy problem for stationary heat equation. We also reformulate the Cauchy problem as an operator equation, with a compact operator, and apply the Landweber iteration method to solve the equation. The case of the degenerate elliptic equation has not been previously studied in this context. For numerical computation, we consider the case where noisy data is present and analyse the convergence.


Introduction
Determination of the temperature on an inaccessible part of the boundary, by given available measurements on another part of the boundary, has many industrial applications related to engineering and science.Some of these applications are metal quenching [23], combustion [19], steel industry [34,38], brakes [40], among many others.They are referred to as inverse heat conduction problems.This area of research have been studied extensively, in 1D setting, by various authors, such as [3,4,11,14,15,21,39].The Cauchy problem for the steady-state anisotropic heat conduction in 2D and 3D has also been considered in [37].
The manuscript belongs to Applications of PDEs edited by Hyeonbae Kang.
In our research paper, we consider two cylindrical containers connected by a thin plate.This situation models a thermos bottle where a liquid is contained in the interior cylinder and the outer cylinder is a protective shell.The thin plate allows for heat conduction between the two cylinders, which allows us to monitor the temperature of the interior cylinder.The thin plate that connects the two cylinders is a 3D volume which is denoted by A, as shown in Fig. 1 We assume that c(x) is a differentiable and a continuous function, that behaves like a linear function near the end points x = 0 or x = b, and that c(0) = c(b) = 0.The space of function that contain the solution to the problem is denoted by V ( ).More details are given in Sect. 2. The domain is illustrated in Fig. 2.
The main goal of this research is to find the temperature in the interior part from the exterior measurements.We approximate the heat conduction in the thin plate by the following two dimensional steady-state heat conduction problem, with no heat generation, as follows where T is subject to the boundary conditions Here k is the thermal conductivity [W/m• • C], h is the convection heat transfer coefficient [W/m 2 • • C], c(x) is the thickness of the plate.We assume that the ambient temperature T ∞ = 0, because the equation is linear.The Eq. (1.1) is a degenerate elliptic equation since the coefficient c(x) vanishes at x = 0 or x = b.We consider the following Cauchy problem for stationary heat equation where f and g are specified temperature and net heat flux at 0 .Our aim is to compute both the temperature ζ and net heat flux η at 1 .This problem is ill-posed [20], i.e. a small perturbation in the Cauchy data results in a large error in the solution.Thus classical numerical methods cannot be used to solve such a problem.Therefore, in previous works regularization methods have been introduced to stabilize the solution.Such methods as the Tikhonov and the Fourier regularization methods [6,13], Methods for filtering data such as Wavelet and Fourier method [16], Boundary element method [29], the Landweber iterative method [31], and the Conjugate gradient method [35] have been used to analysed such problems.The alternating algorithm was initially proposed in [26,27], for solving Cauchy problems for linear elliptic equations.The method consists of solving a sequence of boundary value problems.We alternate between Dirichlet and Neumann conditions on a part of the boundary.Previously, it has been developed in various directions, see [2,9,12,24,30,32].In  1 The description of the engineering situation we are studying with two cylinders connected by a thin plate A. The plate is shown on the right Fig. 2 Description of the two dimension domain with a thickness given by 2c(x), whose boundary is divided into four parts 0 , 1 , 2 and 3 .Here, the solution is unknown on 1 , and 0 is where the solution is known.The boundary data (1.2) is given on 2 and 3 y particular [32], it was shown that the Dirichlet-Neumann algorithm is equivalent to solving Landweber iterations for an operator equation.The alternating algorithm has the advantage that they involve only boundary value problem of the same differential equation with different boundary conditions.The corresponding operator equation is also defined by solving the same differential equation with different boundary conditions.This means that we can easily solve problems with variable coefficients in domains with a complicated geometry.The regularizing properties of the algorithm are achieved by an appropriate choice of the boundary conditions [26].
In this work, we consider several methods for solving the Cauchy problem for the stationary heat equation.All of them involve the following two well-posed problems: where and H 1/2 are trace spaces for the degenerate Sobolev spaces connected with the quadratic form associated with the operator in given in (1.3).The definition of these spaces can be found in Sect.2.2.
The first method is the Dirichlet-Neumann Alternating algorithm which can be described as follows: (1) The first approximate temperature T 0 is obtained by solving (1.4), where η is an arbitrary initial approximation of the net heat flux on 1 .(2) Having constructed T 2 j , we find T 2 j+1 by solving (1.5) with ζ = T 2 j on 1 .
(3) We then determine T 2 j+2 by solving (1.4) ∂ y on 1 .The problems (1.4) and (1.5) are well-posed as it is shown in Sect.2.2.It is proven in Sect.2.4, that the alternating algorithm is convergent.Since the convergence of this method is slow we consider alternative methods.The goal is to develop frame work of operator equations that can enable us to define the scalar products which helps in computation of the adjoint operator and also appropriate norms.Thus, faster convergence method can be considered.
The second method for solving the Cauchy problem (1.3) is based on rewriting the problem as an operator equation, see [4,5,8,10,32] for similar work.We introduce the operator where T is the solution of problem (1.4) with f = 0. Now the Cauchy problem (1.3) can be written as an operator equation, with compact operator where ξ = g − ∂ y T (x, 0) and T solves problem (1.4) with η = 0.It is important to have a good description of the adjoint operator K * for K .The operator K * can be describe as follow where T is the solution to the problem (1.5) with ζ = 0.The operators K and K * are explained in more details in Sect.3.
To solve problem (1.6), the operator K can be approximated by a matrix, and traditional methods such as Conjugate gradient method can be applied.Since K can be approximated by a matrix explicitly we can also use Tikhonov regularization method and truncated singular value decomposition (TSVD).In our previous work [5], the Landweber iteration method has been used to solve different operator equation.This method was initially proposed by Landweber in [28].We wish to use the same method to solve the operator equation (1.6).
In this paper we analyse convergence of the Landweber iteration method, defined as: Let η 0 = 0 be the starting guess and compute for m = 0, 1, 2, . . ., where α is the relaxation parameter.It is well known that Landweber iterations (1.7) are convergent provided that (1.8) The proof is quite similar to that for the Richardson iteration in [36,Example 4.1].
The following is an outline of the paper.In Sect.1.1, we derive the mathematical model for the problem considered.In Sect.2, we analyse the properties of the direct problem and show that it's well-posed.We give a proof for the convergence of the Dirichlet-Neumann alternating algorithm in Sect.2.4.In Sect.3, we reformulate the Cauchy problem (1.3) as an operator equation, and show that it is well-defined.We also introduce the Landweber iterations method for solving the operator equation in the same section.We discretize the differential equation (1.11) using finite difference scheme, and also, present numerical results in Sect. 4. Finally, in Sect.5, we analyse the results, and draw some conclusion.

Derivation of the mathematical model
In this section, we derive the mathematical model of heat transfer in the thin plate, and also it's related boundary conditions.
Let us consider the two dimensional domain shown in Fig. 2. Pick a small volume element of length y and width x whose thickness is 2c(x).To derive the mathematical model, we consider an energy balance equation on the volume element where heat flow in both x− and y− direction as follows; where Q is the rate of temperature change with respect to time, q x and q y is obtained using Fourier's law of heat conduction and finally, q c is the convection heat transfer given in terms of Newton's law of cooling, see [22], we obtain where C p is the specific heat capacity [J/kg• • C], ρ is the density [kg/km 3 ], A x is the cross section area of the thin plate as a function of x, A y is the cross section area of the thin plate as a function of y, and A c is the cross section area for convective heat transfer.Dividing both sides of (1.9) by y • x and taking the limits y → 0 and x → 0, we obtain where k is the thermal conductivity coefficient and h is the convection heat transfer coefficient.
For the steady state heat transfer, and assuming that T ∞ = 0, we obtain Similarly, we derive the boundary conditions on 2 and 3 , as shown in Fig. 2, and we obtain lim We note that their exist a limit at x = 0 and x = b, since c(x) tend to zero at the end points.This does not necessarily mean that the limit of ∂ T ∂ x , as x approaches the boundary, is zero, or even bounded.

Properties of the auxiliary problems
In this section, we compute solutions of (1.11) subject to the boundary condition (1.12) using the method of separation of variables.We define the bilinear form associated to the two auxiliary problems (1.4) and (1.5), and show that the weak formulation of these problems are well-posed.We also introduce two norms for each of the spaces H −1/2 ( 0 ) and H −1/2 ( 1 ), and show that they are equivalence.Finally, we state the theorem for the convergence of the Dirichlet-Neumann alternating algorithm.

Separation of Variables
To find solution of the differential equation (1.11) together with the boundary conditions (1.12), we will use the method of separation variables.Let (2.13) Substituting (2.13) in the differential equation (1.11) and separating variables, we get and where λ is a constant.Taking into account the boundary conditions (1.12) we arrive at the following spectral problem (2.16) Due to the linear behaviour of the coefficient c(x) near the end points the spectrum of problem (2.16) consist of simple positive eigenvalues which are denoted by λ j , j = 1, 2, . .., we enumerate them in the increasing order, i.e.
The corresponding eigenfunction we denote by ψ j and they are chosen to satisfy the biorthogonality conditions (2.17) Multiplying both side of (2.16) by ψ j (x), integrating the first term by parts and then applying bi-orthogonality conditions and the boundary conditions, we obtain ) Next, we solve Eq. (2.15) φ yy − λ j φ = 0 with boundary conditions φ j (0) = M j and φ j (a) = N j .As a result we get It is convenient to have here the factor λ 1/4 j for normalization.The general solution becomes (2.21) The constants M j and N j will be found later by using boundary conditions for y = 0 and y = a.

Solvability of the two auxiliary problems
In this section, we introduce the bilinear form associated to the operator and show that the two auxiliary problems (1.4) and (1.5) are well-posed.
In what follows, we use the Green's type identity where Our aim is to express this form in terms of the coefficients N j and M j .Using the representation (2.21) and bi-orthogonality conditions (2.17)-(2.19),we obtain the following formula In order to evaluate the right-hand side here we start from where Then we integrate each term as follows: The first term is integrated by Similarly, we evaluate the second term Finally, we integrate the last term Combining these formulas, we obtain where This implies the existence of α 1 , α 2 > 0 for each β > 0 such that Using the inequality we obtain where the constant c 1 and c 2 are positive constants, and are independent of M j and N j .Let us evaluate the second term in (2.23).We have We integrate each term as follows: We treat the first term by Similarly, the second term is evaluated as Finally, the last term has the form Combining the above formulas, we obtain where Since the function G(z, X , Y ) is positive, for z > 0 and X 2 + Y 2 > 0, we have Now using (2.27), and the positivity of B(z) for z > 0, we get for some positive constants C 1 and C 2 .Combining (2.28) and (2.32), we obtain Let us introduce the spaces of functions in .The space V ( ) consists of functions defined on having a representation given by (2.21), where φ j are defined by (2.20).We supply this space with the norm Let also H 1/2 (0, b) be the space of functions with the norms , and H −1/2 (0, b) be the space of functions with the norms and and hence both T (x, 0) and T (x, a) belong to H 1/2 ( 0 ) and H 1/2 ( 1 ).In order to emphasis the fact that we are dealing with the trace of function from V ( ) on 0 ( 1 ), we shall use the notation H 1/2 ( 0 ) and Furthermore, and Similarly, and To give a clear notation for the dual spaces on 0 ( 1 ), we shall use H −1/2 ( 0 ) and H −1/2 ( 1 ) instead of H −1/2 (0, b).More detail information about function spaces can be found in [33].
The following lemma shows that the problem (1.4) is well-posed.
Lemma 2.1 Let f ∈ H 1/2 ( 0 ) and η ∈ H −1/2 ( 1 ).Then the problem (1.4) has a unique solution T ∈ V ( ) satisfying the estimate and From (2.21) and the second equation in (1.4), we get (2.41) Similarly using the third relation in (1.4), we obtain The auxiliary problem (1.5) can be analysed similarly.The corresponding solvability results is the following.
. Then the problem (1.5) has a unique solution T ∈ V ( ) satisfying the estimate (2.45)

Equivalence of norms
In [1,5], we use different norms for the trace functions on 0 and 1 .These norms are more suitable for the case when we can not write the solutions explicitly, by using separation of variables.These norms are defined by using only the bilinear form γ (T 1 , T 2 ) given in (2.22).Their definitions can be easily extended to more general degenerate elliptic operators and domains.
Let us introduce this norms and their corresponding inner products.
, and T 1 , T 2 ∈ V ( ) be the solution of problem (1.4) where η 1 = η and η 2 = η on 1 , and f = 0 on 0 .An inner product on the space H −1/2 ( 1 ) is define as The corresponding norm is given by η γ, 1 = η, η 1/2 1 .Similarly, we can introduce the inner product for g on 0 in terms of the solution of the well-posed problem (1.5).

.49)
Proof Applying Green's type identity, we have (2.50) Let T 1 ∈ V ( ) be a solution to problem (1.4) with f = 0 on 0 .Using the representation (2.21) the solution of T 1 (x, y) is given by where M j = 0 at y = 0 on 0 .Computing the derivative of (2.51) with respect to y at y = a, we obtain sinh λ j a ψ j (x). (2.52) The solution obtain in (2.52) is equals to (2.53) Equating (2.52) into (2.53), and identifying the coefficient gives, Substituting (2.54) in (2.51), we get (2.55) The solution of T 1 (x, y) at y = a is given by (2.56) Substituting (2.56) in (2.50), and using the fact that we have (2.57) By applying bi-orthogonality conditions (2.17) in (2.57), we obtain

The convergence of the Dirichlet-Neumann alternating algorithm
We now state the theorem of convergence of the Dirichlet-Neumann algorithm for solving the Cauchy problem (1.3) with the exact data, described in Sect. 1.Let {T j } ∞ j=0 be the iteration describe in the Dirichlet-Neumann algorithm shown in the introduction.This iterations depends on f , g and η.
Proof The proof of this theorem is similar to that in [1,9] and therefore we do not present it here.

The operator equation
In this section, we define the operator equations that are used in the numerical computations.Definition 3.1 Let η ∈ H −1/2 ( 1 ).We define the operator where T ( f ,η) is the solution of problem (1.4).We introduce the notation T ( f ,η) so that we can refer to solutions of the problem (1.4) with different boundary conditions f , η, e.g.T (0,η) solves (1.4) with f = 0.
Let us express the operator K in terms of the basis function ψ j (x).From (2.55), we have Computing the derivative of (3.59) with respect to y at y = 0, we obtain The Cauchy problem (1.3) can be written in terms of an operator equation with respect to η as follows K η = g − ∂ y T ( f ,0) (x, 0).(3.61) Lemma 3.2 Let K be given as in Definition 3.1.The expression (3.61) can be written as follows Proof Using the general solution (2.21), we obtain sinh λ j a + sinh λ j y sinh λ j a ψ j (x), (3.63)where M j = N j cosh λ j a at y = a.From (3.63), we compute at y = 0, and we obtain As a result we get, j N j cosh λ j a. (3.64) From (3.64) we obtain the coefficient N j as follows Evaluating the derivative of (3.66) with respect to y at y = 0, we obtain For more detail information on how to write a Cauchy problem as an operator equation, one can see for instance [5,17].
Let us now state and proof the lemma for the adjoint operator K * .Lemma 3.3 K be defined as in Definition 3.1, the adjoint operator where we use the inner products •, • 1 and •, • 0 is given by where T 0 is the solution to problem (1.5) with ζ = 0 on 1 .
Proof Using the general solution (2.21), we obtain where N j = 0 at y = a.Differentiating (3.69) with respect to y at y = 0, we obtain solution for g(x) as follows Differentiating (3.74) with respect to y at y = a we obtain By applying Lemma 2.5 and (3.60), we evaluate the inner product (1.5).Similarly, applying Lemma 2.5 and (3.75) we compute the inner product where T 0 ∈ V ( ) solves problem (1.4).Clearly from (3.76) and (3.77) we can see that Let us now state and proof a lemma for compactness of the operator K and also, compute its norm.

Lemma 3.4 The operator K defined in Definition 3.1 is compact and its norm in the spaces
Proof From (3.60) and (3.75), we have Therefore, This implies that the operator K is compact.Also, applying Lemma 2.5 and (3.60) we calculate the norm This assertions proves the lemma.

Landweber iterative method
We now define the iterative method for solving the problem where ξ = − ∂ y T ( f ,0) (x, 0) defined in (3.61).We introduce Landweber iterative method for solving the ill-posed problem (3.79).Let η 0 = 0 be the initial approximation, compute for m = 0, 1, 2, . . ., where α is the relaxation parameter.It is well known that the iteration (3.80) is convergent provided that (3.81) The proof is similar to that for the Richardson iteration in [36,Example 4.1].By Lemma 3.4 this inequality can be written as Similarly, the iterations (3.80) converges if we have exact right hand-side in (3.79), see [18,Theorem 6.1].As was shown in [5], the Landweber iterative method with α = 2/||K || 2 is equivalent to the Alternating iterative method.This explains the slow convergence.
As the stopping rule for the Landweber iterations we use the Discrepancy principle [18].The Discrepancy principle is applied as follows: Let ξ ∈ H −1/2 ( 0 ) be the exact data and Consider the Landweber iterations where 0 < α < 2 cosh 2 √ λ 1 a.Let τ > 1.The iterations are terminated when The stopping index is denoted by Note that m(δ, ξ δ ) → ∞ as δ → 0. Thus the following theorem [25, Theorem 2.15 and Theorem 2.19] guarantees convergence.
Then the Landweber iterations, together with the Discrepancy principle, satisfies Note that for the discrepancy principle in (3.84) the residual norm, i.e. the left side can be computed in each iteration.For the numerical experiments we know the exact data ξ .Thus, we can compute δ by evaluating the norm in (3.82).For practical application the exact data is unknown.Since we assume normally distributed random noise it is more natural to assume a bound of (3.82) in the L 2 norm.However, the point of the experiment section is to illustrate the theoretical part of the paper.Thus, we use the norms in H −1/2 .Note however that the L 2 norm is stronger than H −1/2 norms.Thus we could implement the discrepancy principle (3.84) using L 2 norm instead.

Numerical results
In this section, present numerical experiments for the iterative method analysed in the previous section.To conduct the experiments, we implement a finite difference method for solving the two auxiliary problems (1.4) and (1.5).Let a and b be positive numbers and recall that the domain is given as To carry out the tests, we consider the following Cauchy problem for stationary heat equation in , i.e.
.85) To compute the numerical solutions, we discretize the differential equation in (4.85) on a uniform grid of size N × M. Since the grid is uniform the number of grid points M in the y− directions is uniquely determined by N .We consider a standard accuracy of O(dx 2 ) in the finite difference approximation.Let T i, j be the discrete approximation to T (x i , y j ).The finite difference approximation for the equation on each interior grid points (x i , y j ) is given by where We cannot implement the degenerate boundary condition, directly, since c(x) tends to zero near the boundary.This means that we can get unbounded solutions.An alternative boundary condition that allows for similar solutions is considered.We choose the second order derivative with a standard accuracy of O(dx 2 ).The boundary i = N is treated similarly.This allows for functions that behave like a straight line near the boundary.This type of solution also satisfies the original boundary conditions (4.86).
Similarly, for the boundaries corresponding to y = 0, and y = a one can obtain different equations by taking into the consideration the boundary condition given.For the case of Dirichlet boundary conditions, we obtain the discretization as follows: where d j i is the Dirichlet data given at y = 0 or y = a.Furthermore, the discretization of the Neumann boundary conditions is given by We also observe that the rate of convergence is quite slow.The graphs for both methods are identical.The results are displayed in Fig. 4.
Example 4.2 In our second test example, we solve K η δ = ξ δ using the Landweber iteration.The parameter value chosen are α = 1.0, the noise levels are β = 4 • 10 1 and 12 • 10 1 , and τ = 1.1.The exact data and noisy data are displayed in Fig. 5.We see that the noise level is quite significant in both cases.
For the noise level β = 4 • 10 1 , we have seen that the solution obtained is of good quality.This shows that acceptable solutions can be found with realistic noise levels.We also displayed the errors ||η δ m − η|| γ, 1 , and residuals ||ξ δ − K η δ m || γ, 0 .We have observed that minimum norm solution and the smallest total error is reached after 871 iterations.Similarly, the residual is monotonically decreasing, and the discrepancy principle is achieved after 41 iterations.The results are displayed in Fig. 6.In the case of the noise level β = 12 • 10 1 , we can see that the solution obtained is also of good quality.We also presented the errors ||η δ m − η|| γ, 1 , and residuals ||ξ δ − K η δ m || γ, 0 .We have found out that the smallest total error is achieved after 474 iterations.The residual curve is also monotonically decreasing, and the discrepancy principle is satisfied after 21 iterations.The numerical results are displayed in Fig. 7.In both cases, we see that the errors decrease up-to a certain number of iterations and then it starts to increase.This is known as semi-convergences effect.From the residuals m || γ, 0 .The results indicates that the Landweber iteration converges for α 2. We also observed that the rate of converges is faster for smaller values of α as compare to larger one.If we compare with the theoretical result (3.81) we have an indication that the norm of the operator K is equal to 1.By looking at Lemma 3.4 this allows us to compute λ 1 .

Conclusions
In this paper, we consider steady state heat conduction in a thin plate.The thin plate connects two cylindrical containers.Its purpose is to fix their relative positions and at the same time to measure the temperature on the inner cylinder.We derive a two dimensional mathematical model, and use it to approximate the heat conduction in the thin plate.Since the plate has sharp edges on the sides, the resulting equation we obtained is a degenerate elliptic equation.The application of the Alternating iterative algorithm for solving a degenerate elliptic equation is new for this study.
To find the temperature on the interior part from the exterior measurements, we formulate the problem as a Cauchy problem for stationary heat equation.Furthermore, we describe two boundary value problems (1.4) and (1.5), and showed that they are well-posed.Similarly, we We have used here α = 1.9 (black curve), α = 2.0 (dashdot curve) and α = 2.1 (dashed curve) analysed the convergence of the Dirichlet-Neumann alternating algorithm for solving the Cauchy problem (1.3).Moreover, we have reformulated the Cauchy problem as an operator equation, with a compact operator, and applied the Landweber iteration method to solve the problem.We have shown that the Discrepancy principle works for this case and we have a regularization method.
From our numerical experiments, we see that it is possible to obtained good quality solutions, even in the presence of large noise.We also see that in all cases the residual norms are monotonically decreasing.This means that the Discrepancy principle can be applied in practice.We also compute errors during the Landweber iterations.The results seem to indicate that the Discrepancy principle underestimate the optimal number of iterations.We also see the semi-convergence effect in the error curves.From the numerical computation, we have seen that ill-posedness depends strongly on the geometry of the domain [7].Hence for this application to work the inner and outer cylinder has to be relatively close.
In future work, instead of Dirichlet-Neumann Alternating algorithm and the corresponding operator equation, we will consider using Dirichlet-Robin boundary conditions in the alternating algorithm.We will also consider methods with faster convergences rates such as the Conjugate gradient method.In order to implement the Conjugate gradient method we need to reformulate the alternating algorithm as a Landweber iteration for a specific operator equation.In this context, for this to work, we need to develop the scalar product that allows for the computation of the adjoint operator and also appropriate norms.With the adjoint operator in place we can implement the Conjugate gradient method.To improve on the stability of the solutions we will consider using Tikhonov regularization method.
Funding Open access funding provided by Linköping University.and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Fig.1The description of the engineering situation we are studying with two cylinders connected by a thin plate A. The plate is shown on the right .65) Substituting (3.65) in (3.63), we obtain (3.67)    Substituting (3.67) on the second term of the right-hand side in (3.61), we get the results.
.58) This proves the first formula in Lemma 2.5.The second formula can be proven similarly.123