Transient thermal mixed boundary value problems in the half-space

The Wiener-Hopf and Cagniard-de Hoop techniques are employed in order to solve a range of transient thermal mixed boundary value problems on the half-space. The thermal field is determined via a rapidly convergent integral, which can be evaluated straightforwardly and quickly on a desktop PC.


Introduction
Traditionally, great interest has been shown in determining the disturbances that are generated when loads are applied on the surface of a half-space. Lamb [1] obtained the exact solution when an impulsive, concentrated load is applied along a line of the free surface of an isotropic linear elastic medium. de Hoop reappraised this problem [2], modifying the method originally devised by Cagniard [3], [4], leading to the now wellknown Cagniard-de Hoop (CdH) technique. This method has been used widely since, allowing exact solutions to be obtained for a wide range of transient elasticity problems. The method can also be useful in order to render solutions into integral forms that are rapidly convergent when calculated numerically.
Transient thermoelastic half-space problems were considered by Danilovskaya [5], Boley and Tolins [6] and Achenbach [7] but in these problems the forcing was such that the CdH technique was not required. The extension of these problems to inhomogeneous media was considered by Baczynski [8] and Parnell [9]. A purely thermal, transient problem that employed the CdH method was solved in [10]. The thermoelastic Lamb problem was studied by Nayfeh and Nemat-Nasser [11] who used generalized thermoelasticity in order to retain a finite thermal wave speed, employing the CdH technique to determine the solution.
All of the above problems are of fundamental importance in an array of applications where a number of alternative boundary conditions on the surface can arise. What appears to be rather lacking in the literature however are studies of transient problems with mixed boundary conditions, where in the context of the thermal problems the condition takes the form, e.g.
T (0, y, t) = f (y, t) for y > 0, ∂T ∂x (0, y, t) = g(y, t) for y < 0, (1.1) where T (x, y, t) is the temperature field, f (y, t) and g(y, t) are two specified functions, t is time and with reference to Fig. 1 x and y are Cartesian coordinates. The half-space resides in x ≥ 0 and y runs parallel to the surface, which is defined by x = 0. Generally such problems lead to the propagation of a thermal disturbance into the half-space. Indeed, such thermal front problems are of importance in a number of applications including defect sizing [12], transient thermography [13], solar cell manufacturing [14] and thermal insulation [15]. Caflisch and Keller [16], Levine [17] and Satapathy and Sahoo [18] studied front propagation in the thermal context with mixed boundary conditions but in the context of steady problems with applications in quenching. Kozlov et al. [19] considered a transient half-space problem with mixed boundary conditions and made progress by using cylindrical coordinates due to the special form of the boundary condition chosen.
Mixed boundary conditions are generally difficult to handle even in steady problems and analytical or semi-analytical solutions are frequently only possible by the application of the Wiener-Hopf method [20]. This method exploits the analyticity properties of functions in order to yield an explicit or approximate solution in the Fourier transform domain. Contour integration then yields the solution in the physical domain.
Here we shall consider a rather general mixed boundary value problem in the context of thermal front propagation and determine solutions using the Wiener-Hopf method and Cagniard-de Hoop technique. This problem of mixed boundary conditions of the form (1.1) is of particular interest in analyzing the field close to the location of the change in boundary condition type, i.e. x = y = 0 in (1.1).
We obtain a solution in single integral form by using a deformation of the Laplace contour in a similar manner to the Cagniard-de Hoop method. Although it appears that we cannot obtain an explicit solution, the solution determined can be evaluated rapidly on a desktop PC and therefore it is of great utility due to its general form and its ability to circumvent a direct numerical simulation of the problem. Although similar problems, involving a discontinuous temperature boundary condition have been considered in the building insulation literature, see e.g. Claesson and Hegentoft [21] and Hegentoft and Claesson [22] to the authors' knowledge it does not appear that the solution we provide has been written down anywhere in the literature before now.
In this paper we shall first set out the problem description in section 2 before determining the solution in the transform domain in section 3. In section 4 we describe how we deform the Laplace contour onto a steepest descent path, in the manner of the Cagniard-de Hoop technique in order to obtain a solution in terms of a single integral along the deformed contour path, with an integrand that decays exponentially. In section 5 we illustrate the efficacy of the scheme by determining the solution for a number of different boundary conditions, with validation provided by finite element solutions.

Problem description
Assume that the problem under consideration is two-dimensional, being independent of z and define the two dimensional half-space domain D = {(x, y) : 0 ≤ x < ∞, −∞ < y < ∞}. We seek solutions to the anisotropic heat equation: where k and kℓ are the thermal conductivities (ℓ > 0) in the x and y directions respectively, c V is the specific heat at constant volume, ρ is the mass density, t is time and T = T (x, y, t) is the temperature field. We can combine the constants as κ = k/(ρc V ), the thermal diffusivity. It is convenient to non-dimensionalise the governing equation, using coordinates with a "hat" and scale the y coordinate to remove the anisotropy coefficient. Write (x,ŷ,T ,t) = (x/x * , y/y * , (T − T * )/T * , t/t * ), where and T * is the reference temperature in Kelvin. Upon doing so and "dropping hats" we find where ∇ 2 = ∂ 2 /∂x 2 + ∂ 2 /∂y 2 . We wish to solve (2.3) on the (scaled) domain D with boundary ∂D = ∂D − ∪ ∂D + as illustrated in Fig. 1. We consider homogeneous initial conditions of the form and boundary conditions of the form (1.1) but simplify by removing the y dependence, i.e.
where T 0 and T ′ 0 are real constants and f 0 (t) and g 0 (t) are piecewise continuous functions of time.
We therefore have a mixed boundary value problem, which in general are not straightforward to solve even in the steady context so that the time dependence adds an additional element of complexity. Furthermore we allow for the fact that we could have a step change at t = 0 on x = 0, leading to a propagating discontinuity front in the half-space.
In order to determine T it is convenient to introduce an alternative problem (for convergence issues as will be shown), giving rise to a different temperature distribution T , depending on a small parameter ǫ and where T converges to T as ǫ → 0. This problem is described as follows T (x ≥ 0, y, t = 0) = 0 (2.8) As is easily seen, we recover the solution to our original problem by taking the limit as ǫ tends to zero: T (x, y, t).

Solution in the transform domains via the Wiener-Hopf technique
Define the Laplace transform in time for any function φ(x, y, t) by and hence applying this to the governing scaled equations (2.7)-(2.10) we have and boundary conditions becomẽ Although s ∈ C the set of complex numbers, for the sake of the analysis to follow we can assume it to be real and positive. This allows us to scale the (x, y) variables to simplify the governing equations. The derivation goes through retaining explicit dependence on s in the governing equation, but the algebra becomes rather heavy and tedious and does not render any greater understanding of the problem; both approaches lead to the same result. As such we will rescale x and y in order to eliminate s from the governing equation. Thus define and therefore we obtain where Next define the Fourier transform in y 0 as and define Θ + and Θ − as so that Θ = Θ − + Θ + . Applying the Fourier transform to the governing (Laplace transformed) equation (3.6) we find that and to the boundary conditions, we find that With reference to Fig. 2, we note that Θ + is analytic on Ω + = {α ∈ C, ℑ(α) > −ǫ}, Θ − is analytic on Ω − = {α ∈ C, ℑ(α) < ǫ} and Θ + , Θ − and Θ are analytic on the strip S = Ω − ∩ Ω + . The superscript + and − notation thus indicates analyticity in the domains Ω + and Ω − respectively. The solution of (3.11) is where the branch of the square root function in the exponent is chosen so that its real part satisfies ℜ(α 2 + 1) 1/2 > 0, ensuring that the solution decays as x 0 → ∞. For conciseness let us introduce Φ + and Ψ − , analytic on Ω + and Ω − respectively, as Imposing the boundary conditions (3.12)-(3.13), employing (3.15) and eliminating A 1 between the two resulting equations, we arrive at The kernel here K(α) = (α 2 + 1) 1/2 is easily factorized as K(α) = K − /K + where Multiplying (3.16) by K + we obtain where we have indicated that we wish to determine a sum factorization of the function S. One can employ the pole removal method [23] in order to show quite straightforwardly that and Referring to (3.16) we can therefore define a function E(α) such that on Ω − \S, and therefore is analytic on the whole α-complex plane. It can be shown that as |α| → ∞, as |α| → ∞. This together with the analyticity of E(α) and the extended Liouville theorem (see for example [20]), implies that E(α) is constant. However, we also know that E(α) → 0 as |α| → ∞ and α ∈ Ω + . We can then conclude that E(α) = 0 everywhere. Given this, we therefore have from (3.26) From the original expressions for A 1 determined from the boundary conditions we can show that where (3.17), (3.23), (3.24) and (3.22) have all been used. Referring to (3.14), the solution in transform space is therefore Formally inverting the Fourier transform and using (3.5) gives the Laplace transformed solution asT where

Semi-analytical inversion via a Cagniard-de Hoop approach
Motivated by the Cagniard-de Hoop technique, let us introduce polar coordinates r and θ related to x and y in the usual manner, i.e. x = r cos θ, y = r sin θ where θ ∈ [−π/2, π/2] and introduce the parameter β via the expression Inverting for α we therefore determine the two paths B + and B − in the right and left halves of the α-plane In the α-plane, with β ∈ [1, ∞) these paths start at α = −i sin θ and move off to infinity either in the upper half-plane (θ ∈ (−π/2, 0), see Fig. 3) or the lower half-plane (θ ∈ (0, π/2), see Fig. 4).
Since B ± are steepest descent paths for the integrals, the idea is to deform the integrals (3.32) from the real line onto these to aid convergence. In classical Cagniard de Hoop problems this frequently permits one to write the α integral in the form of a Laplace transform of a function that is independent of s and thus we can determine the inverse transform immediately, thus rendering explicit solutions. Here we are not so fortunate, the function will not be independent of s but nevertheless we are able to make significant progress due to the fact that the inverse Laplace transform integral can be determined analytically in many important cases as we shall see shortly. This leaves the solution in a single integral form that is rapidly convergent.
At this point note that  Referring to (3.32) and Fig. 3, we see that the integrand of I 1 has a pole at α = −iǫ and a branch point at α = i. Deforming the contour from the real axis to the B ± contour in the α-plane we do not cross any singularities and hence we find that On B ± , we have α = α ± = −iβ sin θ ± β 2 − 1 cos θ and so we have dα = (−i sin θ ± β cos θ(β 2 − 1) −1/2 )dβ. (4.6) Noting that β ∈ [1, ∞) is a parametrization of the paths B ± we can rewrite I ± 1 as follows, taking the limit as ǫ → 0, Therefore from (4.5) we determine the form where As an aside, we note by deforming into the lower half-plane that (4.11)

Evaluation of I 2
Referring to (3.32) and Fig. 3, we see that the integrand of I 2 has a simple pole at α = iǫ and a branch point at α = i. Deforming the contour from the real axis to the B ± contour in the α-plane we cross the simple pole and pick up its residue. Accounting for this contribution we find that where and in particular, Noting that β ∈ [1, ∞) is a parametrization of the paths B ± we can rewrite I ± 2 as follows, taking the limit as ǫ → 0, Therefore from (4.12) and (4.14) we determine the form

An expression for T (r, θ, t)
We show in Appendix A that where G is a real-valued function. As such, using this together with (3.31), (4.8), (4.16) and (4.17) we find the following expression forT = lim ǫ→0T Finally we recover T by taking the inverse Laplace Transform, where T 1 (r, β, t) = 1 2iπ  This case follows entirely analogously to the negative θ scenario, the difference here being that the paths B ± reside in the lower half of the complex α-plane and as such the pole that leads to a contribution to the integral is at α = −iǫ. We find that As an aside, for θ ∈ (0, π/2], analogously to the derivation of (4.11) we have ∞ 1 F (β, θ) dβ = 0 (4.23) Therefore using, (4.11), (4.17) and (4.23) we have where H(θ) is the Heaviside step function.

A summary of the solution
Combining the results from sections 4.1 and 4.2, we can write down the solution for all values of θ, + H(θ)T 0 T 1 (r, cos θ, t) − (1 − H(θ))T ′ 0 T 2 (r, cos θ, t). Both integrals in (4.25) and the additional term have a discontinuity at θ = 0 and as such this form is not particularly "clean". We are able to improve upon this form, using (4.24) to generate Each term in this expression is now continuous across θ = 0. We shall discuss this aspect further in the context of specific examples in the next section.
5 Some specific boundary conditions 5.1 Perfect insulator on y < 0 Let us now assume that T ′ 0 = 0 so that we have a perfect insulator on y < 0. We shall consider a variety of temperature profiles for y > 0. The solution is therefore obtained by setting T ′ 0 = 0 in (4.25) or equivalently (4.26). As such only T 1 enters the analysis.

Step temperature change
Take the simplest form, f 0 (t) = 1 so thatf 0 (s) = 1/s and we need to determine T 1 defined in (4.20). It transpires that it is convenient to differentiate the expression for T 1 with respect to t, which enables the inverse Laplace integral to be evaluated analytically in this case We then integrate (definitely) with respect to t with a lower limit of t = 0, finding that Appealing to (4.25) we have as our solution noting that G(β, θ) is the real function defined in (A.12). Each term of (5.3) is easily computed numerically, and, in Fig. 5 we plot the resulting temperature profile on the horizontal axis against y running vertically, at time t = 0.02 for two values of x, x = 0.05 and x = 0.2. The circles are results taken from a finite element solution of the same problem in COMSOL and provide validation of the present semi-analytical scheme. Note that the temperature profile is continuous across the x-axis, i.e. θ = 0. Rather interestingly, both terms on the right hand side of (5.3) are discontinuous across θ = 0, as is shown in Fig. 6 but the two discontinuities compensate exactly to yield a continuous temperature profile. If instead of the form (4.25), we use (4.26), the solution is written as This expression is also straightforward to evaluate numerically and (obviously) gives the same results as those presented in Fig. 5, but this time, as shown in Fig. 7, both terms on the right hand side of (5.4) are continuous across θ = 0. Finally, in Fig. 8, we plot the two dimensional temperature contour profile on the (x, y) plane at t = 0.02 illustrating how the distribution spreads out from the upper half-plane.

Imperfect insulator on y < 0
Let us now consider the case when T ′ 0 = 0 and let us take f (t) = g(t) = 1 so that f 0 (s) =g 0 (s) = 1/s and now both T 1 and T 2 play a role. Once again it is convenient to differentiate with respect to t in order to evaluate the inverse Laplace transforms and subsequently integrating these expressions definitely with respect to t with a lower limit Either of the expressions (4.25) or (4.26) then recover the temperature profile. Both formulations are easily computed numerically and give rise to a continuous temperature profile, as seen in Fig. 9 where the profile is plotted at t = 0.02 for x = 0.05 and x = 0.2. The contour plot of the thermal field at t = 0.02 is given in Fig. 10. can be unsteady and in particular when f 0 (t) is a step and superposed general ramp up and down profile given by and illustrated in Fig. 11. Figure 11: Plot of f 0 (t) as defined by (5.7). A ramp up and down is superposed on a step change unit temperature profile.
The Laplace transform of f 0 (t) is        has already been dealt with in Section 5.1.1 and is thus given by (5.2). Moreover, T 1 (r, β, t) and T (4) 1 (r, β, t) have a very similar structure, so we only need to consider only the case of T  twice with respect to time, we obtain which, following the same reasoning as in Section 5.1.1 implies that 5.12) and by definite integration with respect to time, we obtain The functions T 1 (r, β, t) and T 1 (r, β, t) are obtained in the same manner. Using these expressions in the general formulation (4.25) or (4.26) enables the solution to be computed rather rapidly. In Fig. 12 we plot the resulting thermal field at the locations (x, y) = (0.1, 0.03) and (0.2, 0.03) in the cases when T 0 = 1 with T ′ 0 = 0 (left) and T 0 = 1 with T ′ 0 = 1 (right). We also plot the associated solutions determined previously where no ramp up and down is present in the boundary condition.