Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
Non-steady-state heat conduction in composite walls Bernard Deconinck, Beatrice Pelloni and Natalie E. Sheils Proc. R. Soc. A 2014 470, 20130605, published 12 March 2014
Supplementary data
"Data Supplement" http://rspa.royalsocietypublishing.org/content/suppl/2014/03/14/rs pa.2013.0605.DC1.html
References
This article cites 3 articles, 2 of which can be accessed free
http://rspa.royalsocietypublishing.org/content/470/2165/2013060 5.full.html#ref-list-1 Article cited in: http://rspa.royalsocietypublishing.org/content/470/2165/20130605.full.ht ml#related-urls
Subject collections
Articles on similar topics can be found in the following collections applied mathematics (373 articles) differential equations (114 articles) mathematical modelling (122 articles)
Email alerting service
Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here
To subscribe to Proc. R. Soc. A go to: http://rspa.royalsocietypublishing.org/subscriptions
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
Non-steady-state heat conduction in composite walls Bernard Deconinck1 , Beatrice Pelloni2 rspa.royalsocietypublishing.org
and Natalie E. Sheils1 1 Department of Applied Mathematics, University of Washington,
Research Cite this article: Deconinck B, Pelloni B, Sheils NE. 2014 Non-steady-state heat conduction in composite walls. Proc. R. Soc. A 470: 20130605. http://dx.doi.org/10.1098/rspa.2013.0605 Received: 9 September 2013 Accepted: 14 February 2014
Subject Areas: applied mathematics, differential equations, mathematical modelling Keywords: heat equation, interface problem, non-steady-state heat conduction, composite walls, Fokas method Author for correspondence: Natalie Sheils e-mail: [email protected]
Seattle, WA 98195-2420, USA 2 Department of Mathematics, University of Reading, Reading RG6 6AX, UK The problem of heat conduction in one-dimensional piecewise homogeneous composite materials is examined by providing an explicit solution of the one-dimensional heat equation in each domain. The location of the interfaces is known, but neither temperature nor heat flux is prescribed there. Instead, the physical assumptions of their continuity at the interfaces are the only conditions imposed. The problem of two semi-infinite domains and that of two finite-sized domains are examined in detail. We indicate also how to extend the solution method to the setting of one finite-sized domain surrounded on both sides by semi-infinite domains, and on that of three finite-sized domains.
1. Introduction The problem of heat conduction in a composite wall is a classical problem in design and construction. It is usual to restrict to the case of walls whose constitutive parts are in perfect thermal contact and have physical properties that are constant throughout the material and that are considered to be of infinite extent in the directions parallel to the wall. Further, we assume that temperature and heat flux do not vary in these directions. In that case, the mathematical model for heat conduction in each wall layer is given by [1, ch. 10] (j)
(j)
ut = κj uxx , and
Electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2013.0605 or via http://rspa.royalsocietypublishing.org.
(j)
x ∈ (aj , bj )
u(j) (x, t = 0) = u0 (x),
x ∈ (aj , bj ),
(1.1a) (1.1b)
where u(j) (x, t) denotes the temperature in the wall layer indexed by (j), κj > 0 is the heat-conduction coefficient of the jth layer (the inverse of its thermal diffusivity),
2014 The Author(s) Published by the Royal Society. All rights reserved.
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
x = aj is the left extent of the layer and x = bj is its right extent. The sub-indices denote derivatives with respect to the one-dimensional spatial variable x and the temporal variable t. The function
2. Two semi-infinite domains In this section, we consider the problem of heat flow through two walls of semi-infinite width or of two semi-infinite rods. We seek two functions: uL (x, t),
x ∈ (−∞, 0), t ≥ 0
and
uR (x, t),
x ∈ (0, ∞), t ≥ 0,
satisfying the equations uLt (x, t) = σL2 uLxx (x, t), and
2 R uR t (x, t) = σR uxx (x, t),
x ∈ (−∞, 0), t > 0
(2.1a)
x ∈ (0, ∞), t > 0,
(2.1b)
with the initial conditions being uL (x, 0) = uL0 (x), and
x ∈ (−∞, 0)
(2.2a)
(x, 0) = uR 0 (x),
x ∈ (0, ∞),
(2.2b)
lim uL (x, t) = γL ,
t≥0
(2.3a)
R
u
the asymptotic conditions being x→−∞
...................................................
(j)
of its associated heat flux κj ux are imposed across the interface between layers. In what follows, √ it is convenient to use the quantity σj , defined as the positive square root of κj : σj = κj . If the layer is either at the far left or far right of the wall, Dirichlet, Neumann, or Robin boundary conditions can be imposed on its far left or right boundary, respectively, corresponding to prescribing ‘outside’ temperature, heat flux or a combination of these. A derivation of the interface boundary conditions is found in [1, ch. 1]. It should be noted that the set-up presented in (1.1a) also applies to the case of one-dimensional rods in thermal contact. In this paper, we use the Fokas method [2–4] to provide explicit solution formulae for different heat transport interface problems of the type described above. We investigate problems in both finite and infinite domains, and we compare our method with classical solution approaches that can be found in the literature. Throughout, our emphasis is on non-steady-state solutions. Even for the simplest of the problems we consider (§3, two finite walls in thermal contact), the classical approach using separation of variables [1] can provide an explicit answer only implicitly. Indeed, the solution obtained in [1] depends on certain eigenvalues defined through a transcendental equation that can be solved only numerically. By contrast, the Fokas method produces an explicit solution formula involving only known quantities. For other problems we consider, to our knowledge no solution has been derived using classical methods and we believe the solution formulae presented here are new. The representation formulae for the solution can be evaluated numerically, hence the problem can be solved in practice using hybrid analytical–numerical approaches [5] or asymptotic approximations for them may be obtained using standard techniques [3]. The result of such a numerical calculation is shown at the end of §2. The problem of heat conduction through composite walls is discussed in many excellent texts (see for instance [1,6]). References to the treatment of specific problems are given in the sections below where these problems are investigated. In §2, we investigate the problem of two semiinfinite walls. Section 3 discusses the interface problem with two finite walls. Following that, we consider first the problem of one finite wall between two semi-infinite ones and the problem of three finite walls. Both of them are briefly sketched in §4 and full solutions are presented in the electronic supplementary material.
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
(j)
u0 (x) is the prescribed initial condition of the system. The continuity of the temperature u(j) and
2
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
and
3 R
(2.3b)
t>0
(2.4a)
and the continuity interface conditions being uL (0, t) = uR (0, t), and σL2 uLx (0, t) = σR2 uR x (0, t),
t > 0.
(2.4b)
The sub- and super-indices L and R denote the left and right rods, respectively. A special case of this problem is discussed in [1, ch. 10], but only for a specific initial condition. Further, for the problem treated there, both limx→∞ uR (x, t) and limx→−∞ uL (x, t), are assumed to be zero. This assumption is made for mathematical convenience and no physical reason exists to impose it. If constant (in time) limit values are assumed, a simple translation allows one of the limit values to be equated to zero, but not both. As no great advantage is obtained by assuming a zero limit using our approach, we make the more general assumption (2.3). We define v L (x, t) = uL (x, t) − γ L and v R (x, t) = uR (x, t) − γ R . Then, v L (x, t) and v R (x, t) satisfy L (x, t), vtL (x, t) = σL2 vxx
x ∈ (−∞, 0), t ≥ 0,
(2.5a)
R vtR (x, t) = σR2 vxx (x, t),
x ∈ (0, ∞), t ≥ 0,
(2.5b)
lim v (x, t) = 0,
t ≥ 0,
(2.5c)
lim v R (x, t) = 0,
t ≥ 0,
(2.5d)
v L (0, t) + γ L = v R (0, t) + γ R ,
t≥0
(2.5e)
σL2 vxL (0, t) = σR2 vxR (0, t),
t ≥ 0.
(2.5f )
L
x→−∞ x→∞
and
At this point, we start by following the standard steps in the application of the Fokas method [2–4]. We begin with the so-called ‘local relations’ [2] (e−ikx+(σL k) t v L (x, t))t = (σL2 e−ikx+(σL k) t (vxL (x, t) + ikv L (x, t)))x
(2.6a)
(e−ikx+(σR k) t v R (x, t))t = (σR2 e−ikx+(σR k) t (vxR (x, t) + ikv R (x, t)))x .
(2.6b)
2
2
and 2
2
These relations are a one-parameter family obtained by rewriting (2.5a,b). Applying Green’s formula [7] in the strip (−∞, 0) × (0, t) in the left-half plane (figure 1), we find t 0
⇒ +
(e−ikx+(σL k) s v L (x, t))s − (σL2 e−ikx+(σL k) s (vxL (x, s) + ikv L (x, s)))x dx ds = 0 2
0 −∞
0 −∞
t
0
e−ikx v0L (x) dx − 2
2
0 −∞
e−ikx+(σL k) t v L (x, t) dx 2
σL2 e(σL k) s (vxL (0, s) + ikv L (0, s)) ds = 0.
(2.7)
Let C+ = {z ∈ C : Im(z) ≥ 0}. Similarly, let C− = {z ∈ C : Im(z) ≤ 0}. As |x| can become arbitrarily large, we require k ∈ C+ in (2.7) in order to guarantee that the first two integrals are well defined. Let D = {k ∈ C : Re(k2 ) < 0} = D+ ∪ D− . The region D is shown in figure 2.
...................................................
t≥0
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
lim u (x, t) = γR ,
x→∞
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
t
4 ...................................................
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
DL
DR
x
Figure 1. Domains for the application of Green’s formula for v L (x, t) and v R (x, t). (Online version in colour.)
Im(k)
D+
Re(k)
D–
Figure 2. The domains D+ and D− for the heat equation. (Online version in colour.)
For k ∈ C, we define the following transforms: g0 (ω, t) =
t 0
eωs v L (0, s) ds =
t 0
eωs (v R (0, s) + γ R − γ L ) ds
t (γ R − γ L )(eωt − 1) + eωs v R (0, s) ds, ω 0 t t 2 σ g1 (ω, t) = eωs vxL (0, s) ds = R2 eωs vxR (0, s) ds, σL 0 0 0 0 e−ikx v L (x, t) dx, vˆ0L (k) = e−ikx v0L (x) dx vˆ L (k, t) = =
−∞
and
vˆ R (k, t) =
∞ 0
e−ikx v R (x, t) dx,
−∞
vˆ0R (k) =
∞ 0
e−ikx v0R (x) dx.
Using these definitions, the global relation (2.7) is rewritten as 2
vˆ0L (k) − e(σL k) t vˆ L (k, t) + ikσL2 g0 ((σL k)2 , t) + σL2 g1 ((σL k)2 , t) = 0,
k ∈ C+ .
(2.8)
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
As the dispersion relation ωL (k) = (σL k)2 is invariant under k → −k, so are g0 ((σL k)2 , t) and g1 ((σL k)2 , t). Thus, we can supplement (2.8) with its evaluation at −k, namely
This relation is valid on k ∈ C− . Using Green’s formula on (0, ∞) × (0, t) (figure 1), the global relation for v R (x, t) is 2 (γ L − γ R )(e(σR k) t − 1) R (σR k)2 t R 2 2 vˆ (k, t) − ikσR g0 ((σR k) , t) + vˆ0 (k) − e (σR k)2 − σL2 g1 ((σR k)2 , t) = 0,
(2.10)
valid in k ∈ C− . As above, using the invariance of ωR (k) = (σR k)2 , g0 ((σR k)2 , t) and g1 ((σR k)2 , t) under k → −k, we supplement (2.10) with 2 (γ L − γ R )(e(σR k) t − 1) R (σR k)2 t R 2 2 vˆ (−k, t) + ikσR g0 ((σR k) , t) + − σL2 g1 ((σR k)2 , t) = 0, vˆ0 (−k) − e (σR k)2 (2.11) for k ∈ C+ . Inverting the Fourier transforms in (2.8), we have 1 v (x, t) = 2π L
∞ −∞
2 eikx−(σL k) t vˆ0L (k) dk +
σL2 2π
∞ −∞
2
eikx−(σL k) t (ikg0 ((σL k)2 , t) + g1 ((σL k)2 , t)) dk (2.12)
for x ∈ (−∞, 0) and t > 0. The integrand of the second integral in (2.12) is entire and decays as k → ∞ for k ∈ C− \ D− . Using the analyticity of the integrand and applying Jordan’s lemma [7], we can replace the contour of integration of the second integral by − ∂D− : v L (x, t) =
1 2π
∞ −∞
σL2 2π
2
eikx−(σL k) t vˆ0L (k) dk −
2
∂D−
eikx−(σL k) t (ikg0 ((σL k)2 , t) + g1 ((σL k)2 , t)) dk. (2.13)
Proceeding similarly on the right, starting from (2.10), we have 1 ∞ ikx−(σR k)2 t R e vˆ0 (k) dk v R (x, t) = 2π −∞ 2 1 ∞ ikx−(σR k)2 t (γ L − γ R )(e(σR k) t − 1) 2 2 − e ikσR g0 ((σR k) , t) + 2π −∞ (σR k)2 + σL2 g1 ((σR k)2 , t) dk, ⎞⎞ ⎛ ⎛ ∞ x γL − γR ⎝ 2 ⎠⎠ + 1 eikx−(σR k) t vˆ0R (k) dk = 1 − erf ⎝ 2 2π 2 −∞ 2 σ t −
1 2π
∞ −∞
R
2
eikx−(σR k) t (ikσR2 g0 ((σR k)2 , t) + σL2 g1 ((σR k)2 , t)) dk,
⎞⎞ ⎛ ⎛ ∞ x γL − γR ⎝ 2 ⎠⎠ + 1 = eikx−(σR k) t vˆ0R (k) dk 1 − erf ⎝ 2 2π 2 −∞ 2 σ t 1 − 2π
R
2
∂D+
eikx−(σR k) t (ikσR2 g0 ((σR k)2 , t) + σL2 g1 ((σR k)2 , t)) dk,
(2.14)
√ z for x ∈ (0, ∞) and t > 0. Here, erf(·) denotes the error function: erf(z) = (2/ π ) 0 exp(−y2 ) dy. To obtain the second equality above, we integrated the terms that are explicit.
...................................................
(2.9)
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
2
vˆ0L (−k) − e(σL k) t vˆ L (−k, t) − ikσL2 g0 ((σL k)2 , t) + σL2 g1 ((σL k)2 , t) = 0.
5
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
g1 ((σL k)2 , t) =
L R kσL σ v ˆ (−k, t) − σ v ˆ , t R L σR σL2 (σL + σR )
2 kσL i(γL − γR )(1 − e(σL k) t ) − σR vˆ0L (−k) + , + σL vˆ0R σR kσL2 (σL + σR ) 1
2
e(σL k)
t
(2.16)
valid for k ∈ C− . These expressions are substituted into (2.13) and (2.14). This results in expressions for v L (x, t) and v R (x, t) that appear to depend on v L and v R themselves. We examine the contribution of the terms involving vˆ L and vˆ R . Starting with (2.13), we obtain for v L (x, t) the following expression: ⎞⎞ ⎛ ⎛ ∞ x σR (γ R − γ L ) ⎝ 2 L ⎠⎠ + 1 eikx−(σL k) t vˆ0L (k) dk v (x, t) = 1 + erf ⎝ σL + σR 2π 2 −∞ 2 σ t L
+ − +
∂D−
∂D−
∂D−
σR − σL 2 eikx−(σL k) t vˆ0L (−k) dk 2π (σL + σR )
σL ikx−(σL k)2 t R kσL dk vˆ0 e π (σL + σR ) σR
kσL σL − σR σL eikx vˆ L (−k, t) dk + eikx vˆ R , t dk 2π (σL + σR ) σR ∂D− π (σL + σR )
(2.17)
for x ∈ (−∞, 0), t > 0. The first four terms depend only on known functions. To examine the second-to-last term, we note that the integrand is analytic for all k ∈ C− and that vˆ L (−k, t) decays for k → ∞ for k ∈ C− . Thus, by Jordan’s lemma, the integral of exp(ikx)vˆ L (−k, t) along a closed, bounded curve in C− vanishes. In particular, we consider the closed curve L− = L∂D− ∪ L− C, − : |k| = C} (figure 3). = {k ∈ D where L∂D− = ∂D− ∩ {k : |k| < C} and L− C As the integral along L− C vanishes for large C, the fourth integral on the right-hand side of (2.17) must vanish as the contour L∂D− becomes ∂D− as C → ∞. The uniform decay of vˆ L (−k, t) for large k is exactly the condition required for the integral to vanish, using Jordan’s lemma. For the final integral in equation (2.17), we use the fact that uˆ R (kσL /σR , t) is analytic and bounded for k ∈ C− . Using the same argument as above, the fifth integral in (2.17) vanishes and we have an explicit representation for v L (x, t) in terms of initial conditions ⎞⎞ ⎛ ⎛ ∞ R − γ L) x (γ σ 2 R ⎠⎠ + 1 ⎝1 + erf ⎝ v L (x, t) = eikx−(σL k) t vˆ0L (k) dk σL + σR 2π 2 −∞ 2 σ t L
+ −
∂D−
∂D−
σR − σL 2 eikx−(σL k) t vˆ0L (−k) dk 2π (σL + σR )
kσL σL 2 eikx−(σL k) t vˆ0R dk. π (σL + σR ) σR
(2.18)
...................................................
and
6
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
Expressions (2.13) and (2.14) for v L (x, t) and v R (x, t) depend on the unknown functions g0 and g1 , evaluated at different arguments. These functions need to be expressed in terms of known quantities. To obtain a system of two equations for the two unknown functions, we use (2.9) and (2.10) for g0 ((σL k)2 , t) and g1 ((σL k)2 , t). This requires the transformation k → −σL k/σR in (2.10). The − sign is required to ensure that both equations are valid on C− , allowing for their simultaneous solution. We find
kσL i 2 e(σL k) t vˆ L (−k, t) + vˆ R ,t g0 ((σL k)2 , t) = kσL (σL + σR ) σR
2 kσL (γL − γR )(1 − e(σL k) t ) + (2.15) − vˆ0L (−k) − vˆ0R σR (σL k)2 (σL + σR )
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
Im(k)
7
Re(k)
L–∂D LC–
Figure 3. The contour L− is shown as a dashed line. An application of Cauchy’s integral theorem [7] using this contour allows the elimination of the contribution of vˆ L (−k, t) from integral expression (2.17). Similarly, the contour L+ is shown as a solid line and application of Cauchy’s integral theorem using this contour allows the elimination of the contribution of uˆR (−k, t) from integral expression (2.19). (Online version in colour.)
To find an explicit expression for v R (x, t), we need to evaluate g0 and g1 at different arguments, also ensuring that the expressions are valid for k ∈ C+ \ D+ . From (2.15) and (2.16), we find g0 ((σR k)2 , t) =
kσR −i 2 e(σR k) t vˆ L , t + vˆ R (−k, t) kσR (σL + σR ) σL
2 kσR (γ L − γ R )(1 − e(σR k) t ) − vˆ0R (−k) + − vˆ0L σL k2 σR (σL + σR )
and g1 ((σR k)2 , t) =
(σR k)2 t L kσR R σ e v ˆ , t − σ v ˆ (−k, t) + σL vˆ0R (−k) R L σL σL2 (σL + σR )
2 kσR i(γ L − γ R )(1 − e(σR k) t ) − . − σR vˆ0L σL kσL2 (σL + σR ) 1
Substituting these into equation (2.14), we obtain ⎞⎞ ⎛ ⎛ ∞ L − γ R) x (γ σ 2 L ⎠⎠ + 1 ⎝1 − erf ⎝ eikx−(σR k) t vˆ0R (k) dk v R (x, t) = 2π σL + σR 2 −∞ 2 σ t R
+ + +
∂D+
∂D+
∂D+
σR − σL 2 eikx−(σR k) t vˆ0R (−k) dk 2π (σL + σR )
kσR σR 2 eikx−(σR k) t vˆ0L , t dk π (σL + σR ) σL
σL − σR σR ikx R ikx L kσR e vˆ (−k, t) dk − e vˆ , t dk 2π (σL + σR ) σL ∂D+ π (σL + σR )
(2.19)
for x ∈ (0, ∞), t > 0. As before, everything about the first three integrals is known. To compute the fourth integral, we proceed as we did before for v L (x, t) and eliminate integrals that decay in
...................................................
+ L∂D
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
L+C
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
the regions over which we are integrating. The final solution is ⎞⎞ ⎛ ⎛ ∞ L − γ R) (γ σ x 2 L ⎠⎠ + 1 ⎝1 − erf ⎝ v R (x, t) = eikx−(σR k) t vˆ0R (k) dk σL + σR 2π −∞ 2 σ 2t
∂D+
σR − σL 2 eikx−(σR k) t v0R (−k) dk 2π (σL + σR )
kσR σR 2 eikx−(σR k) t v0L dk. π (σL + σR ) σL
(2.20)
Returning to the original variables, we have the following proposition which determines uR and uL fully explicitly in terms of the given initial conditions and the prescribed boundary conditions as |x| → ∞. Proposition 2.1. The solution of the heat transfer problem (2.1)–(2.4) is given by ⎞⎞ ⎛ ⎛ ∞ R − γ L) x (γ σ 2 R ⎠⎠ + 1 ⎝1 − erf ⎝ eikx−(σL k) t vˆ0L (k) dk uL (x, t) = γ L + σL + σR 2π −∞ 2 σ 2t L
+ − and
∂D−
∂D−
σR − σL 2 eikx−(σL k) t vˆ0L (−k) dk 2π (σL + σR )
kσL σL 2 eikx−(σL k) t vˆ0R dk π (σL + σR ) σR
(2.21a)
⎞⎞ ⎛ ⎛ ∞ x σL (γ L − γ R ) ⎝ 2 ⎠⎠ + 1 eikx−(σR k) t vˆ0R (k) dk 1 − erf ⎝ u (x, t) = γ + σL + σR 2π 2 −∞ 2 σ t R
R
R
+ +
∂D+
∂D+
σR − σL 2 eikx−(σR k) t v0R (−k) dk 2π (σL + σR )
kσR σR 2 eikx−(σR k) t v0L dk. π (σL + σR ) σL
(2.21b)
Remarks. — The use of the discrete symmetries of the dispersion relation is an important aspect of the Fokas method [2–4]. When solving the heat equation in a single medium, the only discrete symmetry required is k → −k, which was used here as well to obtain (2.9) and (2.11). Owing to the two media, there are two dispersion relations in this problem: ω1 = (σL k)2 and ω2 = (σR k)2 . The collection of both dispersion relations {ω1 , ω2 } retains the discrete symmetry k → −k but admits two additional ones, namely: k → (σR /σL )k and k → (σL /σR )k, which transform the two dispersion relations to each other. All non-trivial discrete symmetries of {ω1 , ω2 } are needed to derive the final solution representation, and indeed they are used, for example, to obtain relations (2.15) and (2.16). — With σL = σR and γ L = γ R = 0, solution formulae (2.21) in their proper x-domain of definition reduce to the solution of the whole line problem as given in [3]. — Classical approaches to the problem presented in this section can be found in the literature, for the case γL = 0 = γR . For instance, for one special pair of initial conditions, a solution is presented in [1]. No explicit solution formulae using classical methods with general initial conditions exist to our knowledge. At best, one is left with having to find the solution of an equation involving inverse Laplace transforms, where the unknowns are embedded within these inverse transforms. — The steady-state solution to (2.1) with initial conditions, which decay sufficiently fast to the boundary values (2.3) at ±∞, is easily obtained by letting t → ∞ in (2.21). This gives limt→∞ uR (x, t) = limt→∞ uL (x, t) = (γ L σL + γ R σR )/(σL + σR ). This is the weighted
...................................................
+
∂D+
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
R
+
8
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
u(x,t)
9
–0.01
0.01 x
0
u(x, t) 0.02
–0.01
0.01 0
0.01
x u(x, t)
5 × 10–7 0.01
x
0.01
–0.01 t
0
Figure 4. Results for solutions (2.18) and (2.20) with uL0 (x) = x 2 e(25) x , uR0 (x) = x 2 e−(30) x and σL = 0.02, σR = 0.06, γ L = γ R = 0, t ∈ [0, 0.02] using the hybrid analytical–numerical method of [5]. (Online version in colour.) 2
2
average of the boundary conditions at infinity with weights given by σL and σR . To our knowledge, this is a new result. It is consistent with the steady-state limit (γ L + γ R )/2 for the whole-line problem with initial conditions that limit to different values γ L and γ R as x → ±∞. This result is easily obtained from the solution of the heat equation defined on the whole line using the Fokas method, but it can also be observed by employing piecewise-constant initial data in the classical Green’s function solution, as described in Theorem 4-1 on p. 171 (and comments thereafter) of [8]. It should be emphasized that the steady-state problem for (2.1)–(2.4) or even for the heat equation defined on the whole
...................................................
t
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
5 × 10–7
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
line with different boundary conditions at +∞ and −∞ is ill posed in the sense that the steady solution cannot satisfy the boundary conditions.
2
and 2 −αR x , uR 0 (x) = x e 2
with αL = 25 and αR = 30 (figure 4). The Fourier transforms of these initial conditions may be computed explicitly. We choose σL = 0.02 and σR = 0.06. The initial conditions are chosen so as to satisfy interface boundary conditions (2.4) at t = 0. The results clearly illustrate the discontinuity in the first derivative of the temperature at the interface x = 0.
3. Two finite domains Next, we consider the problem of heat conduction through two walls of finite width (or of two finite rods) with Robin boundary conditions. We seek two functions: uL (x, t),
x ∈ (−a, 0), t ≥ 0 and uR (x, t),
x ∈ (0, b), t ≥ 0,
satisfying the equations uLt (x, t) = σL2 uLxx (x, t),
−a < x < 0, t > 0
(3.1a)
0 < x < b, t > 0,
(3.1b)
uL (x, 0) = uL0 (x),
−a < x < 0
(3.2a)
uR (x, 0) = uR 0 (x),
0 < x < b,
(3.2b)
and 2 R uR t (x, t) = σR uxx (x, t),
the initial conditions
and
the boundary conditions f L (t) = α1 uL (−a, t) + α2 uLx (−a, t),
t>0
(3.3a)
and f R (t) = α3 uR (b, t) + α4 uR x (b, t),
t>0
(3.3b)
and the continuity conditions uL (0, t) = uR (0, t),
t>0
(3.4a)
and σL2 uLx (0, t) = σR2 uR x (0, t),
t > 0,
(3.4b)
as illustrated in figure 5, where a > 0, b > 0 and αi , 1 ≤ i ≤ 4 constant. If α1 = α3 = 0, then Neumann boundary conditions are prescribed, whereas if α2 = α4 = 0 then Dirichlet conditions are given. As before, we have the local relations (e−ikx+(σL k) t uL (x, t))t = (σL2 e−ikx+(σL k) t (uLx (x, t) + ikuL (x, t)))x
(3.5a)
R (e−ikx+(σR k) t uR (x, t))t = (σR2 e−ikx+(σR k) t (uR x (x, t) + iku (x, t)))x .
(3.5b)
2
2
and 2
2
...................................................
uL0 (x) = x2 eαL x
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
Using a slight variation on the method presented in [5], one can compute solutions (2.21) numerically with specified initial conditions. We plot solutions for the case of vanishing boundary conditions (γ L = γ R = 0) with
10
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
utL = sL2 uLxx
utR = sR2 uRxx
11
x
Figure 5. The heat equation for two finite rods.
We define the time transforms of the initial and boundary data and the spatial transforms of u for k ∈ C as follows: uˆ L0 (k) = uˆ R 0 (k) = fˆL (ω, t) = hL1 (ω, t) = hR 1 (ω, t) = g1 (ω, t) = g0 (ω, t) =
and
0 −a
b 0
t
0
t
0
t
0
t
0
t
0
e−ikx uL0 (x) dx,
e−ikx uR 0 (x) dx, eωs fL (s) ds,
uˆ L (k, t) = uˆ R (k, t) =
fˆR (ω, t) =
eωs uLx (−a, s) ds, eωs uR x (b, s) ds, eωs uLx (0, s) ds = eωs uL (0, s) ds =
t 0
0 −a
b
σR2 σL2 t
t 0
e−ikx uR (x, t) dx,
0
eωs fR (s) ds,
hL0 (ω, t) = hR 0 (ω, t) =
e−ikx uL (x, t) dx,
t
t
eωs uL (−a, s) ds,
0
eωs uR (b, s) ds,
0
eωs uR x (0, s) ds
eωs uR (0, s) ds.
0
Using Green’s theorem on the domains [−a, 0] × [0, t] and [0, b] × [0, t], respectively, we have the global relations 2
e(σL k) t uˆ L (k, t) = σL2 (g1 ((σL k)2 , t) + ikg0 ((σL k)2 , t)) − eika σL2 (hL1 ((σL k)2 , t) + ikhL0 ((σL k)2 , t)) + uˆ L0 (k)
(3.6a)
and 2 R 2 2 2 e(σR k) t uˆ R (k, t) = e−ikb σR2 (hR 1 ((σR k) , t) + ikh0 ((σR k) ), t) − σL g1 ((σR k) , t) 2
− ikσR2 g0 ((σR k)2 , t) + uˆ R 0 (k).
(3.6b)
Both equations are valid for all k ∈ C, in contrast to (2.8) and (2.10). Using the invariance of ωL (k) = (σL k)2 and ωR (k) = (σR k)2 under k → −k, we obtain 2
e(σL k) t uˆ L (−k, t) = σL2 (g1 ((σL k)2 , t) − ikg0 ((σL k)2 , t)) − e−ika σL2 (hL1 ((σL k)2 , t) − ikhL0 ((σL k)2 , t)) + uˆ L0 (−k)
(3.7a)
and 2
2 R 2 e(σR k) t uˆ R (−k, t) = eikb σR2 (hR 1 ((σR k) , t) − ikh0 ((σR k) ), t)
− σL2 g1 ((σR k)2 , t) + ikσR2 g0 ((σR k)2 , t) + uˆ R 0 (−k).
(3.7b)
...................................................
b
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
0
–a
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
Im(k)
12 ...................................................
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
D0+
Re(k)
D0–
Figure 6. Deformation of the contours in figure 2 away from k = 0. (Online version in colour.)
Inverting the Fourier transform in (3.6a), 1 ∞ ikx−(σL k)2 t 2 uL (x, t) = e σL (g1 ((σL k)2 , t) + ikg0 ((σL k)2 , t)) dk 2π −∞ 1 ∞ ik(a+x)−(σL k)2 t 2 L − e σL (h1 ((σL k)2 , t) + ikhL0 ((σL k)2 , t)) dk 2π −∞ 1 ∞ ikx−(σL k)2 t L uˆ 0 (k) dk. + e 2π −∞
(3.8)
The integrand of the first integral is entire and decays as k → ∞ for k ∈ C− \ D− . The second integral has an integrand that is entire and decays as k → ∞ for k ∈ C+ \ D+ . It is convenient to deform both contours away from k = 0 to avoid singularities in the integrands that become apparent in what follows. Initially, these singularities are removable, as the integrands are entire. Writing integrals of sums as sums of integrals, the singularities may cease to be removable. With the deformations away from k = 0, the apparent singularities are no cause for concern. In other − − words, we deform D+ to D+ 0 and D to D0 as shown in figure 6. Thus, −1 2 uL (x, t) = eikx−(σL k) t σL2 (g1 ((σL k)2 , t) + ikg0 ((σL k)2 , t)) dk 2π ∂D−0 1 2 − eik(a+x)−(σL k) t σL2 (hL1 ((σL k)2 , t) + ikhL0 ((σL k)2 , t)) dk + 2π ∂D0 1 ∞ ikx−(σL k)2 t L uˆ 0 (k) dk. + e (3.9) 2π −∞ To obtain the solution on the right, we apply the inverse Fourier transform to (3.6b) 1 ∞ ik(x−b)−(σR k)2 t 2 R 2 e σR (h1 ((σR k)2 , t) + ikhR uR (x, t) = 0 ((σR k) , t)) dk 2π −∞ 1 ∞ ikx−(σR k)2 t − e (ikσR2 g0 ((σR k)2 , t) + σL2 g1 ((σR k)2 , t)) dk 2π −∞ 1 ∞ ikx−(σR k)2 t R uˆ 0 (k) dk. e + 2π −∞
(3.10)
The integrand of the first integral is entire and decays as k → ∞ for k ∈ C− \ D− . The second integral has an integrand that is entire and decays as k → ∞ for k ∈ C+ \ D+ . We deform the
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
contours as above to obtain 2
∂D− 0
−
1 2π
+
1 2π
2 R 2 eik(x−b)−(σR k) t σR2 (hR 1 ((σR k) , t) + ikh0 ((σR k) , t)) dk
2
∂D+ 0
∞
−∞
eikx−(σR k) t (ikσR2 g0 ((σR k)2 , t) + σL2 g1 ((σR k)2 , t)) dk 2
eikx−(σR k) t uˆ R 0 (k) dk.
(3.11)
Taking the time transform of the boundary conditions results in fˆ L (ω, s) = α1 hL0 (ω, t) + R α2 hL1 (ω, t) and fˆ R (ω, s) = α3 hR 0 (ω, t) + α4 h1 (ω, t). These two equations together with (3.7a,b) are four equations to be solved for the four L R unknowns hL0 (ω, t), hR 0 (ω, t), h1 (ω, t) and h1 (ω, t). The resulting expressions are substituted in (3.9) and (3.11). Although we could solve this problem in its full generality, we restrict to the case of Dirichlet boundary conditions (α2 = α4 = 0), to simplify the already cumbersome formulae below. Then, hL0 (ω, t) and hR 0 (ω, t) are determined and we solve two equations for two unknowns. System (3.7a,b) is not solvable for hL1 (ω, t) and hR 1 (ω, t) if L (k) = 0, where L (k) = π (σL (e2iak + 1)(e2ibkσL /σR − 1) + σR (e2iak − 1)(e2ibkσL /σR + 1))
bkσL 2iak 2ibkσL /σR + σR tan(ak) . + 1)(e + 1) σL tan = iπ (e σR It is easily seen that all values of k satisfying this (including k = 0) are on the real line. Thus on the contours, the equations are solved without problem, resulting in the expressions below. As before, the right-hand sides of these expressions involve uˆ L (k, t) and uˆ L (k, t), evaluated at a variety of arguments. All terms with such dependence are written out explicitly below. Terms that depend on only known quantities are contained in KL and KR , the expressions for which are given later. uL (x, t) = KL + + + + + + + +
∂D− 0
eikx (σL + σR ) + eik(x+2bσL /σR ) (σR − σL ) L uˆ (k, t) dk 2 L (k)
∂D− 0
−eik(x+2a) (σL + σR ) + eik(x+2a+2bσL /σR ) (σL − σR ) L uˆ (−k, t) dk 2 L (k)
σL eik(x+2a+2bσL /σR ) R kσL uˆ , t dk L (k) σR
−σL eik(x+2a) R −kσL uˆ , t dk L (k) σR
∂D+ 0
eik(x+2a) (σR − σL ) + eik(x+2a+2bσL /σR ) (σL + σR ) L uˆ (k, t) dk 2 L (k)
∂D− 0
∂D− 0
∂D+ 0
∂D+ 0
∂D+ 0
−eik(x+2a) (σL + σR ) + eik(x+2a+2bσL /σR ) (σL − σR ) L uˆ (−k, t) dk 2 L (k)
σL eik(x+2a+2bσL /σR ) R kσL ˆu , t dk L (k) σR
−σL eik(x+2a) R −kσL , t dk uˆ L (k) σR
(3.12)
...................................................
−1 2π
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
uR (x, t) =
13
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
and
14
+ + + + + +
kσR ,t σL
dk
−kσR , t dk σL
∂D+ 0
−σR eik(x+2aσR /σL ) L uˆ R (k)
∂D+ 0
eik(2b+x) (σL − σR ) + eik(x+2b+2aσR /σL ) (σL + σR ) R uˆ (k, t) dk 2 R (k)
∂D− 0
eikx (σR − σL ) − eik(x+2aσR /σL ) (σL + σR ) R uˆ (−k, t) dk 2 R (k)
σR eikx L kσR −σR eik(x+2aσR /σL ) L −kσR uˆ uˆ , t dk + , t dk R (k) σL R (k) σL ∂D− 0
∂D− 0
eikx (σL + σR ) − eik(x+2aσR /σL ) (σL − σR ) R uˆ (k, t) dk 2 R (k)
∂D− 0
eikx (σR − σL ) − eik(x+2aσR /σL ) (σL + σR ) R uˆ (−k, t) dk, 2 R (k)
∂D+ 0
(3.13)
where R (k) = L (kσR /σL ). The integrands written explicitly in (3.12) and (3.13) decay in the regions around whose boundaries they are integrated. Thus, using Jordan’s lemma and Cauchy’s theorem, these integrals are shown to vanish. Thus, the final solution is given by KL and KR . Proposition 3.1. The solution of heat transfer problem (3.1)–(3.4) is given by uL (x, t) = KL =
∞ −∞
+ + + + +
+
2
∂D− 0
∂D− 0
∂D− 0
−eikx−(σL k) t (σL + σR ) + eik(x+2bσL /σR )−(σL k) t (σL − σR ) L uˆ 0 (k) dk 2 L (k)
2
2
2
2
∂D− 0
eik(x+2a)−(σL k) t (σL + σR ) + eik(x+2a+2bσL /σR )−(σL k) t (σR − σL ) L uˆ 0 (−k) dk 2 L (k)
2 −σL eik(x+2a+2bσL /σR )−(σL k) t R kσL uˆ 0 dk L (k) σR
2 σL eik(x+2a)−(σL k) t R −kσL uˆ 0 dk L (k) σR
∂D+ 0
ikσL2 eik(x+a)−(σL k) t (σL + σR ) − ikσL2 eik(x+a+2bσL /σR )−(σL k) t (σL − σR ) L fˆ ((σL k)2 , t) dk α1 L (k)
∂D+ 0
−2ikσL2 σR eik(x+2a+bσL /σR )−(σL k) t (1 + σL σR ) R fˆ ((σL k)2 , t) dk α3 L (k)
∂D+ 0
−eik(x+2a)−(σL k) t (σR − σL ) − eik(x+2a+2bσL /σR )−(σL k) t (σL + σR ) L uˆ 0 (k) dk 2 L (k)
2
∂D− 0
∂D− 0
2
2
+
−2ikσL2 σR eik(x+2a+bσL /σR )−(σL k) t R fˆ ((σL k)2 , t) dk α3 L (k)
ikσL2 eik(x+a)−(σL k) t (σL + σR ) − ikσL2 eik(x+a+2bσL /σR )−(σL k) t (σL − σR ) L fˆ ((σL k)2 , t) dk α1 L (k)
+
2
eikx−(σL k) t uˆ L0 (k) dk +
2
2
2
2
...................................................
∂D+ 0
eikx σR L uˆ R (k)
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
uR (x, t) = KR +
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
+
∂D+ 0
∂D+ 0
∂D+ 0
2
15
(3.14)
for −a < x < 0, and for 0 < x < b uR (x, t) = KR =
∞ −∞
+ + + +
+ + +
2
∂D− 0
∂D− 0
∂D− 0
−eik(x+2b)−(σR k) t (σL − σR ) − eik(x+2b+2aσR /σL )−(σR k) t (σL + σR ) R uˆ 0 (k) dk 2 R (k)
∂D− 0
eikx−(σR k) t (σL − σR ) + eik(x+2aσR /σL )−(σR k) t (σL + σR ) R uˆ 0 (−k) dk 2 R (k)
∂D+ 0
2ikσL σR2 eik(x+aσR /σL )−(σR k) t L fˆ ((σR k)2 , t) dk α1 R (k)
2
∂D− 0
2
2
2
2
2
2
+
2ikσL σR2 eik(x+aσR /σL )−(σR k) t L fˆ ((σR k)2 , t) dk α1 R (k)
−ikσR2 eik(x+b)−(σR k) t (σL − σR ) − ikσR2 eik(x+b+2aσR /σL )−(σR k) t R fˆ ((σR k)2 , t) dk α3 R (k)
2 2 −σR eikx−(σR k) t L kσR σR eik(x+2aσR /σL )−(σR k) t L −kσR uˆ 0 uˆ 0 dk + dk R (k) σL R (k) σL ∂D− 0
+
2
eikx−(σR k) t uˆ R 0 (k) dk +
∂D+ 0
−ikσR eik(x+b)−(σR k) t (σL − σR ) − ikσR2 eik(x+b+2aσR /σL )−(σR k) t (σL + σR ) R fˆ ((σR k)2 , t) dk α3 R (k)
2 2 −σR eikx−(σR k) t L kσR σR eik(x+2aσR /σL )−(σR k) t L −kσR uˆ 0 uˆ 0 dk + dk R (k) σL R (k) σL ∂D+ 0
∂D+ 0
−eikx−(σR k) t (σL + σR ) + eik(x+2aσR /σL )−(σR k) t (σR − σL ) R uˆ 0 (k) dk 2 R (k)
∂D+ 0
eikx−(σR k) t (σL − σR ) + eik(x+2aσR /σL )−(σR k) t (σL + σR ) R uˆ 0 (−k) dk. 2 R (k)
2
∂D+ 0
2
2
2
2
2
(3.15)
Remarks. — The solution of the problem posed in (3.1)–(3.4) may be obtained using the classical method of separation of variables and superposition [1]. The solutions uL (x, t) and uR (x, t) are given by a series of eigenfunctions with eigenvalues that satisfy a transcendental equation, closely related to the equation L (k) = 0. This series solution may be obtained + from proposition 3.1 by deforming the contours along ∂D− 0 and ∂D0 to the real line, including small semicircles around each root of either L (k) or R (k), depending on whether uL (x, t) or uR (x, t) is being calculated. Indeed, this is allowed because all integrands decay in the wedges between these contours and the real line, and the zeros of L (k) and R (k) occur only on the real line, as stated above. Careful calculation of all different contributions, following the examples in [2,3], shows that the contributions along the real line cancel, leaving only residue contributions from the small circles. Each residue contribution corresponds to a term in the classical series solution. It is not necessarily beneficial to leave the form of the solution in proposition 3.1 for the series
...................................................
+
2
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
+
eik(x+2a)−(σL k) t (σL + σR ) + eik(x+2a+2bσL /σR )−(σL k) t (σR − σL ) L uˆ 0 (−k) dk 2 L (k)
2 −σL eik(x+2a+2bσL /σR )−(σL k) t R kσL uˆ 0 dk L (k) σR
2 σL eik(x+2a)−(σL k) t R −kσL uˆ 0 dk, L (k) σR
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
σR2 (γ R − γL ) bσL2 + aσR2
x+
bγ L σL2 + aγ R σR2 bσL2 + aσR2
−a < x < 0
,
and uR (x, t) =
σL2 (γ R − γL ) bσL2
+ aσR2
x+
bγ L σL2 + aγ R σR2 bσL2 + aσR2
0 < x < b,
,
which is piecewise linear and continuous. — A more direct way to recover only the steady-state solution to (3.1) with limt→∞ f L (t) = f¯L and limt→∞ f R (t) = f¯R constant is to write the solution as the superposition of two parts: uL (x, t) = u¯ L (x) + uˇ L (x, t) and uR (x, t) = u¯ R (x) + uˇ R (x, t). The first parts u¯ L and u¯ R satisfy the boundary conditions as t → ∞ and the stationary heat equation. In other words, 0 = σL2 u¯ Lxx (x),
−a < x < 0,
0 = σR2 u¯ R xx (x),
0 < x < b,
f¯L = α1 u¯ L (−a) + α2 u¯ Lx (−a) f¯R = α3 u¯ R (b) + α4 u¯ R x (b).
and
A piecewise linear ansatz with the imposition of the interface conditions results in linear algebra for the unknown coefficients [1]. With the steady-state solution in hand, the second (time-dependent) parts uˇ L and uˇ R satisfy the initial conditions modified by the steady-state solution and the boundary conditions minus their value as t → ∞ uˇ Lt (x, t) = σL2 uˇ Lxx (x, t),
a < x < 0, t > 0,
2 ˇR uˇ R xx (x, t), t (x, t) = σR u
0 < x < b, t > 0,
ˇL
L
¯L
L
u (x, 0) = u (x, 0) − u (x),
−a < x < 0,
uˇ R (x, 0) = uR (x, 0) − u¯ R (x),
0 < x < b,
¯L
ˇL
f (t) − f = α1 u and
(−a, t) + α2 uˇ Lx (−a, t),
f R (t) − f¯R = α3 uˇ R (b, t) + α4 uˇ R x (b, t),
t>0
t > 0,
where, as usual, we impose continuity of temperature and heat flux at the interface x = 0. The dynamics of the solution is described by uˇ L and uˇ R , both of which decay to zero as t → ∞. Their explicit form is easily found using the method described in this section.
...................................................
uL (x, t) =
16
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
representation, as the latter depends on the roots of L (k) and R (k), which are not known explicitly. By contrast, the representation of proposition 3.1 depends on known quantities only and may be readily computed, using one’s favourite parameterization of + the contours ∂D− 0 and ∂D0 . — Similarly, the familiar piecewise linear steady-state solution of (3.1)–(3.4) with Dirichlet boundary conditions [1] can be observed from (3.14) and (3.15) by choosing initial conditions that decay appropriately and constant boundary conditions fL (t) = γ L and fR (t) = γ R . It is convenient to choose zero initial conditions, as the initial conditions do not affect the steady state. As above, the contours are deformed so that they are along the real line with semicircular paths around the zeros of L (k) and R (k), including − k = 0. As one of these deformations arises from D+ 0 while the other comes from D0 , the contributions along the real line cancel each other, while the semicircles add to give full residue contributions from the poles associated with the zeros of L (k). All such residues vanish as t → ∞, except at k = 0. It follows that the steady-state behaviour is determined by the residue at the origin. This results in
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
4. Other problems
17
In this section, we consider the heat equation defined on two semi-infinite rods enclosing a single rod of length 2a. We seek three functions: uL (x, t), x ∈ (−∞, −a),
uM (x, t), x ∈ (−a, b)
and uR (x, t), x ∈ (b, ∞),
with t ≥ 0, satisfying the equations
and
uLt (x, t) = σL2 uLxx (x, t),
−∞ < x < −a, t > 0,
(4.1a)
2 M uM t (x, t) = σM u xx (x, t),
−a < x < a, t > 0
(4.1b)
2 R uR t (x, t) = σR uxx (x, t),
a < x < ∞, t > 0,
(4.1c)
uLt (x, 0) = uL0 (x),
−∞ < x < −a,
(4.2a)
uM t (x, 0) = uM 0 (x),
−a < x < a
(4.2b)
with the initial conditions being
R uR t (x, 0) = u0 (x),
and
a < x < ∞,
(4.2c)
the asymptotic conditions being
and
lim uL (x, t) = 0, x→−∞ t
t≥0
(4.3a)
lim uR (x, t) = 0, x→∞ t
t≥0
(4.3b)
and the continuity conditions being uL (−a, t) = uM (−a, t), uM (a, t) = uR (a, t),
t ≥ 0, t ≥ 0,
2 M σL2 uLx (−a, t) = σM ux (−a, t),
and
2 M σM ux (a, t) = σR2 uR x (a, t),
(4.4a) (4.4b)
t≥0
t ≥ 0,
(4.4c) (4.4d)
as shown in figure 7. Asymptotic conditions (4.3) are not the most general possible but are used here to simplify calculations. To our knowledge, no aspect of this problem is addressed in the literature. To solve this problem, one would proceed with the following steps: (i) Write the local relations in each of the three domains: left, middle and right. (ii) Use Green’s theorem to define the three global relations, keeping track of where they are valid in the complex k plane. (iii) Solve the three global relations for uL (x, t), uM (x, t) and uR (x, t) by inverting the Fourier transforms. (iv) Using the k → −k symmetry of the dispersion relationships on the three original global relations, write down three more global relations. This uses the discrete symmetry of each individual dispersion relation. − (v) Deform the integrals of the solution expressions to D+ 0 or D0 as dictated by the region of analyticity of the integrand.
...................................................
(a) Infinite domain with three layers
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
With the basic principles of the method outlined in the previous two sections, problems with more layers, both finite and infinite, may be addressed. We state two additional problems below and include solutions for specific initial conditions. More complete solutions can be found in the electronic supplementary material.
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
utL = sL2 uLxx
utM = sM2 uMxx
utR = sR2 uRxx
18
x
Figure 7. The heat conduction problem for a single rod of length 2a between two semi-infinite rods.
(vi) Use the discrete symmetries of the dispersion relationship to the collection of global relations for g0 (ω, t), g1 (ω, t), h0 (ω, t) and h1 (ω, t). Care should be taken that relations can only be solved simultaneously provided their regions of validity overlap. At this stage, as in the previous two sections, the discrete symmetries from one dispersion relation to another come into play. (vii) Terms containing unknown functions are shown to be zero by examining the regions of analyticity and decay for the relevant integrands and the use of Jordan’s lemma. For brevity of the presentation, we will assume that uM (x, 0) = 0 = uR (x, 0). The solution to the non-restricted problem can be found in the electronic supplementary material. After defining the transforms uˆ L0 (k) = uˆ M (k) =
−a
e−ikx uL0 (x) dx,
−∞ a
e−ikx uM (x, t) dx,
−a
t
uˆ L (k, t) =
uˆ R (k, t) = 2 σM
t
−∞ ∞
e−ikx uL (x, t) dx, e−ikx uR (x, t) dx,
a
eωs uM x (−a, s) ds, σL2 0 t t h0 (ω, t) = eωs uL (−a, s) ds = eωs uM (−a, s) ds,
h1 (ω, t) =
0
eωs uLx (−a, s) ds =
−a
0
0
0
eωs uM x (a, s) ds =
σR2
t
eωs uR x (a, s) ds 2 σM 0 t t g0 (ω, t) = eωs uM (a, s) ds = eωs uR (a, s) ds,
g1 (ω, t) = and
t
0
0
we proceed as outlined above. The solution formulae are given in the following proposition: Proposition 4.1. The solution of the heat transfer problem (4.1)–(4.4) with uM (x, 0) = uR (x, 0) = 0 is uL (x, t) =
1 2π −
∞
1 2
−∞
2
eikx−(σL k) t uˆ L0 (k) dk
2
eik(x+2a)−(σL k) t ((σL + σM )(σM − σR ) L (k)
∂D− 0
+ e4ikaσL /σM (σL − σM )(σM + σR ))uˆ L0 (−k) dk, for −∞ < x < −a with L (k) = π ((σL − σM )(σM − σR ) + e4iakσL /σM (σL + σM )(σM + σR )), uM (x, t) = −σM (σM − σR )
2
∂D− 0
+ σM (σM + σR )
−kσM dk σL
2 eik(x+a−aσM /σL )−(σM k) t L kσM uˆ 0 dk, R (kσM /σR ) σL
eik(x+a+aσM /σL )−(σM k) t L uˆ 0 L (kσM /σL )
∂D+ 0
...................................................
a
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
–a
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
utL = sL2 uLxx
utM = sM2 uM xx
utR = sR2 uRxx
19
x c
Figure 8. The heat equation for three finite layers.
for −a < x < b with R (k) = π ((σL + σM )(σM + σR ) + e4iakσR /σM (σL − σM )(σM − σR )). Lastly, the expression for uR (x, t), valid for x > a, is uR (x, t) = 2σM σR
eik(x−a−aσR /σL +2aσR /σM )−(σR k) t L uˆ 0 R (k) 2
∂D+ 0
kσR σL
dk.
(b) Finite domain with three layers We consider the heat conduction problem in three rods of finite length as shown in figure 8. We seek three functions: uL (x, t),
x ∈ (−a, 0), t ≥ 0,
uM (x, t),
x ∈ (0, b), t ≥ 0 and uR (x, t),
x ∈ (b, c), t ≥ 0,
satisfying the equations uLt (x, t) = σL2 uLxx (x, t), M
u and
−a < x < 0, t > 0,
(4.5a)
0 < x < b, t > 0
(4.5b)
2 M t (x, t) = σM u xx (x, t),
2 R uR t (x, t) = σR uxx (x, t),
b < x < c, t > 0,
(4.5c)
uLt (x, 0) = uL0 (x),
−∞ < x < −a,
(4.6a)
uM t (x, 0) = uM 0 (x),
−a < x < a
(4.6b)
with the initial conditions being
R uR t (x, 0) = u0 (x),
and
a < x < ∞,
(4.6c)
the boundary conditions being f L (t) = α1 uL (−a, t) + α2 uLx (−a, t), and
f R (t) = α3 uR (c, t) + α4 uR x (c, t),
t>0
t>0
(4.7a) (4.7b)
and the continuity conditions being uL (0, t) = uM (0, t), M
R
u (b, t) = u (b, t),
t ≥ 0,
(4.8a)
t ≥ 0,
(4.8b)
2 M σL2 uLx (0, t) = σM ux (0, t),
and
2 M σM ux (b, t) = σR2 uR x (b, t),
t≥0
(4.8c)
t ≥ 0.
(4.8d)
The solution process is as before, following the steps outlined in the previous section. For simplicity, we assume Neumann boundary data (α1 = α3 = 0), zero boundary conditions (f L (t) = f R (t) = 0) and uM (x, 0) = 0 = uR (x, 0). The solution with uM (x, 0) = 0 and uR (x, 0) = 0 is given in the
...................................................
b
0
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
–a
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
electronic supplementary material. We define
gL1 (ω, t) =
b
e−ikx uL0 (x) dx,
uˆ L (k, t) =
e−ikx uM (x, t) dx,
0
t
eωs uLx (−a, s) ds,
2 σM
t
t
e−ikx uR (x, t) dx,
b
eωs uL (−a, s) ds,
0
eωs uM x (0, s) ds, σL2 0 t t hL0 (ω, t) = eωs uL (0, s) ds = eωs uM (0, s) ds, hL1 (ω, t) =
0
eωs uLx (0, s) ds =
gL0 (ω, t) =
e−ikx uL (x, t) dx,
c
uˆ R (k, t) =
0
t
−a
0
t
0
0
eωs uM x (b, s) ds =
σR2
t
eωs uR x (b, s) ds, 2 σM 0 t t gR eωs uM (b, s) ds = eωs uR (b, s) ds 0 (ω, t) = gR 1 (ω, t) =
0
hR 1 (ω, t) =
and
t
0
0
eωs uR x (c, s) ds,
hR 0 (ω, t) =
t
eωs uR (c, s) ds.
0
The solution is given by the following proposition: Proposition 4.2. The solution to the heat transfer problem (4.5)–(4.8) with α1 = α3 = 0, f L (t) = f R (t) = 0 and uM (x, 0) = uR (x, 0) = 0 is 1 ∞ ikx−(σL k)2 t L uˆ 0 (k) dk uL (x, t) = e 2π −∞ 2 eikx−(σL k) t 2ibkσL /σR e (σL − σM )(σM − σR ) + 2 L (k) ∂D− 0 + e2ikσL (c/σR +b/σM ) (σL + σM )(σM − σR ) + e2ickσL /σR (σL − σM )(σM + σR ) +e2ibkσL (1/σR +1/σM ) (σL + σM )(σM + σR ) uˆ L0 (k) dk +
2
∂D− 0
eik(x+2a)−(σL k) 2 L (k)
t
(e2ibkσL /σR (σL − σM )(σM − σR )
+ e2ikσL (c/σR +b/σM ) (σL + σM )(σM − σR ) + e2ickσL /σR (σL − σM )(σM + σR ) +e2ibkσL (1/σR +1/σM ) (σL + σM )(σM + σR ) uˆ L0 (−k) dk +
∂D+ 0
2 eik(x+2a)−(σL k) t 2ikσL (c/σR +b/σM ) e (σL − σM )(σM − σR ) 2 L (k)
+ e2ibkσL /σR (σL + σM )(σM − σR ) + e2ibkσL (1/σR +1/σM ) (σL − σM )(σM + σR ) +e2ickσL /σR (σL + σM )(σM + σR ) uˆ L0 (k) dk +
2
∂D+ 0
eik(x+2a)−(σL k) 2 L (k)
t
(e2ibkσL /σR (σL − σM )(σM − σR )
+ e2ikσL (c/σR +b/σM ) (σL + σM )(σM − σR ) + e2ickσL /σR (σL − σM )(σM + σR ) +e2ibkσL (1/σR +1/σM ) (σL + σM )(σM + σR ) )uˆ L0 (−k) dk,
...................................................
uˆ M (k) =
−a
20
0
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
uˆ L0 (k) =
0
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
for −a < x < 0 with
21
+ e2ickσL /σR (σL − σM )(σM + σR ) + e2ik(a+bσL /σR +bσL /σM ) (σM − σL )(σM + σR ) + e2ibk(σL /σR +σL /σM ) (σL + σM )(σM + σR ) − e2ik(cσL /σR +a) (σL + σM )(σM + σR )). Next, uM (x, t) =
−eik(x+b)−(σM k) t σM 2ickσM /σR e (σM − σR ) − M (k) ∂D0 kσ M dk +e2ibkσM /σR (σM + σR ) uˆ L0 σL 2 −eik(x+b+2aσM /σL )−(σM k) t σM 2ickσM /σR e + (σM − σR ) M (k) ∂D− 0 −σ M dk +e2ibkσM /σR (σM + σR ) uˆ L0 σL 2 −eik(x+b)−(σM k) t σM 2ickσM /σR e + (σM − σR ) M (k) ∂D+ 0 kσ M dk +e2ibkσM /σR (σM + σR ) uˆ L0 σL 2 eik(x+b+2aσM /σL )−(σM k) t σM 2ickσM /σR e − (σM − σR ) M (k) ∂D+ 0 −kσ M dk, +e2ibkσM /σR (σM + σR ) uˆ L0 σL 2
for 0 < x < b, with M (k) = π (eik(aσM /σL +b+cσM /σR ) (σL − σM )(σM − σR ) + e2bikσM /σR (σM − σL )(σM − σR ) + e2ikσM (b/σR +a/σL ) (σL + σM )(σM − σR ) + e2ik(cσM /σR +b) (σL + σM )(σR − σM ) + e2ik(aσM /σL +bσM /σR +b) (σL − σM )(σM + σR ) + e2ickσM /σR (σM − σL )(σM + σR ) − e2ibk(σM /σR +1) (σL + σM )(σM + σR ) + e2ikσM (c/σR +a/σL ) (σL + σM )(σM + σR )) and uR (x, t) =
−2eik(x+b+bσR /σM )−(σR k) t (σM σR )uˆ L0 R (k) 2
∂D− 0
+ + +
∂D− 0
∂D+ 0
∂D+ 0
kσR σL
dk
2 −kσR −2 eik(x+b+bσR /σM +2aσR /σL )−(σR k) t (σM σR )uˆ L0 dk R (k) σL
2 kσR 2 eik(x+b+bσR /σM )−(σR k) t (σM σR )uˆ L0 dk R (k) σL
2 −kσR 2 eik(x+b+bσR /σM +2aσR /σL )−(σR k) t (σM σR )uˆ L0 dk, R (k) σL
...................................................
+ e2ikσL (c/σR +b/σM ) (σL + σM )(σM − σR ) + e2ik(bσL /σR +a) (σL + σM )(σR − σM )
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
L (k) = π (e2ibkσL /σR (σL − σM )(σM − σR ) + e2ik(cσL /σR +bσL /σM +a) (σM − σL )(σM − σR )
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
for b < x < c with
22
R (k) = π (e2ibk (σL − σM )(σM − σR ) + e2ik(aσR /σL +bσR /σM +c) (σM − σL )(σM − σR )
+ e2ibk(1+σR /σM ) (σL + σM )(σM + σR ) − e2ik(aσR /σL +c) (σL + σM )(σM + σR )). Acknowledgements. Any opinions, findings and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the funding sources. We thank the American Institute for Mathematics (AIM, Palo Alto, CA, USA) for its hospitality, while some of this work was being conducted. Funding statement. This work was generously supported by the National Science Foundation under grant no. NSF-DMS-1008001 (B.D., N.E.S.). N.E.S. also acknowledges support from the National Science Foundation under grant no. NSF-DGE-0718124.
References 1. Hahn D, Özisik M. 2012 Heat conduction, 3rd edn. Hoboken, NJ: John Wiley & Sons, Inc. 2. Deconinck B, Trogdon T, Vasan V. 2014 The method of Fokas for solving linear partial differential equations. SIAM Rev. 56, 159–186. (doi:10.1137/110821871) 3. Fokas AS. 2008 A unified approach to boundary value problems. CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 78. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM). 4. Fokas AS, Pelloni B. 2005 A transform method for linear evolution PDEs on a finite interval. IMA J. Appl. Math. 70, 564–587. (doi:10.1093/imamat/hxh047) 5. Flyer N, Fokas AS. 2008 A hybrid analytical-numerical method for solving evolution partial differential equations. I. The half-line. Proc. R. Soc. A 464, 1823–1849. (doi:10.1098/rspa. 2008.0041) 6. Carslaw HS, Jaeger JC. 1959 Conduction of heat in solids, 2nd edn. New York, NY: Oxford University Press. 7. Ablowitz MJ, Fokas AS. 2003 Complex variables: introduction and applications, 2nd edn, Cambridge Texts in Applied Mathematics. Cambridge, UK: Cambridge University Press. 8. Guenther RB, Lee JW. 1996 Partial differential equations of mathematical physics and integral equations. Mineola, NY: Dover Publications Inc. (Corrected reprint of the 1988 original.)
...................................................
+ e2ick (σL − σM )(σM + σR ) + e2ik(a/σL +b+bσR /σM ) (σM − σL )(σM + σR )
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
+ e2ik(bσR /σM +c) (σL + σM )(σM − σR ) + e2ik(aσR /σL +b) (σL + σM )(σR − σM )
Non-steady-state heat conduction in composite walls Bernard Deconinck, Beatrice Pelloni and Natalie E. Sheils Proc. R. Soc. A 2014 470, 20130605, published 12 March 2014
Supplementary data
"Data Supplement" http://rspa.royalsocietypublishing.org/content/suppl/2014/03/14/rs pa.2013.0605.DC1.html
References
This article cites 3 articles, 2 of which can be accessed free
http://rspa.royalsocietypublishing.org/content/470/2165/2013060 5.full.html#ref-list-1 Article cited in: http://rspa.royalsocietypublishing.org/content/470/2165/20130605.full.ht ml#related-urls
Subject collections
Articles on similar topics can be found in the following collections applied mathematics (373 articles) differential equations (114 articles) mathematical modelling (122 articles)
Email alerting service
Receive free email alerts when new articles cite this article - sign up in the box at the top right-hand corner of the article or click here
To subscribe to Proc. R. Soc. A go to: http://rspa.royalsocietypublishing.org/subscriptions
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
Non-steady-state heat conduction in composite walls Bernard Deconinck1 , Beatrice Pelloni2 rspa.royalsocietypublishing.org
and Natalie E. Sheils1 1 Department of Applied Mathematics, University of Washington,
Research Cite this article: Deconinck B, Pelloni B, Sheils NE. 2014 Non-steady-state heat conduction in composite walls. Proc. R. Soc. A 470: 20130605. http://dx.doi.org/10.1098/rspa.2013.0605 Received: 9 September 2013 Accepted: 14 February 2014
Subject Areas: applied mathematics, differential equations, mathematical modelling Keywords: heat equation, interface problem, non-steady-state heat conduction, composite walls, Fokas method Author for correspondence: Natalie Sheils e-mail: [email protected]
Seattle, WA 98195-2420, USA 2 Department of Mathematics, University of Reading, Reading RG6 6AX, UK The problem of heat conduction in one-dimensional piecewise homogeneous composite materials is examined by providing an explicit solution of the one-dimensional heat equation in each domain. The location of the interfaces is known, but neither temperature nor heat flux is prescribed there. Instead, the physical assumptions of their continuity at the interfaces are the only conditions imposed. The problem of two semi-infinite domains and that of two finite-sized domains are examined in detail. We indicate also how to extend the solution method to the setting of one finite-sized domain surrounded on both sides by semi-infinite domains, and on that of three finite-sized domains.
1. Introduction The problem of heat conduction in a composite wall is a classical problem in design and construction. It is usual to restrict to the case of walls whose constitutive parts are in perfect thermal contact and have physical properties that are constant throughout the material and that are considered to be of infinite extent in the directions parallel to the wall. Further, we assume that temperature and heat flux do not vary in these directions. In that case, the mathematical model for heat conduction in each wall layer is given by [1, ch. 10] (j)
(j)
ut = κj uxx , and
Electronic supplementary material is available at http://dx.doi.org/10.1098/rspa.2013.0605 or via http://rspa.royalsocietypublishing.org.
(j)
x ∈ (aj , bj )
u(j) (x, t = 0) = u0 (x),
x ∈ (aj , bj ),
(1.1a) (1.1b)
where u(j) (x, t) denotes the temperature in the wall layer indexed by (j), κj > 0 is the heat-conduction coefficient of the jth layer (the inverse of its thermal diffusivity),
2014 The Author(s) Published by the Royal Society. All rights reserved.
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
x = aj is the left extent of the layer and x = bj is its right extent. The sub-indices denote derivatives with respect to the one-dimensional spatial variable x and the temporal variable t. The function
2. Two semi-infinite domains In this section, we consider the problem of heat flow through two walls of semi-infinite width or of two semi-infinite rods. We seek two functions: uL (x, t),
x ∈ (−∞, 0), t ≥ 0
and
uR (x, t),
x ∈ (0, ∞), t ≥ 0,
satisfying the equations uLt (x, t) = σL2 uLxx (x, t), and
2 R uR t (x, t) = σR uxx (x, t),
x ∈ (−∞, 0), t > 0
(2.1a)
x ∈ (0, ∞), t > 0,
(2.1b)
with the initial conditions being uL (x, 0) = uL0 (x), and
x ∈ (−∞, 0)
(2.2a)
(x, 0) = uR 0 (x),
x ∈ (0, ∞),
(2.2b)
lim uL (x, t) = γL ,
t≥0
(2.3a)
R
u
the asymptotic conditions being x→−∞
...................................................
(j)
of its associated heat flux κj ux are imposed across the interface between layers. In what follows, √ it is convenient to use the quantity σj , defined as the positive square root of κj : σj = κj . If the layer is either at the far left or far right of the wall, Dirichlet, Neumann, or Robin boundary conditions can be imposed on its far left or right boundary, respectively, corresponding to prescribing ‘outside’ temperature, heat flux or a combination of these. A derivation of the interface boundary conditions is found in [1, ch. 1]. It should be noted that the set-up presented in (1.1a) also applies to the case of one-dimensional rods in thermal contact. In this paper, we use the Fokas method [2–4] to provide explicit solution formulae for different heat transport interface problems of the type described above. We investigate problems in both finite and infinite domains, and we compare our method with classical solution approaches that can be found in the literature. Throughout, our emphasis is on non-steady-state solutions. Even for the simplest of the problems we consider (§3, two finite walls in thermal contact), the classical approach using separation of variables [1] can provide an explicit answer only implicitly. Indeed, the solution obtained in [1] depends on certain eigenvalues defined through a transcendental equation that can be solved only numerically. By contrast, the Fokas method produces an explicit solution formula involving only known quantities. For other problems we consider, to our knowledge no solution has been derived using classical methods and we believe the solution formulae presented here are new. The representation formulae for the solution can be evaluated numerically, hence the problem can be solved in practice using hybrid analytical–numerical approaches [5] or asymptotic approximations for them may be obtained using standard techniques [3]. The result of such a numerical calculation is shown at the end of §2. The problem of heat conduction through composite walls is discussed in many excellent texts (see for instance [1,6]). References to the treatment of specific problems are given in the sections below where these problems are investigated. In §2, we investigate the problem of two semiinfinite walls. Section 3 discusses the interface problem with two finite walls. Following that, we consider first the problem of one finite wall between two semi-infinite ones and the problem of three finite walls. Both of them are briefly sketched in §4 and full solutions are presented in the electronic supplementary material.
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
(j)
u0 (x) is the prescribed initial condition of the system. The continuity of the temperature u(j) and
2
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
and
3 R
(2.3b)
t>0
(2.4a)
and the continuity interface conditions being uL (0, t) = uR (0, t), and σL2 uLx (0, t) = σR2 uR x (0, t),
t > 0.
(2.4b)
The sub- and super-indices L and R denote the left and right rods, respectively. A special case of this problem is discussed in [1, ch. 10], but only for a specific initial condition. Further, for the problem treated there, both limx→∞ uR (x, t) and limx→−∞ uL (x, t), are assumed to be zero. This assumption is made for mathematical convenience and no physical reason exists to impose it. If constant (in time) limit values are assumed, a simple translation allows one of the limit values to be equated to zero, but not both. As no great advantage is obtained by assuming a zero limit using our approach, we make the more general assumption (2.3). We define v L (x, t) = uL (x, t) − γ L and v R (x, t) = uR (x, t) − γ R . Then, v L (x, t) and v R (x, t) satisfy L (x, t), vtL (x, t) = σL2 vxx
x ∈ (−∞, 0), t ≥ 0,
(2.5a)
R vtR (x, t) = σR2 vxx (x, t),
x ∈ (0, ∞), t ≥ 0,
(2.5b)
lim v (x, t) = 0,
t ≥ 0,
(2.5c)
lim v R (x, t) = 0,
t ≥ 0,
(2.5d)
v L (0, t) + γ L = v R (0, t) + γ R ,
t≥0
(2.5e)
σL2 vxL (0, t) = σR2 vxR (0, t),
t ≥ 0.
(2.5f )
L
x→−∞ x→∞
and
At this point, we start by following the standard steps in the application of the Fokas method [2–4]. We begin with the so-called ‘local relations’ [2] (e−ikx+(σL k) t v L (x, t))t = (σL2 e−ikx+(σL k) t (vxL (x, t) + ikv L (x, t)))x
(2.6a)
(e−ikx+(σR k) t v R (x, t))t = (σR2 e−ikx+(σR k) t (vxR (x, t) + ikv R (x, t)))x .
(2.6b)
2
2
and 2
2
These relations are a one-parameter family obtained by rewriting (2.5a,b). Applying Green’s formula [7] in the strip (−∞, 0) × (0, t) in the left-half plane (figure 1), we find t 0
⇒ +
(e−ikx+(σL k) s v L (x, t))s − (σL2 e−ikx+(σL k) s (vxL (x, s) + ikv L (x, s)))x dx ds = 0 2
0 −∞
0 −∞
t
0
e−ikx v0L (x) dx − 2
2
0 −∞
e−ikx+(σL k) t v L (x, t) dx 2
σL2 e(σL k) s (vxL (0, s) + ikv L (0, s)) ds = 0.
(2.7)
Let C+ = {z ∈ C : Im(z) ≥ 0}. Similarly, let C− = {z ∈ C : Im(z) ≤ 0}. As |x| can become arbitrarily large, we require k ∈ C+ in (2.7) in order to guarantee that the first two integrals are well defined. Let D = {k ∈ C : Re(k2 ) < 0} = D+ ∪ D− . The region D is shown in figure 2.
...................................................
t≥0
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
lim u (x, t) = γR ,
x→∞
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
t
4 ...................................................
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
DL
DR
x
Figure 1. Domains for the application of Green’s formula for v L (x, t) and v R (x, t). (Online version in colour.)
Im(k)
D+
Re(k)
D–
Figure 2. The domains D+ and D− for the heat equation. (Online version in colour.)
For k ∈ C, we define the following transforms: g0 (ω, t) =
t 0
eωs v L (0, s) ds =
t 0
eωs (v R (0, s) + γ R − γ L ) ds
t (γ R − γ L )(eωt − 1) + eωs v R (0, s) ds, ω 0 t t 2 σ g1 (ω, t) = eωs vxL (0, s) ds = R2 eωs vxR (0, s) ds, σL 0 0 0 0 e−ikx v L (x, t) dx, vˆ0L (k) = e−ikx v0L (x) dx vˆ L (k, t) = =
−∞
and
vˆ R (k, t) =
∞ 0
e−ikx v R (x, t) dx,
−∞
vˆ0R (k) =
∞ 0
e−ikx v0R (x) dx.
Using these definitions, the global relation (2.7) is rewritten as 2
vˆ0L (k) − e(σL k) t vˆ L (k, t) + ikσL2 g0 ((σL k)2 , t) + σL2 g1 ((σL k)2 , t) = 0,
k ∈ C+ .
(2.8)
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
As the dispersion relation ωL (k) = (σL k)2 is invariant under k → −k, so are g0 ((σL k)2 , t) and g1 ((σL k)2 , t). Thus, we can supplement (2.8) with its evaluation at −k, namely
This relation is valid on k ∈ C− . Using Green’s formula on (0, ∞) × (0, t) (figure 1), the global relation for v R (x, t) is 2 (γ L − γ R )(e(σR k) t − 1) R (σR k)2 t R 2 2 vˆ (k, t) − ikσR g0 ((σR k) , t) + vˆ0 (k) − e (σR k)2 − σL2 g1 ((σR k)2 , t) = 0,
(2.10)
valid in k ∈ C− . As above, using the invariance of ωR (k) = (σR k)2 , g0 ((σR k)2 , t) and g1 ((σR k)2 , t) under k → −k, we supplement (2.10) with 2 (γ L − γ R )(e(σR k) t − 1) R (σR k)2 t R 2 2 vˆ (−k, t) + ikσR g0 ((σR k) , t) + − σL2 g1 ((σR k)2 , t) = 0, vˆ0 (−k) − e (σR k)2 (2.11) for k ∈ C+ . Inverting the Fourier transforms in (2.8), we have 1 v (x, t) = 2π L
∞ −∞
2 eikx−(σL k) t vˆ0L (k) dk +
σL2 2π
∞ −∞
2
eikx−(σL k) t (ikg0 ((σL k)2 , t) + g1 ((σL k)2 , t)) dk (2.12)
for x ∈ (−∞, 0) and t > 0. The integrand of the second integral in (2.12) is entire and decays as k → ∞ for k ∈ C− \ D− . Using the analyticity of the integrand and applying Jordan’s lemma [7], we can replace the contour of integration of the second integral by − ∂D− : v L (x, t) =
1 2π
∞ −∞
σL2 2π
2
eikx−(σL k) t vˆ0L (k) dk −
2
∂D−
eikx−(σL k) t (ikg0 ((σL k)2 , t) + g1 ((σL k)2 , t)) dk. (2.13)
Proceeding similarly on the right, starting from (2.10), we have 1 ∞ ikx−(σR k)2 t R e vˆ0 (k) dk v R (x, t) = 2π −∞ 2 1 ∞ ikx−(σR k)2 t (γ L − γ R )(e(σR k) t − 1) 2 2 − e ikσR g0 ((σR k) , t) + 2π −∞ (σR k)2 + σL2 g1 ((σR k)2 , t) dk, ⎞⎞ ⎛ ⎛ ∞ x γL − γR ⎝ 2 ⎠⎠ + 1 eikx−(σR k) t vˆ0R (k) dk = 1 − erf ⎝ 2 2π 2 −∞ 2 σ t −
1 2π
∞ −∞
R
2
eikx−(σR k) t (ikσR2 g0 ((σR k)2 , t) + σL2 g1 ((σR k)2 , t)) dk,
⎞⎞ ⎛ ⎛ ∞ x γL − γR ⎝ 2 ⎠⎠ + 1 = eikx−(σR k) t vˆ0R (k) dk 1 − erf ⎝ 2 2π 2 −∞ 2 σ t 1 − 2π
R
2
∂D+
eikx−(σR k) t (ikσR2 g0 ((σR k)2 , t) + σL2 g1 ((σR k)2 , t)) dk,
(2.14)
√ z for x ∈ (0, ∞) and t > 0. Here, erf(·) denotes the error function: erf(z) = (2/ π ) 0 exp(−y2 ) dy. To obtain the second equality above, we integrated the terms that are explicit.
...................................................
(2.9)
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
2
vˆ0L (−k) − e(σL k) t vˆ L (−k, t) − ikσL2 g0 ((σL k)2 , t) + σL2 g1 ((σL k)2 , t) = 0.
5
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
g1 ((σL k)2 , t) =
L R kσL σ v ˆ (−k, t) − σ v ˆ , t R L σR σL2 (σL + σR )
2 kσL i(γL − γR )(1 − e(σL k) t ) − σR vˆ0L (−k) + , + σL vˆ0R σR kσL2 (σL + σR ) 1
2
e(σL k)
t
(2.16)
valid for k ∈ C− . These expressions are substituted into (2.13) and (2.14). This results in expressions for v L (x, t) and v R (x, t) that appear to depend on v L and v R themselves. We examine the contribution of the terms involving vˆ L and vˆ R . Starting with (2.13), we obtain for v L (x, t) the following expression: ⎞⎞ ⎛ ⎛ ∞ x σR (γ R − γ L ) ⎝ 2 L ⎠⎠ + 1 eikx−(σL k) t vˆ0L (k) dk v (x, t) = 1 + erf ⎝ σL + σR 2π 2 −∞ 2 σ t L
+ − +
∂D−
∂D−
∂D−
σR − σL 2 eikx−(σL k) t vˆ0L (−k) dk 2π (σL + σR )
σL ikx−(σL k)2 t R kσL dk vˆ0 e π (σL + σR ) σR
kσL σL − σR σL eikx vˆ L (−k, t) dk + eikx vˆ R , t dk 2π (σL + σR ) σR ∂D− π (σL + σR )
(2.17)
for x ∈ (−∞, 0), t > 0. The first four terms depend only on known functions. To examine the second-to-last term, we note that the integrand is analytic for all k ∈ C− and that vˆ L (−k, t) decays for k → ∞ for k ∈ C− . Thus, by Jordan’s lemma, the integral of exp(ikx)vˆ L (−k, t) along a closed, bounded curve in C− vanishes. In particular, we consider the closed curve L− = L∂D− ∪ L− C, − : |k| = C} (figure 3). = {k ∈ D where L∂D− = ∂D− ∩ {k : |k| < C} and L− C As the integral along L− C vanishes for large C, the fourth integral on the right-hand side of (2.17) must vanish as the contour L∂D− becomes ∂D− as C → ∞. The uniform decay of vˆ L (−k, t) for large k is exactly the condition required for the integral to vanish, using Jordan’s lemma. For the final integral in equation (2.17), we use the fact that uˆ R (kσL /σR , t) is analytic and bounded for k ∈ C− . Using the same argument as above, the fifth integral in (2.17) vanishes and we have an explicit representation for v L (x, t) in terms of initial conditions ⎞⎞ ⎛ ⎛ ∞ R − γ L) x (γ σ 2 R ⎠⎠ + 1 ⎝1 + erf ⎝ v L (x, t) = eikx−(σL k) t vˆ0L (k) dk σL + σR 2π 2 −∞ 2 σ t L
+ −
∂D−
∂D−
σR − σL 2 eikx−(σL k) t vˆ0L (−k) dk 2π (σL + σR )
kσL σL 2 eikx−(σL k) t vˆ0R dk. π (σL + σR ) σR
(2.18)
...................................................
and
6
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
Expressions (2.13) and (2.14) for v L (x, t) and v R (x, t) depend on the unknown functions g0 and g1 , evaluated at different arguments. These functions need to be expressed in terms of known quantities. To obtain a system of two equations for the two unknown functions, we use (2.9) and (2.10) for g0 ((σL k)2 , t) and g1 ((σL k)2 , t). This requires the transformation k → −σL k/σR in (2.10). The − sign is required to ensure that both equations are valid on C− , allowing for their simultaneous solution. We find
kσL i 2 e(σL k) t vˆ L (−k, t) + vˆ R ,t g0 ((σL k)2 , t) = kσL (σL + σR ) σR
2 kσL (γL − γR )(1 − e(σL k) t ) + (2.15) − vˆ0L (−k) − vˆ0R σR (σL k)2 (σL + σR )
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
Im(k)
7
Re(k)
L–∂D LC–
Figure 3. The contour L− is shown as a dashed line. An application of Cauchy’s integral theorem [7] using this contour allows the elimination of the contribution of vˆ L (−k, t) from integral expression (2.17). Similarly, the contour L+ is shown as a solid line and application of Cauchy’s integral theorem using this contour allows the elimination of the contribution of uˆR (−k, t) from integral expression (2.19). (Online version in colour.)
To find an explicit expression for v R (x, t), we need to evaluate g0 and g1 at different arguments, also ensuring that the expressions are valid for k ∈ C+ \ D+ . From (2.15) and (2.16), we find g0 ((σR k)2 , t) =
kσR −i 2 e(σR k) t vˆ L , t + vˆ R (−k, t) kσR (σL + σR ) σL
2 kσR (γ L − γ R )(1 − e(σR k) t ) − vˆ0R (−k) + − vˆ0L σL k2 σR (σL + σR )
and g1 ((σR k)2 , t) =
(σR k)2 t L kσR R σ e v ˆ , t − σ v ˆ (−k, t) + σL vˆ0R (−k) R L σL σL2 (σL + σR )
2 kσR i(γ L − γ R )(1 − e(σR k) t ) − . − σR vˆ0L σL kσL2 (σL + σR ) 1
Substituting these into equation (2.14), we obtain ⎞⎞ ⎛ ⎛ ∞ L − γ R) x (γ σ 2 L ⎠⎠ + 1 ⎝1 − erf ⎝ eikx−(σR k) t vˆ0R (k) dk v R (x, t) = 2π σL + σR 2 −∞ 2 σ t R
+ + +
∂D+
∂D+
∂D+
σR − σL 2 eikx−(σR k) t vˆ0R (−k) dk 2π (σL + σR )
kσR σR 2 eikx−(σR k) t vˆ0L , t dk π (σL + σR ) σL
σL − σR σR ikx R ikx L kσR e vˆ (−k, t) dk − e vˆ , t dk 2π (σL + σR ) σL ∂D+ π (σL + σR )
(2.19)
for x ∈ (0, ∞), t > 0. As before, everything about the first three integrals is known. To compute the fourth integral, we proceed as we did before for v L (x, t) and eliminate integrals that decay in
...................................................
+ L∂D
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
L+C
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
the regions over which we are integrating. The final solution is ⎞⎞ ⎛ ⎛ ∞ L − γ R) (γ σ x 2 L ⎠⎠ + 1 ⎝1 − erf ⎝ v R (x, t) = eikx−(σR k) t vˆ0R (k) dk σL + σR 2π −∞ 2 σ 2t
∂D+
σR − σL 2 eikx−(σR k) t v0R (−k) dk 2π (σL + σR )
kσR σR 2 eikx−(σR k) t v0L dk. π (σL + σR ) σL
(2.20)
Returning to the original variables, we have the following proposition which determines uR and uL fully explicitly in terms of the given initial conditions and the prescribed boundary conditions as |x| → ∞. Proposition 2.1. The solution of the heat transfer problem (2.1)–(2.4) is given by ⎞⎞ ⎛ ⎛ ∞ R − γ L) x (γ σ 2 R ⎠⎠ + 1 ⎝1 − erf ⎝ eikx−(σL k) t vˆ0L (k) dk uL (x, t) = γ L + σL + σR 2π −∞ 2 σ 2t L
+ − and
∂D−
∂D−
σR − σL 2 eikx−(σL k) t vˆ0L (−k) dk 2π (σL + σR )
kσL σL 2 eikx−(σL k) t vˆ0R dk π (σL + σR ) σR
(2.21a)
⎞⎞ ⎛ ⎛ ∞ x σL (γ L − γ R ) ⎝ 2 ⎠⎠ + 1 eikx−(σR k) t vˆ0R (k) dk 1 − erf ⎝ u (x, t) = γ + σL + σR 2π 2 −∞ 2 σ t R
R
R
+ +
∂D+
∂D+
σR − σL 2 eikx−(σR k) t v0R (−k) dk 2π (σL + σR )
kσR σR 2 eikx−(σR k) t v0L dk. π (σL + σR ) σL
(2.21b)
Remarks. — The use of the discrete symmetries of the dispersion relation is an important aspect of the Fokas method [2–4]. When solving the heat equation in a single medium, the only discrete symmetry required is k → −k, which was used here as well to obtain (2.9) and (2.11). Owing to the two media, there are two dispersion relations in this problem: ω1 = (σL k)2 and ω2 = (σR k)2 . The collection of both dispersion relations {ω1 , ω2 } retains the discrete symmetry k → −k but admits two additional ones, namely: k → (σR /σL )k and k → (σL /σR )k, which transform the two dispersion relations to each other. All non-trivial discrete symmetries of {ω1 , ω2 } are needed to derive the final solution representation, and indeed they are used, for example, to obtain relations (2.15) and (2.16). — With σL = σR and γ L = γ R = 0, solution formulae (2.21) in their proper x-domain of definition reduce to the solution of the whole line problem as given in [3]. — Classical approaches to the problem presented in this section can be found in the literature, for the case γL = 0 = γR . For instance, for one special pair of initial conditions, a solution is presented in [1]. No explicit solution formulae using classical methods with general initial conditions exist to our knowledge. At best, one is left with having to find the solution of an equation involving inverse Laplace transforms, where the unknowns are embedded within these inverse transforms. — The steady-state solution to (2.1) with initial conditions, which decay sufficiently fast to the boundary values (2.3) at ±∞, is easily obtained by letting t → ∞ in (2.21). This gives limt→∞ uR (x, t) = limt→∞ uL (x, t) = (γ L σL + γ R σR )/(σL + σR ). This is the weighted
...................................................
+
∂D+
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
R
+
8
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
u(x,t)
9
–0.01
0.01 x
0
u(x, t) 0.02
–0.01
0.01 0
0.01
x u(x, t)
5 × 10–7 0.01
x
0.01
–0.01 t
0
Figure 4. Results for solutions (2.18) and (2.20) with uL0 (x) = x 2 e(25) x , uR0 (x) = x 2 e−(30) x and σL = 0.02, σR = 0.06, γ L = γ R = 0, t ∈ [0, 0.02] using the hybrid analytical–numerical method of [5]. (Online version in colour.) 2
2
average of the boundary conditions at infinity with weights given by σL and σR . To our knowledge, this is a new result. It is consistent with the steady-state limit (γ L + γ R )/2 for the whole-line problem with initial conditions that limit to different values γ L and γ R as x → ±∞. This result is easily obtained from the solution of the heat equation defined on the whole line using the Fokas method, but it can also be observed by employing piecewise-constant initial data in the classical Green’s function solution, as described in Theorem 4-1 on p. 171 (and comments thereafter) of [8]. It should be emphasized that the steady-state problem for (2.1)–(2.4) or even for the heat equation defined on the whole
...................................................
t
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
5 × 10–7
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
line with different boundary conditions at +∞ and −∞ is ill posed in the sense that the steady solution cannot satisfy the boundary conditions.
2
and 2 −αR x , uR 0 (x) = x e 2
with αL = 25 and αR = 30 (figure 4). The Fourier transforms of these initial conditions may be computed explicitly. We choose σL = 0.02 and σR = 0.06. The initial conditions are chosen so as to satisfy interface boundary conditions (2.4) at t = 0. The results clearly illustrate the discontinuity in the first derivative of the temperature at the interface x = 0.
3. Two finite domains Next, we consider the problem of heat conduction through two walls of finite width (or of two finite rods) with Robin boundary conditions. We seek two functions: uL (x, t),
x ∈ (−a, 0), t ≥ 0 and uR (x, t),
x ∈ (0, b), t ≥ 0,
satisfying the equations uLt (x, t) = σL2 uLxx (x, t),
−a < x < 0, t > 0
(3.1a)
0 < x < b, t > 0,
(3.1b)
uL (x, 0) = uL0 (x),
−a < x < 0
(3.2a)
uR (x, 0) = uR 0 (x),
0 < x < b,
(3.2b)
and 2 R uR t (x, t) = σR uxx (x, t),
the initial conditions
and
the boundary conditions f L (t) = α1 uL (−a, t) + α2 uLx (−a, t),
t>0
(3.3a)
and f R (t) = α3 uR (b, t) + α4 uR x (b, t),
t>0
(3.3b)
and the continuity conditions uL (0, t) = uR (0, t),
t>0
(3.4a)
and σL2 uLx (0, t) = σR2 uR x (0, t),
t > 0,
(3.4b)
as illustrated in figure 5, where a > 0, b > 0 and αi , 1 ≤ i ≤ 4 constant. If α1 = α3 = 0, then Neumann boundary conditions are prescribed, whereas if α2 = α4 = 0 then Dirichlet conditions are given. As before, we have the local relations (e−ikx+(σL k) t uL (x, t))t = (σL2 e−ikx+(σL k) t (uLx (x, t) + ikuL (x, t)))x
(3.5a)
R (e−ikx+(σR k) t uR (x, t))t = (σR2 e−ikx+(σR k) t (uR x (x, t) + iku (x, t)))x .
(3.5b)
2
2
and 2
2
...................................................
uL0 (x) = x2 eαL x
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
Using a slight variation on the method presented in [5], one can compute solutions (2.21) numerically with specified initial conditions. We plot solutions for the case of vanishing boundary conditions (γ L = γ R = 0) with
10
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
utL = sL2 uLxx
utR = sR2 uRxx
11
x
Figure 5. The heat equation for two finite rods.
We define the time transforms of the initial and boundary data and the spatial transforms of u for k ∈ C as follows: uˆ L0 (k) = uˆ R 0 (k) = fˆL (ω, t) = hL1 (ω, t) = hR 1 (ω, t) = g1 (ω, t) = g0 (ω, t) =
and
0 −a
b 0
t
0
t
0
t
0
t
0
t
0
e−ikx uL0 (x) dx,
e−ikx uR 0 (x) dx, eωs fL (s) ds,
uˆ L (k, t) = uˆ R (k, t) =
fˆR (ω, t) =
eωs uLx (−a, s) ds, eωs uR x (b, s) ds, eωs uLx (0, s) ds = eωs uL (0, s) ds =
t 0
0 −a
b
σR2 σL2 t
t 0
e−ikx uR (x, t) dx,
0
eωs fR (s) ds,
hL0 (ω, t) = hR 0 (ω, t) =
e−ikx uL (x, t) dx,
t
t
eωs uL (−a, s) ds,
0
eωs uR (b, s) ds,
0
eωs uR x (0, s) ds
eωs uR (0, s) ds.
0
Using Green’s theorem on the domains [−a, 0] × [0, t] and [0, b] × [0, t], respectively, we have the global relations 2
e(σL k) t uˆ L (k, t) = σL2 (g1 ((σL k)2 , t) + ikg0 ((σL k)2 , t)) − eika σL2 (hL1 ((σL k)2 , t) + ikhL0 ((σL k)2 , t)) + uˆ L0 (k)
(3.6a)
and 2 R 2 2 2 e(σR k) t uˆ R (k, t) = e−ikb σR2 (hR 1 ((σR k) , t) + ikh0 ((σR k) ), t) − σL g1 ((σR k) , t) 2
− ikσR2 g0 ((σR k)2 , t) + uˆ R 0 (k).
(3.6b)
Both equations are valid for all k ∈ C, in contrast to (2.8) and (2.10). Using the invariance of ωL (k) = (σL k)2 and ωR (k) = (σR k)2 under k → −k, we obtain 2
e(σL k) t uˆ L (−k, t) = σL2 (g1 ((σL k)2 , t) − ikg0 ((σL k)2 , t)) − e−ika σL2 (hL1 ((σL k)2 , t) − ikhL0 ((σL k)2 , t)) + uˆ L0 (−k)
(3.7a)
and 2
2 R 2 e(σR k) t uˆ R (−k, t) = eikb σR2 (hR 1 ((σR k) , t) − ikh0 ((σR k) ), t)
− σL2 g1 ((σR k)2 , t) + ikσR2 g0 ((σR k)2 , t) + uˆ R 0 (−k).
(3.7b)
...................................................
b
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
0
–a
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
Im(k)
12 ...................................................
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
D0+
Re(k)
D0–
Figure 6. Deformation of the contours in figure 2 away from k = 0. (Online version in colour.)
Inverting the Fourier transform in (3.6a), 1 ∞ ikx−(σL k)2 t 2 uL (x, t) = e σL (g1 ((σL k)2 , t) + ikg0 ((σL k)2 , t)) dk 2π −∞ 1 ∞ ik(a+x)−(σL k)2 t 2 L − e σL (h1 ((σL k)2 , t) + ikhL0 ((σL k)2 , t)) dk 2π −∞ 1 ∞ ikx−(σL k)2 t L uˆ 0 (k) dk. + e 2π −∞
(3.8)
The integrand of the first integral is entire and decays as k → ∞ for k ∈ C− \ D− . The second integral has an integrand that is entire and decays as k → ∞ for k ∈ C+ \ D+ . It is convenient to deform both contours away from k = 0 to avoid singularities in the integrands that become apparent in what follows. Initially, these singularities are removable, as the integrands are entire. Writing integrals of sums as sums of integrals, the singularities may cease to be removable. With the deformations away from k = 0, the apparent singularities are no cause for concern. In other − − words, we deform D+ to D+ 0 and D to D0 as shown in figure 6. Thus, −1 2 uL (x, t) = eikx−(σL k) t σL2 (g1 ((σL k)2 , t) + ikg0 ((σL k)2 , t)) dk 2π ∂D−0 1 2 − eik(a+x)−(σL k) t σL2 (hL1 ((σL k)2 , t) + ikhL0 ((σL k)2 , t)) dk + 2π ∂D0 1 ∞ ikx−(σL k)2 t L uˆ 0 (k) dk. + e (3.9) 2π −∞ To obtain the solution on the right, we apply the inverse Fourier transform to (3.6b) 1 ∞ ik(x−b)−(σR k)2 t 2 R 2 e σR (h1 ((σR k)2 , t) + ikhR uR (x, t) = 0 ((σR k) , t)) dk 2π −∞ 1 ∞ ikx−(σR k)2 t − e (ikσR2 g0 ((σR k)2 , t) + σL2 g1 ((σR k)2 , t)) dk 2π −∞ 1 ∞ ikx−(σR k)2 t R uˆ 0 (k) dk. e + 2π −∞
(3.10)
The integrand of the first integral is entire and decays as k → ∞ for k ∈ C− \ D− . The second integral has an integrand that is entire and decays as k → ∞ for k ∈ C+ \ D+ . We deform the
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
contours as above to obtain 2
∂D− 0
−
1 2π
+
1 2π
2 R 2 eik(x−b)−(σR k) t σR2 (hR 1 ((σR k) , t) + ikh0 ((σR k) , t)) dk
2
∂D+ 0
∞
−∞
eikx−(σR k) t (ikσR2 g0 ((σR k)2 , t) + σL2 g1 ((σR k)2 , t)) dk 2
eikx−(σR k) t uˆ R 0 (k) dk.
(3.11)
Taking the time transform of the boundary conditions results in fˆ L (ω, s) = α1 hL0 (ω, t) + R α2 hL1 (ω, t) and fˆ R (ω, s) = α3 hR 0 (ω, t) + α4 h1 (ω, t). These two equations together with (3.7a,b) are four equations to be solved for the four L R unknowns hL0 (ω, t), hR 0 (ω, t), h1 (ω, t) and h1 (ω, t). The resulting expressions are substituted in (3.9) and (3.11). Although we could solve this problem in its full generality, we restrict to the case of Dirichlet boundary conditions (α2 = α4 = 0), to simplify the already cumbersome formulae below. Then, hL0 (ω, t) and hR 0 (ω, t) are determined and we solve two equations for two unknowns. System (3.7a,b) is not solvable for hL1 (ω, t) and hR 1 (ω, t) if L (k) = 0, where L (k) = π (σL (e2iak + 1)(e2ibkσL /σR − 1) + σR (e2iak − 1)(e2ibkσL /σR + 1))
bkσL 2iak 2ibkσL /σR + σR tan(ak) . + 1)(e + 1) σL tan = iπ (e σR It is easily seen that all values of k satisfying this (including k = 0) are on the real line. Thus on the contours, the equations are solved without problem, resulting in the expressions below. As before, the right-hand sides of these expressions involve uˆ L (k, t) and uˆ L (k, t), evaluated at a variety of arguments. All terms with such dependence are written out explicitly below. Terms that depend on only known quantities are contained in KL and KR , the expressions for which are given later. uL (x, t) = KL + + + + + + + +
∂D− 0
eikx (σL + σR ) + eik(x+2bσL /σR ) (σR − σL ) L uˆ (k, t) dk 2 L (k)
∂D− 0
−eik(x+2a) (σL + σR ) + eik(x+2a+2bσL /σR ) (σL − σR ) L uˆ (−k, t) dk 2 L (k)
σL eik(x+2a+2bσL /σR ) R kσL uˆ , t dk L (k) σR
−σL eik(x+2a) R −kσL uˆ , t dk L (k) σR
∂D+ 0
eik(x+2a) (σR − σL ) + eik(x+2a+2bσL /σR ) (σL + σR ) L uˆ (k, t) dk 2 L (k)
∂D− 0
∂D− 0
∂D+ 0
∂D+ 0
∂D+ 0
−eik(x+2a) (σL + σR ) + eik(x+2a+2bσL /σR ) (σL − σR ) L uˆ (−k, t) dk 2 L (k)
σL eik(x+2a+2bσL /σR ) R kσL ˆu , t dk L (k) σR
−σL eik(x+2a) R −kσL , t dk uˆ L (k) σR
(3.12)
...................................................
−1 2π
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
uR (x, t) =
13
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
and
14
+ + + + + +
kσR ,t σL
dk
−kσR , t dk σL
∂D+ 0
−σR eik(x+2aσR /σL ) L uˆ R (k)
∂D+ 0
eik(2b+x) (σL − σR ) + eik(x+2b+2aσR /σL ) (σL + σR ) R uˆ (k, t) dk 2 R (k)
∂D− 0
eikx (σR − σL ) − eik(x+2aσR /σL ) (σL + σR ) R uˆ (−k, t) dk 2 R (k)
σR eikx L kσR −σR eik(x+2aσR /σL ) L −kσR uˆ uˆ , t dk + , t dk R (k) σL R (k) σL ∂D− 0
∂D− 0
eikx (σL + σR ) − eik(x+2aσR /σL ) (σL − σR ) R uˆ (k, t) dk 2 R (k)
∂D− 0
eikx (σR − σL ) − eik(x+2aσR /σL ) (σL + σR ) R uˆ (−k, t) dk, 2 R (k)
∂D+ 0
(3.13)
where R (k) = L (kσR /σL ). The integrands written explicitly in (3.12) and (3.13) decay in the regions around whose boundaries they are integrated. Thus, using Jordan’s lemma and Cauchy’s theorem, these integrals are shown to vanish. Thus, the final solution is given by KL and KR . Proposition 3.1. The solution of heat transfer problem (3.1)–(3.4) is given by uL (x, t) = KL =
∞ −∞
+ + + + +
+
2
∂D− 0
∂D− 0
∂D− 0
−eikx−(σL k) t (σL + σR ) + eik(x+2bσL /σR )−(σL k) t (σL − σR ) L uˆ 0 (k) dk 2 L (k)
2
2
2
2
∂D− 0
eik(x+2a)−(σL k) t (σL + σR ) + eik(x+2a+2bσL /σR )−(σL k) t (σR − σL ) L uˆ 0 (−k) dk 2 L (k)
2 −σL eik(x+2a+2bσL /σR )−(σL k) t R kσL uˆ 0 dk L (k) σR
2 σL eik(x+2a)−(σL k) t R −kσL uˆ 0 dk L (k) σR
∂D+ 0
ikσL2 eik(x+a)−(σL k) t (σL + σR ) − ikσL2 eik(x+a+2bσL /σR )−(σL k) t (σL − σR ) L fˆ ((σL k)2 , t) dk α1 L (k)
∂D+ 0
−2ikσL2 σR eik(x+2a+bσL /σR )−(σL k) t (1 + σL σR ) R fˆ ((σL k)2 , t) dk α3 L (k)
∂D+ 0
−eik(x+2a)−(σL k) t (σR − σL ) − eik(x+2a+2bσL /σR )−(σL k) t (σL + σR ) L uˆ 0 (k) dk 2 L (k)
2
∂D− 0
∂D− 0
2
2
+
−2ikσL2 σR eik(x+2a+bσL /σR )−(σL k) t R fˆ ((σL k)2 , t) dk α3 L (k)
ikσL2 eik(x+a)−(σL k) t (σL + σR ) − ikσL2 eik(x+a+2bσL /σR )−(σL k) t (σL − σR ) L fˆ ((σL k)2 , t) dk α1 L (k)
+
2
eikx−(σL k) t uˆ L0 (k) dk +
2
2
2
2
...................................................
∂D+ 0
eikx σR L uˆ R (k)
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
uR (x, t) = KR +
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
+
∂D+ 0
∂D+ 0
∂D+ 0
2
15
(3.14)
for −a < x < 0, and for 0 < x < b uR (x, t) = KR =
∞ −∞
+ + + +
+ + +
2
∂D− 0
∂D− 0
∂D− 0
−eik(x+2b)−(σR k) t (σL − σR ) − eik(x+2b+2aσR /σL )−(σR k) t (σL + σR ) R uˆ 0 (k) dk 2 R (k)
∂D− 0
eikx−(σR k) t (σL − σR ) + eik(x+2aσR /σL )−(σR k) t (σL + σR ) R uˆ 0 (−k) dk 2 R (k)
∂D+ 0
2ikσL σR2 eik(x+aσR /σL )−(σR k) t L fˆ ((σR k)2 , t) dk α1 R (k)
2
∂D− 0
2
2
2
2
2
2
+
2ikσL σR2 eik(x+aσR /σL )−(σR k) t L fˆ ((σR k)2 , t) dk α1 R (k)
−ikσR2 eik(x+b)−(σR k) t (σL − σR ) − ikσR2 eik(x+b+2aσR /σL )−(σR k) t R fˆ ((σR k)2 , t) dk α3 R (k)
2 2 −σR eikx−(σR k) t L kσR σR eik(x+2aσR /σL )−(σR k) t L −kσR uˆ 0 uˆ 0 dk + dk R (k) σL R (k) σL ∂D− 0
+
2
eikx−(σR k) t uˆ R 0 (k) dk +
∂D+ 0
−ikσR eik(x+b)−(σR k) t (σL − σR ) − ikσR2 eik(x+b+2aσR /σL )−(σR k) t (σL + σR ) R fˆ ((σR k)2 , t) dk α3 R (k)
2 2 −σR eikx−(σR k) t L kσR σR eik(x+2aσR /σL )−(σR k) t L −kσR uˆ 0 uˆ 0 dk + dk R (k) σL R (k) σL ∂D+ 0
∂D+ 0
−eikx−(σR k) t (σL + σR ) + eik(x+2aσR /σL )−(σR k) t (σR − σL ) R uˆ 0 (k) dk 2 R (k)
∂D+ 0
eikx−(σR k) t (σL − σR ) + eik(x+2aσR /σL )−(σR k) t (σL + σR ) R uˆ 0 (−k) dk. 2 R (k)
2
∂D+ 0
2
2
2
2
2
(3.15)
Remarks. — The solution of the problem posed in (3.1)–(3.4) may be obtained using the classical method of separation of variables and superposition [1]. The solutions uL (x, t) and uR (x, t) are given by a series of eigenfunctions with eigenvalues that satisfy a transcendental equation, closely related to the equation L (k) = 0. This series solution may be obtained + from proposition 3.1 by deforming the contours along ∂D− 0 and ∂D0 to the real line, including small semicircles around each root of either L (k) or R (k), depending on whether uL (x, t) or uR (x, t) is being calculated. Indeed, this is allowed because all integrands decay in the wedges between these contours and the real line, and the zeros of L (k) and R (k) occur only on the real line, as stated above. Careful calculation of all different contributions, following the examples in [2,3], shows that the contributions along the real line cancel, leaving only residue contributions from the small circles. Each residue contribution corresponds to a term in the classical series solution. It is not necessarily beneficial to leave the form of the solution in proposition 3.1 for the series
...................................................
+
2
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
+
eik(x+2a)−(σL k) t (σL + σR ) + eik(x+2a+2bσL /σR )−(σL k) t (σR − σL ) L uˆ 0 (−k) dk 2 L (k)
2 −σL eik(x+2a+2bσL /σR )−(σL k) t R kσL uˆ 0 dk L (k) σR
2 σL eik(x+2a)−(σL k) t R −kσL uˆ 0 dk, L (k) σR
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
σR2 (γ R − γL ) bσL2 + aσR2
x+
bγ L σL2 + aγ R σR2 bσL2 + aσR2
−a < x < 0
,
and uR (x, t) =
σL2 (γ R − γL ) bσL2
+ aσR2
x+
bγ L σL2 + aγ R σR2 bσL2 + aσR2
0 < x < b,
,
which is piecewise linear and continuous. — A more direct way to recover only the steady-state solution to (3.1) with limt→∞ f L (t) = f¯L and limt→∞ f R (t) = f¯R constant is to write the solution as the superposition of two parts: uL (x, t) = u¯ L (x) + uˇ L (x, t) and uR (x, t) = u¯ R (x) + uˇ R (x, t). The first parts u¯ L and u¯ R satisfy the boundary conditions as t → ∞ and the stationary heat equation. In other words, 0 = σL2 u¯ Lxx (x),
−a < x < 0,
0 = σR2 u¯ R xx (x),
0 < x < b,
f¯L = α1 u¯ L (−a) + α2 u¯ Lx (−a) f¯R = α3 u¯ R (b) + α4 u¯ R x (b).
and
A piecewise linear ansatz with the imposition of the interface conditions results in linear algebra for the unknown coefficients [1]. With the steady-state solution in hand, the second (time-dependent) parts uˇ L and uˇ R satisfy the initial conditions modified by the steady-state solution and the boundary conditions minus their value as t → ∞ uˇ Lt (x, t) = σL2 uˇ Lxx (x, t),
a < x < 0, t > 0,
2 ˇR uˇ R xx (x, t), t (x, t) = σR u
0 < x < b, t > 0,
ˇL
L
¯L
L
u (x, 0) = u (x, 0) − u (x),
−a < x < 0,
uˇ R (x, 0) = uR (x, 0) − u¯ R (x),
0 < x < b,
¯L
ˇL
f (t) − f = α1 u and
(−a, t) + α2 uˇ Lx (−a, t),
f R (t) − f¯R = α3 uˇ R (b, t) + α4 uˇ R x (b, t),
t>0
t > 0,
where, as usual, we impose continuity of temperature and heat flux at the interface x = 0. The dynamics of the solution is described by uˇ L and uˇ R , both of which decay to zero as t → ∞. Their explicit form is easily found using the method described in this section.
...................................................
uL (x, t) =
16
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
representation, as the latter depends on the roots of L (k) and R (k), which are not known explicitly. By contrast, the representation of proposition 3.1 depends on known quantities only and may be readily computed, using one’s favourite parameterization of + the contours ∂D− 0 and ∂D0 . — Similarly, the familiar piecewise linear steady-state solution of (3.1)–(3.4) with Dirichlet boundary conditions [1] can be observed from (3.14) and (3.15) by choosing initial conditions that decay appropriately and constant boundary conditions fL (t) = γ L and fR (t) = γ R . It is convenient to choose zero initial conditions, as the initial conditions do not affect the steady state. As above, the contours are deformed so that they are along the real line with semicircular paths around the zeros of L (k) and R (k), including − k = 0. As one of these deformations arises from D+ 0 while the other comes from D0 , the contributions along the real line cancel each other, while the semicircles add to give full residue contributions from the poles associated with the zeros of L (k). All such residues vanish as t → ∞, except at k = 0. It follows that the steady-state behaviour is determined by the residue at the origin. This results in
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
4. Other problems
17
In this section, we consider the heat equation defined on two semi-infinite rods enclosing a single rod of length 2a. We seek three functions: uL (x, t), x ∈ (−∞, −a),
uM (x, t), x ∈ (−a, b)
and uR (x, t), x ∈ (b, ∞),
with t ≥ 0, satisfying the equations
and
uLt (x, t) = σL2 uLxx (x, t),
−∞ < x < −a, t > 0,
(4.1a)
2 M uM t (x, t) = σM u xx (x, t),
−a < x < a, t > 0
(4.1b)
2 R uR t (x, t) = σR uxx (x, t),
a < x < ∞, t > 0,
(4.1c)
uLt (x, 0) = uL0 (x),
−∞ < x < −a,
(4.2a)
uM t (x, 0) = uM 0 (x),
−a < x < a
(4.2b)
with the initial conditions being
R uR t (x, 0) = u0 (x),
and
a < x < ∞,
(4.2c)
the asymptotic conditions being
and
lim uL (x, t) = 0, x→−∞ t
t≥0
(4.3a)
lim uR (x, t) = 0, x→∞ t
t≥0
(4.3b)
and the continuity conditions being uL (−a, t) = uM (−a, t), uM (a, t) = uR (a, t),
t ≥ 0, t ≥ 0,
2 M σL2 uLx (−a, t) = σM ux (−a, t),
and
2 M σM ux (a, t) = σR2 uR x (a, t),
(4.4a) (4.4b)
t≥0
t ≥ 0,
(4.4c) (4.4d)
as shown in figure 7. Asymptotic conditions (4.3) are not the most general possible but are used here to simplify calculations. To our knowledge, no aspect of this problem is addressed in the literature. To solve this problem, one would proceed with the following steps: (i) Write the local relations in each of the three domains: left, middle and right. (ii) Use Green’s theorem to define the three global relations, keeping track of where they are valid in the complex k plane. (iii) Solve the three global relations for uL (x, t), uM (x, t) and uR (x, t) by inverting the Fourier transforms. (iv) Using the k → −k symmetry of the dispersion relationships on the three original global relations, write down three more global relations. This uses the discrete symmetry of each individual dispersion relation. − (v) Deform the integrals of the solution expressions to D+ 0 or D0 as dictated by the region of analyticity of the integrand.
...................................................
(a) Infinite domain with three layers
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
With the basic principles of the method outlined in the previous two sections, problems with more layers, both finite and infinite, may be addressed. We state two additional problems below and include solutions for specific initial conditions. More complete solutions can be found in the electronic supplementary material.
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
utL = sL2 uLxx
utM = sM2 uMxx
utR = sR2 uRxx
18
x
Figure 7. The heat conduction problem for a single rod of length 2a between two semi-infinite rods.
(vi) Use the discrete symmetries of the dispersion relationship to the collection of global relations for g0 (ω, t), g1 (ω, t), h0 (ω, t) and h1 (ω, t). Care should be taken that relations can only be solved simultaneously provided their regions of validity overlap. At this stage, as in the previous two sections, the discrete symmetries from one dispersion relation to another come into play. (vii) Terms containing unknown functions are shown to be zero by examining the regions of analyticity and decay for the relevant integrands and the use of Jordan’s lemma. For brevity of the presentation, we will assume that uM (x, 0) = 0 = uR (x, 0). The solution to the non-restricted problem can be found in the electronic supplementary material. After defining the transforms uˆ L0 (k) = uˆ M (k) =
−a
e−ikx uL0 (x) dx,
−∞ a
e−ikx uM (x, t) dx,
−a
t
uˆ L (k, t) =
uˆ R (k, t) = 2 σM
t
−∞ ∞
e−ikx uL (x, t) dx, e−ikx uR (x, t) dx,
a
eωs uM x (−a, s) ds, σL2 0 t t h0 (ω, t) = eωs uL (−a, s) ds = eωs uM (−a, s) ds,
h1 (ω, t) =
0
eωs uLx (−a, s) ds =
−a
0
0
0
eωs uM x (a, s) ds =
σR2
t
eωs uR x (a, s) ds 2 σM 0 t t g0 (ω, t) = eωs uM (a, s) ds = eωs uR (a, s) ds,
g1 (ω, t) = and
t
0
0
we proceed as outlined above. The solution formulae are given in the following proposition: Proposition 4.1. The solution of the heat transfer problem (4.1)–(4.4) with uM (x, 0) = uR (x, 0) = 0 is uL (x, t) =
1 2π −
∞
1 2
−∞
2
eikx−(σL k) t uˆ L0 (k) dk
2
eik(x+2a)−(σL k) t ((σL + σM )(σM − σR ) L (k)
∂D− 0
+ e4ikaσL /σM (σL − σM )(σM + σR ))uˆ L0 (−k) dk, for −∞ < x < −a with L (k) = π ((σL − σM )(σM − σR ) + e4iakσL /σM (σL + σM )(σM + σR )), uM (x, t) = −σM (σM − σR )
2
∂D− 0
+ σM (σM + σR )
−kσM dk σL
2 eik(x+a−aσM /σL )−(σM k) t L kσM uˆ 0 dk, R (kσM /σR ) σL
eik(x+a+aσM /σL )−(σM k) t L uˆ 0 L (kσM /σL )
∂D+ 0
...................................................
a
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
–a
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
utL = sL2 uLxx
utM = sM2 uM xx
utR = sR2 uRxx
19
x c
Figure 8. The heat equation for three finite layers.
for −a < x < b with R (k) = π ((σL + σM )(σM + σR ) + e4iakσR /σM (σL − σM )(σM − σR )). Lastly, the expression for uR (x, t), valid for x > a, is uR (x, t) = 2σM σR
eik(x−a−aσR /σL +2aσR /σM )−(σR k) t L uˆ 0 R (k) 2
∂D+ 0
kσR σL
dk.
(b) Finite domain with three layers We consider the heat conduction problem in three rods of finite length as shown in figure 8. We seek three functions: uL (x, t),
x ∈ (−a, 0), t ≥ 0,
uM (x, t),
x ∈ (0, b), t ≥ 0 and uR (x, t),
x ∈ (b, c), t ≥ 0,
satisfying the equations uLt (x, t) = σL2 uLxx (x, t), M
u and
−a < x < 0, t > 0,
(4.5a)
0 < x < b, t > 0
(4.5b)
2 M t (x, t) = σM u xx (x, t),
2 R uR t (x, t) = σR uxx (x, t),
b < x < c, t > 0,
(4.5c)
uLt (x, 0) = uL0 (x),
−∞ < x < −a,
(4.6a)
uM t (x, 0) = uM 0 (x),
−a < x < a
(4.6b)
with the initial conditions being
R uR t (x, 0) = u0 (x),
and
a < x < ∞,
(4.6c)
the boundary conditions being f L (t) = α1 uL (−a, t) + α2 uLx (−a, t), and
f R (t) = α3 uR (c, t) + α4 uR x (c, t),
t>0
t>0
(4.7a) (4.7b)
and the continuity conditions being uL (0, t) = uM (0, t), M
R
u (b, t) = u (b, t),
t ≥ 0,
(4.8a)
t ≥ 0,
(4.8b)
2 M σL2 uLx (0, t) = σM ux (0, t),
and
2 M σM ux (b, t) = σR2 uR x (b, t),
t≥0
(4.8c)
t ≥ 0.
(4.8d)
The solution process is as before, following the steps outlined in the previous section. For simplicity, we assume Neumann boundary data (α1 = α3 = 0), zero boundary conditions (f L (t) = f R (t) = 0) and uM (x, 0) = 0 = uR (x, 0). The solution with uM (x, 0) = 0 and uR (x, 0) = 0 is given in the
...................................................
b
0
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
–a
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
electronic supplementary material. We define
gL1 (ω, t) =
b
e−ikx uL0 (x) dx,
uˆ L (k, t) =
e−ikx uM (x, t) dx,
0
t
eωs uLx (−a, s) ds,
2 σM
t
t
e−ikx uR (x, t) dx,
b
eωs uL (−a, s) ds,
0
eωs uM x (0, s) ds, σL2 0 t t hL0 (ω, t) = eωs uL (0, s) ds = eωs uM (0, s) ds, hL1 (ω, t) =
0
eωs uLx (0, s) ds =
gL0 (ω, t) =
e−ikx uL (x, t) dx,
c
uˆ R (k, t) =
0
t
−a
0
t
0
0
eωs uM x (b, s) ds =
σR2
t
eωs uR x (b, s) ds, 2 σM 0 t t gR eωs uM (b, s) ds = eωs uR (b, s) ds 0 (ω, t) = gR 1 (ω, t) =
0
hR 1 (ω, t) =
and
t
0
0
eωs uR x (c, s) ds,
hR 0 (ω, t) =
t
eωs uR (c, s) ds.
0
The solution is given by the following proposition: Proposition 4.2. The solution to the heat transfer problem (4.5)–(4.8) with α1 = α3 = 0, f L (t) = f R (t) = 0 and uM (x, 0) = uR (x, 0) = 0 is 1 ∞ ikx−(σL k)2 t L uˆ 0 (k) dk uL (x, t) = e 2π −∞ 2 eikx−(σL k) t 2ibkσL /σR e (σL − σM )(σM − σR ) + 2 L (k) ∂D− 0 + e2ikσL (c/σR +b/σM ) (σL + σM )(σM − σR ) + e2ickσL /σR (σL − σM )(σM + σR ) +e2ibkσL (1/σR +1/σM ) (σL + σM )(σM + σR ) uˆ L0 (k) dk +
2
∂D− 0
eik(x+2a)−(σL k) 2 L (k)
t
(e2ibkσL /σR (σL − σM )(σM − σR )
+ e2ikσL (c/σR +b/σM ) (σL + σM )(σM − σR ) + e2ickσL /σR (σL − σM )(σM + σR ) +e2ibkσL (1/σR +1/σM ) (σL + σM )(σM + σR ) uˆ L0 (−k) dk +
∂D+ 0
2 eik(x+2a)−(σL k) t 2ikσL (c/σR +b/σM ) e (σL − σM )(σM − σR ) 2 L (k)
+ e2ibkσL /σR (σL + σM )(σM − σR ) + e2ibkσL (1/σR +1/σM ) (σL − σM )(σM + σR ) +e2ickσL /σR (σL + σM )(σM + σR ) uˆ L0 (k) dk +
2
∂D+ 0
eik(x+2a)−(σL k) 2 L (k)
t
(e2ibkσL /σR (σL − σM )(σM − σR )
+ e2ikσL (c/σR +b/σM ) (σL + σM )(σM − σR ) + e2ickσL /σR (σL − σM )(σM + σR ) +e2ibkσL (1/σR +1/σM ) (σL + σM )(σM + σR ) )uˆ L0 (−k) dk,
...................................................
uˆ M (k) =
−a
20
0
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
uˆ L0 (k) =
0
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
for −a < x < 0 with
21
+ e2ickσL /σR (σL − σM )(σM + σR ) + e2ik(a+bσL /σR +bσL /σM ) (σM − σL )(σM + σR ) + e2ibk(σL /σR +σL /σM ) (σL + σM )(σM + σR ) − e2ik(cσL /σR +a) (σL + σM )(σM + σR )). Next, uM (x, t) =
−eik(x+b)−(σM k) t σM 2ickσM /σR e (σM − σR ) − M (k) ∂D0 kσ M dk +e2ibkσM /σR (σM + σR ) uˆ L0 σL 2 −eik(x+b+2aσM /σL )−(σM k) t σM 2ickσM /σR e + (σM − σR ) M (k) ∂D− 0 −σ M dk +e2ibkσM /σR (σM + σR ) uˆ L0 σL 2 −eik(x+b)−(σM k) t σM 2ickσM /σR e + (σM − σR ) M (k) ∂D+ 0 kσ M dk +e2ibkσM /σR (σM + σR ) uˆ L0 σL 2 eik(x+b+2aσM /σL )−(σM k) t σM 2ickσM /σR e − (σM − σR ) M (k) ∂D+ 0 −kσ M dk, +e2ibkσM /σR (σM + σR ) uˆ L0 σL 2
for 0 < x < b, with M (k) = π (eik(aσM /σL +b+cσM /σR ) (σL − σM )(σM − σR ) + e2bikσM /σR (σM − σL )(σM − σR ) + e2ikσM (b/σR +a/σL ) (σL + σM )(σM − σR ) + e2ik(cσM /σR +b) (σL + σM )(σR − σM ) + e2ik(aσM /σL +bσM /σR +b) (σL − σM )(σM + σR ) + e2ickσM /σR (σM − σL )(σM + σR ) − e2ibk(σM /σR +1) (σL + σM )(σM + σR ) + e2ikσM (c/σR +a/σL ) (σL + σM )(σM + σR )) and uR (x, t) =
−2eik(x+b+bσR /σM )−(σR k) t (σM σR )uˆ L0 R (k) 2
∂D− 0
+ + +
∂D− 0
∂D+ 0
∂D+ 0
kσR σL
dk
2 −kσR −2 eik(x+b+bσR /σM +2aσR /σL )−(σR k) t (σM σR )uˆ L0 dk R (k) σL
2 kσR 2 eik(x+b+bσR /σM )−(σR k) t (σM σR )uˆ L0 dk R (k) σL
2 −kσR 2 eik(x+b+bσR /σM +2aσR /σL )−(σR k) t (σM σR )uˆ L0 dk, R (k) σL
...................................................
+ e2ikσL (c/σR +b/σM ) (σL + σM )(σM − σR ) + e2ik(bσL /σR +a) (σL + σM )(σR − σM )
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
L (k) = π (e2ibkσL /σR (σL − σM )(σM − σR ) + e2ik(cσL /σR +bσL /σM +a) (σM − σL )(σM − σR )
Downloaded from rspa.royalsocietypublishing.org on September 29, 2014
for b < x < c with
22
R (k) = π (e2ibk (σL − σM )(σM − σR ) + e2ik(aσR /σL +bσR /σM +c) (σM − σL )(σM − σR )
+ e2ibk(1+σR /σM ) (σL + σM )(σM + σR ) − e2ik(aσR /σL +c) (σL + σM )(σM + σR )). Acknowledgements. Any opinions, findings and conclusions or recommendations expressed in this paper are those of the authors and do not necessarily reflect the views of the funding sources. We thank the American Institute for Mathematics (AIM, Palo Alto, CA, USA) for its hospitality, while some of this work was being conducted. Funding statement. This work was generously supported by the National Science Foundation under grant no. NSF-DMS-1008001 (B.D., N.E.S.). N.E.S. also acknowledges support from the National Science Foundation under grant no. NSF-DGE-0718124.
References 1. Hahn D, Özisik M. 2012 Heat conduction, 3rd edn. Hoboken, NJ: John Wiley & Sons, Inc. 2. Deconinck B, Trogdon T, Vasan V. 2014 The method of Fokas for solving linear partial differential equations. SIAM Rev. 56, 159–186. (doi:10.1137/110821871) 3. Fokas AS. 2008 A unified approach to boundary value problems. CBMS-NSF Regional Conference Series in Applied Mathematics, vol. 78. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM). 4. Fokas AS, Pelloni B. 2005 A transform method for linear evolution PDEs on a finite interval. IMA J. Appl. Math. 70, 564–587. (doi:10.1093/imamat/hxh047) 5. Flyer N, Fokas AS. 2008 A hybrid analytical-numerical method for solving evolution partial differential equations. I. The half-line. Proc. R. Soc. A 464, 1823–1849. (doi:10.1098/rspa. 2008.0041) 6. Carslaw HS, Jaeger JC. 1959 Conduction of heat in solids, 2nd edn. New York, NY: Oxford University Press. 7. Ablowitz MJ, Fokas AS. 2003 Complex variables: introduction and applications, 2nd edn, Cambridge Texts in Applied Mathematics. Cambridge, UK: Cambridge University Press. 8. Guenther RB, Lee JW. 1996 Partial differential equations of mathematical physics and integral equations. Mineola, NY: Dover Publications Inc. (Corrected reprint of the 1988 original.)
...................................................
+ e2ick (σL − σM )(σM + σR ) + e2ik(a/σL +b+bσR /σM ) (σM − σL )(σM + σR )
rspa.royalsocietypublishing.org Proc. R. Soc. A 470: 20130605
+ e2ik(bσR /σM +c) (σL + σM )(σM − σR ) + e2ik(aσR /σL +b) (σL + σM )(σR − σM )
