Steady Two-Dimensional Conduction: Simple and Double Layer Potentials, Corner Singularities, and Induced Heat Flux

Lienhard, John H.

doi:10.1115/1.4055833

Abstract

Two-dimensional steady-state heat conduction is possible outside closed boundaries on which two isothermal segments at different temperatures are separated by two adiabatic segments. Remarkably, previous research showed that the conduction shape factor for the region exterior to the boundary is equal to that for the interior region, despite the asymmetry and singularity of the boundary heat flux distributions. In this study, classical potential theory is used for the temperature and heat flux distributions as a combination of simple-layer and double-layer potentials, including the relationships between the values inside and outside the boundary curve. Isothermal boundaries exhibit an induced heat flux that varies from point-to-point on the boundary. The induced flux integrates to zero over each isothermal edge. Singularities of the heat flux are identified and resolved. Computations that validate the theory are provided for mixed boundary conditions on a disk and a square. Numerical fits to both the simple-layer and double-layer densities are given for the disk and the square. The analysis explains why the interior and exterior conduction shape factors are equal despite wildly differing heat flux distributions, and the results are compared to a previous study of this configuration. This paper also develops fundamental concepts of potential theory and can serve as a tutorial on the subject.

1 Introduction

In a previous paper [1], I studied the relationship between two-dimensional (2D) heat conduction inside and outside closed curves, showing that the conduction shape factor had the same value for the exterior region and the interior region. The primary method of that study was conformal mapping, showing that the total heat flow between the isothermal sections of the boundary was invariant under mapping. In this study, I address the local heat flux distributions on the inside and outside of the boundary, which can differ widely even though the total heat transfer on either side is the same. I develop a boundary integral formulation of the interior/exterior conduction problem, using methods from classical potential theory, and I apply that formulation to identify the relationship between interior and exterior heat flux, as well as their connection to the shape factors.

Consider a region R bounded by a simple closed curve,

C

⁠, that contains two isothermal segments,

Γ^{h}

and

Γ^{c}

⁠, separated by two adiabatic segments,

Γ_{1}^{a}

and

Γ_{2}^{a}

⁠. Let s be a point on

C

⁠. Steady heat conduction in R is governed by the Laplace equation for temperature, θ, scaled to

- 1 \leq θ \leq 1

\nabla^{2} θ = 0

(1)

θ = {\begin{matrix} + 1 & s \in Γ^{h} \\ - 1 & s \in Γ^{c} \end{matrix}

(2)

\frac{\partial θ}{\partial n} = 0 s \in Γ_{1}^{a} \cup Γ_{2}^{a}

(3)

where n is the outward normal vector of $C$ (Fig. 1(a)).

Fig. 1

View large Download slide

Coordinates for the interior and exterior conduction problems: (a) interior problem and (b) exterior problem

Steady conduction in the region R^e exterior to $C$ (Fig. 1(b)) is governed by the same equations, upon replacing n by the exterior’s outward normal vector, $n_{e}$ ⁠, with the addition of a boundary condition on the far field. In two dimensions, bounded steady-state solutions exist only if no net heat transfer is transferred to the region far from the object, so that heat flows only from the hot segment, $Γ^{h}$ ⁠, to the cold segment of $Γ^{c}$ ⁠. In this case, the far field temperature $θ^{\infty}$ has a specific value between –1 and 1, as discussed in Appendix A and the final section of Ref. [1].

Classically, the Laplace equation in 2D can be solved by a wide variety of techniques [2], including conformal mapping [1,3], potential theory [4,5], Green’s functions [6], and separation of variables. In addition, modern numerical methods, such as the finite element method (FEM), allow quick solutions to the Laplace equation even in relatively complicated domains. The goal of this paper is to develop theoretical understanding, and analytical boundary integral methods are the focus, with FEM used to examine model cases. In particular, we use integral equations based on classical potential theory, rather than the Green’s function integrals attempted in Ref. [1]. As Martinsson notes [7], the Green’s function is known explicitly only on trivial domains, whereas potential theory leads to integral equations with known kernels, for which excellent numerical techniques now exist [8,9].

The conduction shape factor, S, in 2D provides the heat transfer rate,

\dot{Q}

(W/m), between the hot and cold boundaries in terms of the dimensional temperature difference,

Δ T = T_{h} - T_{c}

⁠, and thermal conductivity, λ, as

\dot{Q} = S λ Δ T

⁠. For either the interior or the exterior of

C

⁠, the shape factor may be written as an integral on

Γ^{h}

of the appropriate outward normal derivative of θ [1]

S^{i} = \int_{Γ^{h}} \frac{\partial θ}{\partial n} d l, interior to C

(4a)

S^{e} = \int_{Γ^{h}} \frac{\partial θ^{e}}{\partial n^{e}} d l, exterior to C

(4b)

This paper develops a single equation that describes both the interior and the exterior heat flux on $C$ in terms of potential densities. With that equation, the interior and exterior shape factors are both shown to equal an integral of the simple-layer potential density.

2 Preliminaries: Sources, Dipoles, and Gradients

This section briefly introduces several relationships that will be used in other parts of the paper.

2.1 Fundamental Solution: $E (x | t)$ ⁠.

The fundamental solution of the Laplace equation gives the potential at a point x caused by a unit point source at a point t [2]

- \nabla_{x}^{2} E (x | t) = δ (x - t)

(5)

E (x | t) = \frac{1}{2 π} log \frac{1}{| x - t |} in 2 D

(6)

Obviously, $E (x | t) = E (t | x)$ ⁠.

We may find the integral of the normal derivative of the fundamental solution over the smooth boundary

C

of a region R using the divergence theorem

\int_{R} \nabla_{x}^{2} E (x | t) d x = - \int_{R} δ (x - t) d x

(7)

\int_{C} \frac{\partial E (x | t)}{\partial n_{x}} d σ_{x} = {\begin{matrix} - 1 & t \in R \\ - 1 / 2 & t \in C \\ 0 & t \in R^{e} \end{matrix}

(8)

The physical interpretation is straightforward. For example, if we replace the potential E with the dimensional temperature at x caused by a heat source of power P at t, then $- \nabla_{x}^{2} T = (P / λ) δ (x - t)$ and the leaving heat flux is $- λ \partial T / \partial n_{x}$ ⁠. The net heat flow out of R through $C$ is P if $t \in R$ and is zero for t outside $C$ ⁠. When t is on $C$ ⁠, half the power is released outside $C$ ⁠, with the other half released inside σ and exiting $C$ at other points.

In polar coordinates, consider

x = (r, ϕ)

and

t = (r_{0}, ϕ_{0})

⁠. Then, with reference to Fig. 2, the law of cosines gives

| x - t |^{2} = r_{0}^{2} - 2 r r_{0} cos (ϕ - ϕ_{0}) + r^{2}

(9)

Fig. 2

View large Download slide

Polar coordinates for the fundamental solution

Thus, the fundamental solution in polar coordinates is

E (x | t) = \frac{1}{2 π} log \frac{1}{| x - t |} = \frac{1}{4 π} log \frac{1}{[r_{0}^{2} - 2 r r_{0} cos (ϕ - ϕ_{0}) + r^{2}]}

(10)

The radial derivative is

\frac{\partial E}{\partial r} = - \frac{1}{4 π} \frac{2 r - 2 r_{0} cos (ϕ - ϕ_{0})}{[r_{0}^{2} - 2 r r_{0} cos (ϕ - ϕ_{0}) + r^{2}]}

(11)

The fundamental solution also has this series representation (or multipole expansion)

4 π E (x | t) = - 2 log r_{>} + 2 \sum_{m = 1}^{\infty} r_{<}^{m} r_{>}^{- m} \frac{cos m (ϕ - ϕ_{0})}{m}

(12)

where $r_{>} = \max (r, r_{0})$ and $r_{<} = \min (r, r_{0})$ ⁠.

Letting

| x | = r \to \infty

holding t fixed, with Eqs. or

E (x | t) \sim - \frac{1}{2 π} log r + \frac{1}{2 π} \frac{r_{0}}{r} cos (ϕ - ϕ_{0}) + \dots

(13)

\frac{\partial E}{\partial r} \sim - \frac{1}{2 π r} - \frac{1}{2 π} \frac{r_{0}}{r^{2}} cos (ϕ - ϕ_{0}) - \dots

(14)

2.2 Dipoles.

A dipole is a source-sink pair of equal strength brought to infinitesimal separation. If the source and sink are of strength

1 / ε

at a separation of ε, the result of letting

ε \to 0

is a unit-strength dipole. The dipole has a distinct orientation set by the axis on which the source and sink are aligned. We may think of a source at

t + ε l

approaching a sink at t (Fig. 3). Call the potential of this dipole

D_{l} (x | t)

⁠. In 2D

D_{l} (x | t) = \lim_{ε \to 0} \frac{1}{2 π ε} (log \frac{1}{| x - t - ε l |} - log \frac{1}{| x - t |})

(15a)

= \dots = \frac{cos γ}{2 π | x - t |}

(15b)

where γ is the angle between the vectors $x - t$ and l. Note that $D_{l} (x | t) = - D_{l} (t | x)$ ⁠.

Fig. 3

View large Download slide

Dipole coordinates

2.3 Dipoles as Induced Gradients in Two Dimensions.

The gradient of potential at t produced by a source at x is simply

\nabla_{t} E (t | x) = \nabla_{t} E (x | t)

⁠, and this gradient is collinear with

x - t

(pointing from t to x). The derivative of

E (x | t)

at t in the direction l is

\frac{\partial E}{\partial l} = l \cdot \nabla_{t} E = \frac{cos γ}{2 π | x - t |} = D_{l} (x | t)

(16)

Thus, we may consider $D_{l} (x | t)$ as the gradient of potential (temperature gradient) at t in the l-direction caused by a point source at x.

If t is on a surface, l may be selected as the normal vector at t, $n_{t}$ ⁠. If both x and t are on this surface and $l = n_{t}, cos γ \to 0$ as $| x - t | \to 0$ ⁠. Analysis shows that $D_{n_{t}} (x | t)$ is continuous at x = t when x approaches t along the surface (as opposed to approaching from another direction).

Similarly, the outward normal gradient at x produced by a source at t is

D_{n_{x}} (t | x) = \frac{\partial E (t | x)}{\partial n_{x}}

(17)

From Eqs. and , the outward normal gradient induced at a point x on a 2D boundary by a source density a(t) at another point t is

a (t) D_{n_{x}} (t | x) d l_{t} = a (t) \frac{cos (t - x, n_{x})}{2 π | t - x |} d l_{t}

(18)

where $cos (t - x, n_{x})$ denotes the cosine of the angle between vectors $t - x$ and $n_{x}$ ⁠. This expression may be integrated over the surface that contains the source distribution, a(t), holding the vector $n_{x}$ fixed, to find the outward normal gradient at the single point x caused by the entire source distribution.

2.4 Gradient of a Dipole Potential in Two Dimensions.

Later, we need the gradient of a dipole potential in 2D. Suppose that we have a unit dipole of orientation l at t and that we wish to know the potential’s gradient in the direction k at point x (Fig. 4). We find

\frac{\partial D_{l} (x | t)}{\partial k} = k \cdot \nabla_{x} D_{l} = k \cdot \nabla_{x} (\frac{cos γ}{2 π | x - t |})

(19)

where γ is the angle between the vectors

x - t

and l. For convenience, we may use polar coordinates

(r, γ)

with origin at t and

r = | x - t |

⁠. Then

\nabla_{x} D_{l} = \frac{\partial}{\partial r} (\frac{cos γ}{2 π r}) e_{r} + \frac{1}{r} \frac{\partial}{\partial γ} (\frac{cos γ}{2 π r}) e_{γ}

(20)

= - (\frac{cos γ}{2 π r^{2}}) e_{r} - (\frac{sin γ}{2 π r^{2}}) e_{γ}

(21)

Fig. 4

View large Download slide

Dipole gradient coordinates

Suppose that the unit vector k makes an included angle ψ with respect to

x - t

⁠:

k = cos ψ e_{r} - sin ψ e_{γ}

⁠. We find

\frac{\partial D_{l} (x | t)}{\partial k} = \frac{- cos γ cos ψ + sin γ sin ψ}{2 π r^{2}} = - \frac{cos (γ + ψ)}{2 π r^{2}}

(22)

2.5 Some Properties of Harmonic Functions.

For r > 0, the separation of variables in polar coordinates gives

θ (r, ϕ) = A + B log r + \sum_{m = 1}^{\infty} (c_{m}^{-} r^{- m} + c_{m}^{+} r^{m}) cos (m ϕ + α_{m})

(23)

To have a bounded solution as $r \to \infty$ ⁠, we require B = 0 and $c_{m}^{+} = 0$ ⁠. As consequences, $θ \to A = θ^{\infty}$ as $r \to \infty$ ⁠; and for every circle $r = constant$ ⁠, the average of θ around the circumference is $θ^{\infty}$ ⁠.

With these restrictions, the radial gradient is

\frac{\partial θ (r, ϕ)}{\partial r} = - \sum_{m = 1}^{\infty} \frac{m c_{m}^{-}}{r^{m + 1}} cos (m ϕ + α_{m})

(24)

This gradient $\to 0$ as $r \to \infty$ ⁠.

In particular, for any circle centered on r = 0 surrounding $C$ from the outside, the temperature on the circle is continuous, bounded, can be represented by a Fourier series (cf. Sec. 3.3), has the average temperature $θ^{\infty}$ ⁠, and has zero average heat flux through it.

3 Solution for the Interior and Exterior Dirichlet Problems by a Simple-Layer Potential

In his famous essay of 1828, George Green showed that the solution to the Dirichlet problem in the regions inside (R) and outside (R^e) a 3D surface σ may be written as a simple-layer potential [10]¹

\int_{σ} a (ξ) E (x | ξ) d σ_{ξ} = {\begin{matrix} θ (x) & x \in R \\ θ^{e} (x) & x \in R^{e} \end{matrix}

(25)

where

a (ξ)

is the solution of a Fredholm equation of the first kind

\int_{σ} a (ξ) E (s | ξ) d σ_{ξ} = h (s) s on σ

(26)

Here, h(s) is the temperature distribution on σ, which is the same for the interior and exterior problem. No constraint is placed on the normal derivatives to either side of the boundary, which need not be equal.

The derivation of Eq. defines

a (ξ)

as the difference between the normal derivatives inside and outside σ

a (ξ) = \frac{\partial θ}{\partial n_{ξ}} - \frac{\partial θ^{e}}{\partial n_{ξ}} ξ on σ

(27)

The function $a (ξ)$ may be considered to be a heat source distribution on σ, with heat flowing out to each side.

In 2D, we need a minor modification of Eqs. (25) and (26) to account for the constraint on the far-field temperature. The derivation follows.

3.1 Modifications of Equations (25) and (26) for Two Dimensions.

Consider the exterior region (Fig. 1), R^e, and suppose it to be bounded by the contour

C

and a very distant circular contour of radius r₀, on which

ξ = (r_{0}, ϕ_{0})

⁠. We treat

E (x | ξ)

and

θ^{e} (ξ)

as functions of ξ and apply Green’s second identity to these functions

\begin{matrix} \int_{R^{e}} [E \nabla_{ξ}^{2} θ^{e} - θ^{e} \nabla_{ξ}^{2} E] d R_{ξ} = \int_{C} [E \frac{\partial θ^{e}}{\partial n_{ξ}^{e}} - θ^{e} \frac{\partial E}{\partial n_{ξ}^{e}}] d l_{ξ} \\ + \int_{0}^{2 π} [E \frac{\partial θ^{e}}{\partial r_{0}} - θ^{e} \frac{\partial E}{\partial r_{0}}] r_{0} d ϕ_{0} \end{matrix}

(28)

Here the subscript ξ indicates that differentiation or integration is with respect to the ξ coordinate. The left-hand side simplifies upon substituting Eqs. and

\begin{matrix} θ^{e} (x) = \int_{C} [E \frac{\partial θ^{e}}{\partial n_{ξ}^{e}} - θ^{e} \frac{\partial E}{\partial n_{ξ}^{e}}] d l_{ξ} \\ + \int_{0}^{2 π} [E \frac{\partial θ^{e}}{\partial r_{0}} - θ^{e} \frac{\partial E}{\partial r_{0}}] r_{0} d ϕ_{0} \end{matrix}

(29)

Now we let

r_{0} \to \infty

⁠. We know from Eq. that

θ^{e} \sim θ^{\infty} + c_{1} cos (ϕ_{0} + α_{1}) / r_{0} + c_{2} cos (2 ϕ_{0} + α_{2}) / r_{0}^{2} + \dots

and that

\partial θ^{e} / \partial r_{0} \sim - c_{1} cos (ϕ_{0} + α_{1}) / r_{0}^{2}

⁠. Additionally, from Eqs. and ,

E \sim - log r_{0} / 2 π

and

\partial E / \partial r_{0} \sim - 1 / 2 π r_{0}

⁠. Thus

\int_{0}^{2 π} [E \frac{\partial θ^{e}}{\partial r_{0}} - θ^{e} \frac{\partial E}{\partial r_{0}}] r_{0} d ϕ_{0} \sim O (- \frac{log r_{0}}{r_{0}}) + θ^{\infty} \frac{2 π r_{0}}{2 π r_{0}} \to θ^{\infty}

(30)

Combining with Eq. and setting

n_{e} = - n

\int_{C} [- E \frac{\partial θ^{e}}{\partial n_{ξ}} + θ^{e} \frac{\partial E}{\partial n_{ξ}}] d l_{ξ} + θ^{\infty} = {\begin{matrix} 0 & x \in R \\ θ^{e} (x) & x \in R^{e} \end{matrix}

(31)

Upon adding this result to the analogous interior result (which involves only the contour

C

⁠)

\int_{C} E (x | ξ) a (ξ) d l_{ξ} + θ^{\infty} = {\begin{matrix} θ (x) & x \in R \\ θ^{e} (x) & x \in R^{e} \end{matrix}

(32)

The solution by a simple-layer potential is known to be continuous at

C

⁠, so it follows that:

\int_{C} E (s | ξ) a (ξ) d l_{ξ} + θ^{\infty} = h (s) s on σ

(33)

These results show that $θ^{\infty}$ is imposed by the far field in 2D.

In addition, as

r = | x | \to \infty

in Eq. , with Eq.

\int_{C} E (x | ξ) a (ξ) d l_{ξ} \to - \frac{log | x |}{2 π} \int_{C} a (ξ) d l_{ξ} = 0 for x \in R^{e}

(34)

so that $θ^{e} (x) \to θ^{\infty}$ as expected.

3.2 Two-Dimensional Contours With Zero Net Heat Transfer.

As noted, a 2D exterior conduction problem only has a solution when there is no net heat transfer to the far field. Thus

\underset{= 0}{\underset{︸}{\int_{C} \frac{\partial θ}{\partial n_{ξ}} d l_{ξ}}} - \underset{= 0}{\underset{︸}{\int_{C} \frac{\partial θ^{e}}{\partial n_{ξ}} d l_{ξ}}} = \int_{C} a (ξ) d l_{ξ} = 0

(35)

In the case of an entirely isothermal 2D boundary, $h (s) = θ^{\infty}$ ⁠. Then, from Eq. (33), $a (ξ) = 0$ ⁠, so that $θ (x) = θ^{e} (x) = θ^{\infty}$ from Eq. (32), and $\partial θ / \partial n = \partial θ^{e} / \partial n = 0$ ⁠.

3.3 Example: Simple-Layer Potential on a Unit Disk.

If

h (s) = h (ϕ)

is continuous, we can represent it by a harmonic series. To find the potential function, assume that

a (ϕ_{0}) = cos (m ϕ_{0} + α_{m})

⁠. Integration by parts, with Eq. , gives us

\begin{matrix} \int_{0}^{2 π} E (1, ϕ | 1, ϕ_{0}) a (ϕ_{0}) d ϕ_{0} + θ^{\infty} \\ = - \frac{1}{4 π m} \int_{0}^{2 π} \frac{sin (ϕ - ϕ_{0}) sin (m ϕ_{0} + α_{m})}{[1 - cos (ϕ - ϕ_{0})]} d ϕ_{0} + θ^{\infty} \end{matrix}

(36a)

and with

υ = ϕ_{0} - ϕ

= - \frac{1}{4 π m} \int_{0 - ϕ}^{2 π - ϕ} \frac{sin υ sin [m (υ + ϕ) + α_{m}]}{(1 - cos υ)} d υ + θ^{\infty}

(36b)

and letting

z = e^{i υ}, d υ = d z / i z, cos υ = (z + \bar{z}) / 2, sin υ = (z - \bar{z}) / 2 i

⁠, and doing several lines of algebra, we obtain a contour integral

= ℑ m {\frac{e^{i (m ϕ + α_{m})}}{4 π m} \int_{C} \frac{z^{m - 1} (z + 1)}{z - 1} d z} + θ^{\infty}

(36c)

= ℑ m {\frac{e^{i (m ϕ + α_{m})}}{4 π m} (π i) (2)} + θ^{\infty}

(36d)

= θ^{\infty} + \frac{1}{2 m} cos (m ϕ + α_{m})

(36e)

Hence, when

h (ϕ) = cos (m ϕ + α_{m})

⁠, we have

a (ϕ) = 2 m cos (m ϕ + α_{m})

⁠. Additional calculation gives θ and θ^e just as would the separation of variables solution, Eq.

\begin{matrix} θ (r, ϕ) = \int_{0}^{2 π} E (r, ϕ | 1, ϕ_{0}) 2 m cos (m ϕ + α_{m}) d ϕ_{0} + θ^{\infty} = \dots \\ = ℑ m {\frac{e^{i (m ϕ + α_{m})}}{2 π} \int_{C} \frac{z^{m - 1} (z^{2} - 1)}{(z - r) (z - 1 / r)} d z} + θ^{\infty} \\ = θ^{\infty} + {\begin{matrix} r^{m} cos (m ϕ + α_{m}) & for r < 1 \\ r^{- m} cos (m ϕ + α_{m}) & for r > 1 \end{matrix} \end{matrix}

(37)

4 Heat Flux at the Boundaries

To evaluate the heat flux inside and outside $C$ ⁠, we need to find the local normal derivative of temperature on each side of the isothermal boundaries. The derivative need not be the same inside and outside because the isothermal section behaves as a perfect heat conductor: heat can be redistributed within the boundary without a temperature gradient. This strange behavior also leads to the singularity of the boundary integral. As such, particular care is needed in finding the limiting values of the temperature gradient as the boundary is approached.

The boundary derivatives of the simple-layer potential have been studied for at least two centuries, mainly in connection with the electrostatic and gravitational theories. We develop the relevant relationships briefly in this section, for their importance to what follows, because the literature focuses on 3D rather than 2D, and because classical potential theory is rarely used in the study of heat conduction. We will work with complex variables to take advantage of strong analytical tools.

Classical developments of the boundary derivatives, in real variables, are given by: Poincaré [11], Kellogg [4], MacMillan [5], and Stakgold [2].

4.1 The Simple-Layer Normal Derivative as a Cauchy Integral.

To evaluate the normal derivative of temperature at the contour, we consider the directional derivative of Eq. (32) in the complex z-plane, following Xu et al. [12]. We rename integrated point ξ as $z \in C$ ⁠, and consider the fixed point $x \notin C$ ⁠, near s. We can differentiate at under the integral with respect to x in the direction $n_{s}$ ⁠. The angle between $n_{s}$ and the real axis is $ϕ_{0}$ (Fig. 5). Here, x, s, and z are points in the complex plane.

Fig. 5

View large Download slide

Coordinates for directional derivatives in the complex plane

From Sec. 2.2

D_{n_{s}} (z, x) = \frac{1}{2 π} \frac{\partial}{\partial n_{s}} log \frac{1}{| z - x |} = \frac{cos (z - x, n_{s})}{2 π | z - x |}

(38a)

= R e {\frac{e^{i [ϕ_{0} - \arg (z - x)]}}{2 π | z - x |}} = R e {\frac{e^{i ϕ_{0}}}{2 π (z - x)}}

(38b)

where the angle for the cosine in Eq. (38a) is $ϕ_{0} - \arg (z - x)$ ⁠.

With reference to Fig. 5,

d z = e^{i α} | d z |

with

ϕ = α (z) - π / 2

⁠. Continuing in the complex plane

\frac{\partial θ}{\partial n_{s}} = \int_{C} a (z) \frac{cos (z - x, n_{s})}{2 π | z - x |} | d z |

(39a)

= \int_{C} a (z) R e {\frac{e^{i ϕ_{0}}}{2 π (z - x)}} | d z |

(39b)

= \int_{C} R e {a (z) \frac{e^{i α_{0}} e^{- i α} d z}{2 π i (z - x)}}

(39c)

= R e {e^{i α_{0}} d (x)}

(39d)

for d(x) defined as

d (x) = \frac{1}{2 π i} \int_{C} \frac{a (z) e^{- i α}}{z - x} d z

(40)

where a(z) is real valued.

4.2 Sokhotski–Plemelj Theorem.

This theorem is well-known in complex analysis and widely used in physics. Following Muskhelishvili’s development [13], consider the contour integral:

Φ (x) = \frac{1}{2 π i} \int_{C} \frac{φ (z) d z}{z - x}

(41)

where $φ (z)$ is defined on a smooth contour $C$ and bounded except at a finite number of points. The integral defines an analytic function for $x \notin C$ ⁠; further, $Φ = O (| x |^{- 1})$ as $| x | \to \infty$ ⁠. When $x = s \in C$ ⁠, the integral exists only in the Cauchy principal value sense.

The Sokhotski–Plemelj theorem relates the limiting boundary values of

Φ

at a point s on

C

to the principal value of the integral (/

\int

⁠)

\lim_{x \to s^{\pm}} Φ (x) = \mp \frac{1}{2} φ (x) + \frac{1}{2 π i} / \int_{C} \frac{φ (z) d z}{z - x}

(42)

where s⁺ indicates approaches from outside the contour and s^– indicates approach from inside. The analysis requires that

φ

be Hölder continuous

| φ (z_{1}) - φ (z_{2}) | \leq C | z_{1} - z_{2} |^{μ}

(43)

for all z in the domain of $φ$ and non-negative real constants C and μ [13].

4.3 Plemelj Formulæ for the Normal Derivative.

Equation is a Cauchy integral, and its value is discontinuous across

C

⁠. Assuming that the contour is smooth at s, Eq. gives

\lim_{x \to s^{\pm}} d (x) = \mp \frac{1}{2} a (s) e^{- i α_{0}} + \frac{1}{2 π i} / \int_{C} \frac{a (z) e^{- i α} d z}{z - s}

(44)

Therefore

\lim_{x \to s^{\pm}} \frac{\partial θ}{\partial n_{s}} = \mp \frac{a (s)}{2} + R e {\frac{e^{i α_{0}}}{2 π i} / \int_{C} \frac{a (z) e^{- i α} d z}{z - s}}

(45a)

and by comparison to Eq. , the real variable expression, for

ξ

and s on

C

⁠, is

\lim_{x \to s^{\pm}} \frac{\partial θ}{\partial n_{s}} = \underset{source gradient}{\underset{︸}{\mp \frac{1}{2} a (s)}} + \underset{induced gradient}{\underset{︸}{\int_{C} a (ξ) \frac{cos (ξ - s, n_{s})}{2 π | ξ - s |} d l_{ξ}}}

(45b)

As noted in Sec. 2.3, this kernel does not have a singularity at $ξ = s$ ⁠, so a principal value calculation is not required. In cases such as corner points, a(z) may be singular, but those singularities are generally integrable (see Secs. 6.1 and 7).

The normal derivative, Eq. (45b), includes a source term and an integral representing gradients in the normal direction at $s$ induced by sources at other points $ξ$ around the contour (cf. Eq. (18)). The source a(s) is bisected by the surface. Half the source power exits one side of the surface and half exits the opposite side, in the opposite direction.² The induced gradient term depends on the distance to the source point, $| ξ - s |$ ⁠, and the orientation of $ξ - s$ relative to the fixed point’s outward normal vector, $n_{s}$ ⁠. The integral contains the contributions of the entire boundary. This induced gradient has the same sign and is on either side of the boundary, $x \to s^{\pm}$ ⁠. The induced temperature gradient indicates the presence of an induced heat flux.

4.3.1 More About Induced Heat Flux.

The properties of the induced heat flux are needed in Sec. 5.1, and so further discussion is helpful.

The heat conduction problem has a strong analogy in electrostatics: a surface heat source density, a(s), is analogous to a surface charge density, and a temperature gradient (or heat flux) is analogous to an electric potential gradient (or electric field). In electrostatics, a surface charge distribution may be induced by externally generated electric fields, through the redistribution of charge over the conducting surface (see Kellogg [4] and Jackson [14]). If the surface is perfectly conducting, it remains isopotential even with a nonuniform distribution of surface charge. If the conductor is not grounded, no net additional charge is acquired when the surface charge is redistributed. Thus, the integral of the induced charge density over the surface is zero.

Similarly, a surface heat flux distribution may be induced on a heat conductor by heat sources elsewhere. This induced heat flux adds no additional energy to the system, and consequently, the induced flux integrates to zero over the conductor’s surface. Further, a perfectly conducting surface has no resistance to heat flow and remains isothermal even with a nonuniform surface heat flux. If this conductor is simply one isothermal boundary of a larger object, say with two nonconducting boundaries separating two isothermal boundaries, the induced flux on each isothermal boundary will integrate to zero separately because the induced redistribution of flux is confined to each isothermal boundary separately.

For visualization, Fig. 6 shows the temperature and magnitude of heat flux outside a perfectly conducting disk placed between isothermal vertical plates. For this disk, $a (ξ) = 0$ (the boundary contains no heat source distribution), and so Eq. (45b) shows that the heat flux is composed only of induced heat flux. The induced heat flux has high magnitude (yellow) at some points and tends to zero (blue) at other points. The integral of this induced heat flux around the conductor is zero, consistent with steady energy conservation. The conductor’s temperature depends on its shape and position relative to the external source and sink but is easily seen to be zero in the present case. If the disk were moved off center, the temperature could be found by solving Eq. (26) with $a (ξ) = 0$ on the disk.

Fig. 6

Perfectly conducting disk between vertical plates at θ=±1 with adiabatic horizontal sides. MATLAB FEM simulation for disk of diameter 0.5 inside 2 × 2 square. (a) Temperature distribution; (b) magnitude of flux vector (lightest yellow is largest flux, darkest blue is zero flux); plus/minus signs indicate where the induced heat flux entering the disk is more positive and more negative.

View large Download slide

Perfectly conducting disk between vertical plates at $θ = \pm 1$ with adiabatic horizontal sides. MATLAB FEM simulation for disk of diameter 0.5 inside 2 × 2 square. (a) Temperature distribution; (b) magnitude of flux vector (lightest yellow is largest flux, darkest blue is zero flux); plus/minus signs indicate where the induced heat flux entering the disk is more positive and more negative.

To fix the perfect conductor to a specific temperature, we can imagine it to be grounded in to a heat reservoir at that temperature. The reservoir adds or removes heat as needed to maintain the temperature. This heat flow is distributed over the surface of the conductor and appears mathematically as the heat source distribution, a(s).

The preceding discussion was for perfect conductors. In the case of nonconductors, external heat sources (or electric charges) induce dipoles on the surface of the nonconductor. These dipoles balance out the local potential gradients (temperature gradients or electric fields) caused by external sources. In the electrostatic case, the net charge density induced is zero, even locally, because dipoles have no net charge. In heat conduction, the induced dipoles add no net energy; and, they cancel the heat flux that sources elsewhere would otherwise induce, leaving the surface adiabatic. A nonconducting object will develop internal temperature (or potential) gradients. As a result, the temperature (or potential) is discontinuous across a nonconducting boundary (cf. Eq. (49)).

We return induced flux in Sec. 5.1.

4.4 Special Cases of the Normal Derivative.

In some special cases, the integral in Eq. (45b) can vanish, as the following examples show.

4.4.1 Straight Lines and Parallel Plates.

For a straight-line segment of $C$ with outward normal n, $cos (ξ - s, n) = 0$ if $ξ$ and $s$ are both on that segment. The reason is that a dipole creates no potential perpendicular to its axis. If we generalize to hot and cold parallel plates of infinite length, the integral along the plate containing s is therefore zero, and that along the opposing plate is easily shown to vanish in the limit of an infinitely long plate.

4.4.2 Circular Disk.

For the case of a circular disk of radius $R$ ⁠, some trigonometry shows that $cos (ξ - s, n_{s}) / | ξ - s | = - 1 / 2 R$ ⁠, and with Eq. (35), the integral in Eq. (45b) will vanish at every s, leaving $\partial θ / \partial n_{s} = - \partial θ^{e} / \partial n_{s}$ ⁠.

In this unusual case, the isothermal boundary has the same temperature and temperature gradient on both the inside and outside. This behavior is usually not seen in other geometries. (We establish a similar result for the adiabatic boundary of a disk in Sec. 5.2.)

4.5 Plemelj Formulæ for the Integral of the Normal Derivative (Double-Layer Potential).

In Sec. 5, we will need the double-layer potential and its normal derivative. The double-layer potential is an integral of the normal derivative of

E (x | z)

⁠, now with a normal

n_{z}

that moves along the curve during integration in z (Fig. 5). We proceed as in the simple-layer analysis

D_{n_{z}} (x, z) = \frac{1}{2 π} \frac{\partial}{\partial n_{z}} log \frac{1}{| x - z |} = \frac{cos (x - z, n_{z})}{2 π | x - z |}

(46a)

= R e {\frac{e^{i [\arg (x - z) - ϕ]}}{2 π | x - z |}} = - R e {\frac{e^{- i ϕ}}{2 π (z - x)}}

(46b)

since

\arg (x - z) = π - \arg (z - x)

and

e^{i π} = - 1

⁠. Continuing

θ_{d} (x) = \int_{C} b (z) D_{n_{z}} (x, z) d z

(47a)

= - \int_{C} b (z) R e {\frac{e^{- i ϕ}}{2 π (z - x)}} | d z |

(47b)

= - \int_{C} R e {b (z) \frac{e^{+ i π / 2} d z}{2 π (z - x)}} = R e {g (x)}

(47c)

for real-valued b(z),

- e^{i π / 2} = 1 / i

⁠, and g(x) defined as

g (x) = \frac{1}{2 π i} \int_{C} \frac{b (z)}{z - x} d z

(48)

Equation is a Cauchy integral, and it is discontinuous across

C

⁠. Assuming that the contour is smooth at s, and with b(s) real valued, Eq. gives

\begin{matrix} \lim_{x \to s^{\pm}} θ_{d} (x) = \mp \frac{b (s)}{2} + R e {\frac{1}{2 π i} / \int_{C} \frac{b (z) d z}{z - s}} \\ = \mp \frac{b (s)}{2} + \int_{C} b (ξ) \frac{cos (s - ξ, n_{ξ})}{2 π | s - ξ |} d l_{ξ} \end{matrix}

(49)

for $n_{ξ}$ the outward normal at the moving point ξ.

The double-layer potential θ_d is discontinuous at $C$ ⁠, but the normal derivative $\partial θ_{d} / \partial n_{x}$ has the same limit from either side (see Poincaré [11], Kellogg [4], MacMillan [5], Günter [15]). Kellogg shows that, if the double layer potential b(z) has continuous second derivatives, the difference of the normal derivatives to either side tends to zero and that the limits of the derivatives exist at the surface; Günter provides a similar result.

Discontinuities of the higher derivatives of the simple and double-layer potentials are discussed by Schmidt [16] and Hinrichsen [17].

5 Mixed Boundary Conditions as Simple and Double Layer Potentials

We now consider the configuration of primary interest, Fig. 1. The boundary

C

consists of two isothermal segments,

Γ^{h}

and

Γ^{c}

⁠, separated by two adiabatic segments,

Γ_{1}^{a}

and

Γ_{2}^{a}

⁠. Let

Γ^{a} = Γ_{1}^{a} \cup Γ_{2}^{a}

and

Γ^{i} = Γ^{h} \cup Γ^{c}

⁠. The initial results follow from Sec. 3.1. For the exterior region, we have as before

\int_{C} [- E \frac{\partial θ^{e}}{\partial n_{ξ}} - θ^{e} \frac{\partial E}{\partial n_{ξ}}] d l_{ξ} + θ^{\infty} = {\begin{matrix} 0 & x \in R \\ θ^{e} (x) & x \in R^{e} \end{matrix}

(50)

For the interior region, we have

\int_{C} [E \frac{\partial θ}{\partial n_{ξ}} - θ \frac{\partial E}{\partial n_{ξ}}] d l_{ξ} = {\begin{matrix} θ (x) & x \in R \\ 0 & x \in R^{e} \end{matrix}

(51)

By adding the interior and exterior results, we get

\int_{C} [E (\frac{\partial θ}{\partial n_{ξ}} - \frac{\partial θ^{e}}{\partial n_{ξ}}) - (θ - θ^{e}) \frac{\partial E}{\partial n_{ξ}}] d l_{ξ} + θ^{\infty} = {\begin{matrix} θ (x) & x \in R \\ θ^{e} (x) & x \in R^{e} \end{matrix}

(52)

For

ξ \in Γ^{a}

the first term inside the integral is zero, and for

ξ \in Γ^{i}

the second term inside the integral is zero. We define the single-layer potential density as in Eq. , so that

a (ξ) = {\begin{matrix} \partial θ / \partial n_{ξ} - \partial θ^{e} / \partial n_{ξ} & ξ \in Γ^{i} \\ 0 & ξ \in Γ^{a} \end{matrix}

(53)

We define the double-layer potential density as

b \equiv θ^{e} - θ

⁠, so that

b (ξ) = {\begin{matrix} 0 & ξ \in Γ^{i} \\ θ^{e} - θ & ξ \in Γ^{a} \end{matrix}

(54)

Collecting all this

\int_{Γ^{i}} a (ξ) E (x | ξ) d l_{ξ} + \int_{Γ^{a}} b (ξ) \frac{\partial E (x | ξ)}{\partial n_{ξ}} d l_{ξ} + θ^{\infty} = {\begin{matrix} θ (x) & x \in R \\ θ^{e} (x) & x \in R^{e} \end{matrix}

(55)

As before, these results show that $θ^{\infty}$ is imposed by the far field in 2D. The relationship of $θ^{\infty}$ to the average temperature on $C$ is discussed in Appendix A.

We now form boundary conditions on

C

⁠. For the isothermal boundaries, with

s \in Γ^{i}

⁠, the first integral in Eq. is continuous, but the second is not. We apply the Plemelj formulæ to the second integral with Eq.

\begin{matrix} \int_{Γ^{i}} a (ξ) E (s | ξ) d l_{ξ} \\ \underset{= 0 on Γ^{i}}{\underset{︸}{\mp \frac{1}{2} b (s)}} + / \int_{Γ^{a}} b (ξ) \frac{\partial E (s | ξ)}{\partial n_{ξ}} d l_{ξ} + θ^{\infty} = θ (s) \end{matrix}

(56)

For the adiabatic boundaries,

s \in Γ^{a}

⁠, with Eqs. , , and

\begin{matrix} \underset{= 0 on Γ^{a}}{\underset{︸}{\mp \frac{1}{2} a (s)}} + / \int_{Γ^{i}} a (ξ) \frac{\partial E (ξ | s)}{\partial n_{s}} d l_{ξ} \\ + \lim_{x \to s^{\pm}} \frac{\partial}{\partial n_{x}} \int_{Γ^{a}} b (ξ) \frac{\partial E (x | ξ)}{\partial n_{ξ}} d l_{ξ} = 0 \end{matrix}

(57)

The double layer potential has the same normal derivative for either $x \to s^{\pm}$ (see Sec. 4.5), as is clearly required to satisfy Eq. (57).

Equation (57) shows that, locally on the adiabatic edges, the gradient of the dipole distribution exactly cancels the gradient imposed by the sources on the isothermal boundaries. We may think of the dipole distributions on the adiabatic edges as being induced by the potential gradients of the sources on the isothermal edges.

In Eq. (56), the kernel of the integral on $Γ^{i}$ has a logarithmic singularity, while the second integral has a nonsingular kernel (Sec. 2.3).³ In Eq. (57), the first integral also has a nonsingular kernel. The second term in Eq. (57), in general, is hypersingular; but as previously noted, under appropriate conditions on b [4,15], this limits of this normal derivative exist and are equal inside and outside $C$ ⁠.

The temperature gradient in the direction outwardly normal to the isothermal contours,

s \in Γ^{i}

⁠, is

\begin{matrix} \frac{\partial θ}{\partial n_{s}} |_{s^{\pm}} = \mp \frac{1}{2} a (s) + \int_{Γ^{i}} a (ξ) \frac{\partial E (ξ | s)}{\partial n_{s}} d l_{ξ} \\ + \frac{\partial}{\partial n_{x}} \int_{Γ^{a}} b (ξ) \frac{\partial E (x | ξ)}{\partial n_{ξ}} d l_{ξ} |_{x \to s^{\pm}} \end{matrix}

(58)

The normal derivative of the second integral has no singularity for $s \in Γ^{i}$ ⁠, in contrast to $s \in Γ^{a}$ ⁠. Integration of this equation over all of $C$ will show that the first and second terms are zero, as a result of the zero net heat transfer requirement. Thus, the other two terms, integrated over $C$ ⁠, will sum to zero. Further, since these last two terms cancel locally on the adiabatic boundaries, by Eq. (57), we see that the integral of the two terms over the hot boundary is equal and opposite to the integral of the pair on the cold boundary (but note the stronger conclusion given in Sec. 5.1).

The temperatures on the adiabatic edges,

s \in Γ^{a}

⁠, are

θ (s) = θ^{\infty} + \int_{Γ^{i}} a (ξ) E (s | ξ) d l_{ξ} \mp \frac{1}{2} b (s) + \int_{Γ^{a}} b (ξ) \frac{\partial E (s | ξ)}{\partial n_{ξ}} d l_{ξ}

(59)

Equations (58) and (59) have two very important implications:

On isothermal boundaries, the interior and exterior temperatures are by definition equal. Equation (58) shows that the interior and exterior heat fluxes need not have locally equal magnitudes on the isothermal boundaries.
On adiabatic boundaries, the interior and exterior fluxes are both zero. Equation (59) shows that the interior and exterior temperatures need not be locally equal on the adiabatic boundaries, as expected for a dipole boundary condition.

For convenience, the integral kernels are summarized in Table 1.

Table 1

Integral kernels

Meaning	Primary symbol	Full expression	Equivalent symbol
Potential at s caused by unit source density at ξ	$E (s \| ξ)$	$\frac{1}{2 π} log \frac{1}{\| s - ξ \|}$	$E (ξ \| s)$
Outward normal gradient at s caused by unit source density at ξ	$\frac{\partial E (ξ \| s)}{\partial n_{s}}$	$\frac{cos (ξ - s, n_{s})}{2 π \| ξ - s \|}$	$D_{n_{s}} (ξ \| s)$
Potential at s caused by unit outward-normal dipole density at ξ	$\frac{\partial E (s \| ξ)}{\partial n_{ξ}}$	$\frac{cos (s - ξ, n_{ξ})}{2 π \| s - ξ \|}$	$D_{n_{ξ}} (s \| ξ)$

Meaning	Primary symbol	Full expression	Equivalent symbol
Potential at s caused by unit source density at ξ	$E (s \| ξ)$	$\frac{1}{2 π} log \frac{1}{\| s - ξ \|}$	$E (ξ \| s)$
Outward normal gradient at s caused by unit source density at ξ	$\frac{\partial E (ξ \| s)}{\partial n_{s}}$	$\frac{cos (ξ - s, n_{s})}{2 π \| ξ - s \|}$	$D_{n_{s}} (ξ \| s)$
Potential at s caused by unit outward-normal dipole density at ξ	$\frac{\partial E (s \| ξ)}{\partial n_{ξ}}$	$\frac{cos (s - ξ, n_{ξ})}{2 π \| s - ξ \|}$	$D_{n_{ξ}} (s \| ξ)$

The cosine of the angle between vectors $a$ and $b$ is denoted by $cos (a, b)$ ⁠.

5.1 Induced Heat Flux.

Equation provides the local temperature gradient (taken in the outward normal direction) on either side of an isothermal edge. The gradient consists of three parts: (i) one-half of the heat source density at s (half to either side of the boundary); (ii) the gradient induced at s by the sources at all other points on the boundary; and (iii) the gradient induced at s by dipoles at all other points on the boundary. As discussed in Sec. 4.3.1, the net induced heat flux on each conducting (isothermal) boundary is zero. Therefore, for the hot isothermal boundary

\begin{array}{l} \int_{Γ_{s}^{h}} \int_{Γ_{ξ}^{i}} a (ξ) \frac{\partial E (ξ | s)}{\partial n_{s}} d l_{ξ} d l_{s} \\ + \int_{Γ_{s}^{h}} \frac{\partial}{\partial n_{x}} \int_{Γ_{ξ}^{a}} b (ξ) \frac{\partial E (x | ξ)}{\partial n_{ξ}} d l_{ξ} |_{x \to s} d l_{s} = 0 for Γ^{h} isothermal \end{array}

(60)

The limit should be continuous, as the integral is nonsingular for $x \to s \in Γ^{h}$ and $ξ \in Γ^{a}$ ⁠, except at the point where $Γ^{h}$ and $Γ^{a}$ intersect. The same conclusion applies to the net induced flux on the cold isothermal boundary, $Γ^{c}$ ⁠.

We may pair this result with that previously established for the nonconducting (adiabatic) edges, Eq.

\begin{matrix} \int_{Γ^{i}} a (ξ) \frac{\partial E (ξ | s)}{\partial n_{s}} d l_{ξ} \\ + \frac{\partial}{\partial n_{x}} \int_{Γ^{a}} b (ξ) \frac{\partial E (x | ξ)}{\partial n_{ξ}} d l_{ξ} |_{x \to s^{\pm}} = 0 for s \in Γ^{a} adiabatic \end{matrix}

(61)

As noted in Sec. 4.5, the normal derivative of second integral has the same limit for $x \to s^{\pm}$ ⁠, but is undefined for $x = s \in Γ^{a}$ ⁠.

5.2 Mixed Boundary Conditions on a Disk.

Consider a disk with mixed boundary conditions (Fig. 7): the boundary in the first quadrant is isothermal at $θ (s) = 1$ ⁠, the boundary in the third quadrant is isothermal at $θ (s) = - 1$ ⁠, and the other two boundaries are adiabatic. This problem was solved by Schwarz in 1869 using conformal mapping [18].

Fig. 7

View large Download slide

Isotherms and adiabats inside a disk with two isothermal edges and two adiabatic edges (black) [1]

For a circular disk, as in Sec. 4.4.2, the kernel of the integral on $a (ξ)$ in Eq. (57) reduces to a constant and the integral vanishes. The integral in $b (ξ)$ retains an x dependence, but the heat flux boundary condition on the adiabatic sides is satisfied by $b (ξ) = 0$ ⁠, and Eq. (59) then indicates equal interior and exterior temperatures on the adiabatic edges. Further, with $b (ξ) = 0$ ⁠, Eq. (56) reduces to a Fredholm problem for the potential on the isothermal sections of the boundary, with $a (ξ) = 0$ on the adiabatic sides. In addition, in this case, the heat flux outside and inside the isothermal edges are equal and opposite, given by $\mp a (s) / 2$ according to Eq. (58). All these findings are consistent with the solution by Poisson’s formula given in Ref. [1].

5.2.1 Fredholm Equation for Schwarz’s Problem.

We may use these ideas to set a Fredholm problem for

a (ψ)

⁠, where ψ is the angle around the perimeter of the a disk measured from the x-axis. We know that

b (ψ) = 0

and that the heat flux boundary conditions on the adiabatic edges are met. So, Eq. reduces to the following equation on the isothermal edges of the unit disk:

\int_{Γ^{i}} \frac{a (ψ_{0})}{4 π} log \frac{1}{2 [1 - cos (ψ - ψ_{0})]} d ψ_{0} = {\begin{matrix} + 1 & 0 < ψ < π / 2 \\ - 1 & π < ψ < 3 π / 2 \end{matrix}

(62)

where we obtain

E (s | ξ) = E (1, ψ | 1, ψ_{0})

from Eq. , and

Γ^{i} = Γ^{h} \cup Γ^{c}

⁠. Continuing

\begin{matrix} \int_{0}^{π / 2} \frac{a (ψ_{0})}{4 π} log \frac{1}{2 [1 - cos (ψ - ψ_{0})]} d ψ_{0} \\ + \int_{π}^{3 π / 2} \frac{a (ψ_{0})}{4 π} log \frac{1}{2 [1 - cos (ψ - ψ_{0})]} d ψ_{0} \\ = \int_{0}^{π / 2} \frac{a (ψ_{0})}{4 π} log \frac{1}{2 [1 - cos (ψ - ψ_{0})]} d ψ_{0} \\ + \int_{0}^{π / 2} \frac{a (ψ_{0} + π)}{4 π} log \frac{1}{2 [1 - cos (ψ - ψ_{0} + π)]} d ψ_{0} \\ = {\begin{matrix} + 1 & 0 < ψ < π / 2 \\ - 1 & π < ψ < 3 π / 2 \end{matrix} \end{matrix}

(63)

From symmetry, and because the cold edge is a sink, we expect that

a (ψ_{0}) = - a (ψ_{0} + π)

⁠. If we define

4 π K (ψ, ψ_{0}) \equiv log [2 - 2 cos (ψ - ψ_{0})] - log [2 - 2 cos (ψ - ψ_{0} + π)]

(64)

we may observe that

K (ψ, ψ_{0}) = - K (ψ + π, ψ_{0})

⁠. Hence, it will be sufficient to solve the following first-kind Fredholm problem on the hot edge:

\int_{0}^{π / 2} a (ψ_{0}) K (ψ, ψ_{0}) d ψ_{0} = 1 for 0 < ψ < π / 2

(65)

Note that the kernel is weakly singular at $ψ = ψ_{0}$ ⁠. We return to this equation in Sec. 7.1.

6 Singularities at Corners and Dirichlet-to-Neumann Transitions

The FEM results in Sec. 7 contain singularities in the temperature distributions and simple-layer potential densities. This section provides the theoretical form of those singularities.

6.1 Corner Singularities.

Consider a corner or wedge shaped region, and place the origin of polar coordinates

(r, ϕ)

at the vertex. If the included angle is β, and we have one edge that is isothermal and one that is adiabatic, the solutions of the Laplace equation, from Eq. , will be [19]

θ (r, ϕ) = θ_{0} + \sum_{m = 1}^{\infty} c_{m} r^{η_{m}} cos (η_{m} ϕ) for η_{m} = \frac{π}{2 β} (2 m - 1)

(66)

where the adiabatic edge is at

ϕ = 0

and the isothermal edge is at

ϕ = β

⁠. The derivative of θ normal to a line of constant

ϕ

is

\frac{\partial θ}{\partial n} = \frac{1}{r} \frac{\partial θ}{\partial ϕ} = - \sum_{m = 1}^{\infty} c_{m} η_{m} r^{η_{m} - 1} sin (η_{m} ϕ)

(67)

and for

ϕ = β

⁠, the normal derivative becomes

\frac{\partial θ}{\partial n} |_{ϕ = β} = - \sum_{m = 1}^{\infty} {(- 1)}^{m + 1} c_{m} η_{m} r^{η_{m} - 1}

(68)

When $η_{1} < 1$ (⁠ $β > π / 2$ ⁠), the derivative contains an integrable singularity as $r \to 0$ ⁠. For smaller β, the derivative is nonsingular.

For the specific case of the interior of a square corner,

β = π / 2

and

η_{m} = (2 m - 1)

\frac{\partial θ}{\partial n} |_{ϕ = π / 2} = - \sum_{m = 1}^{\infty} {(- 1)}^{m + 1} (2 m - 1) c_{m} r^{2 (m - 1)} = - 1

(69)

where the final equality (⁠ $c_{1} = 1$ ⁠, c_m = 0 for m > 1) is from my knowledge of the solution: $θ (r, ϕ) = 1 - r cos ϕ$ (note that n is outward normal here).

For the region exterior to a square corner,

β = 3 π / 2

and

η_{m} = (2 m - 1) / 3

\frac{\partial θ}{\partial n} |_{ϕ = 3 π / 2} = - \sum_{m = 1}^{\infty} \frac{(2 m - 1) c_{m}}{3} {(- 1)}^{m + 1} r^{2 (m - 2) / 3}

(70a)

= - \frac{c_{1}}{3 r^{2 / 3}} + c_{2} - \frac{5 c_{3} r^{2 / 3}}{3} + \dots

(70b)

We may also consider a corner with a continuous Dirichlet boundary condition,

θ = θ_{0}

(see Fig. 10(b) in Ref. [1]). Now

θ (r, ϕ) = θ_{0} + \sum_{m = 1}^{\infty} c_{m} r^{η_{m}} sin (η_{m} ϕ) for η_{m} = \frac{π m}{β}

(71)

so that for

ϕ = β

\frac{\partial θ}{\partial n} |_{ϕ = β} = - \sum_{m = 1}^{\infty} {(- 1)}^{m} c_{m} η_{m} r^{η_{m} - 1}

(72)

For the interior of a square corner,

β = π / 2

and

η_{1} - 1 = 1

⁠; for the exterior,

β = 3 π / 2

and

η_{1} - 1 = - 1 / 3

⁠, which is again an integrable singularity.

6.2 Change From Dirichlet to Neumann Boundary Condition on a Straight Section.

For a straight segment that goes from Dirichlet to Neumann at r = 0 (e.g., locally on the rim of a disk),

β = π

and

η_{m} = (m - 1 / 2)

⁠, so that

\frac{\partial θ}{\partial n} |_{ϕ = π} = - \sum_{m = 1}^{\infty} \frac{(2 m - 1) c_{m}}{2} {(- 1)}^{m + 1} r^{m - 3 / 2}

(73a)

= - \frac{c_{1}}{2 r^{1 / 2}} + \frac{3 c_{2} r^{1 / 2}}{2} - \frac{5 c_{3} r^{3 / 2}}{2} + \dots

(73b)

Note that β will have the same value inside and outside this segment, so that locally the interior and exterior gradients will have the same dependence on r. Likewise, from Eq. (66), on the adiabatic segment (⁠ $ϕ = 0$ ⁠) the temperature distribution inside and outside will be the same locally.

An extensive literature deals with boundary integrals and Fredholm problems near corners [8,19–21].

7 Finite Element Method Studies of Potential and Induced Flux

In this section, we use FEM calculations to independently explore the conclusions of Secs. 5 and 6 for mixed boundary conditions on a disk and a square, with a focus on the integrals of the simple and double layer densities.

7.1 Disk With Isothermal and Adiabatic Edges.

We may use MATLAB’s FEM codes to analyze the quarter-sectored disk of radius $R = 1$ (Schwarz’s problem, Fig. 7). The resultant temperature gradient along the hot isothermal side is shown in Fig. 8. Convergence of the FEM model is discussed in Appendix B.

Fig. 8

Normal gradient, ∂θ/∂n, along interior of isothermal disk edge in upper left quadrant. The calculated shape factor is Si=0.9986. Least-squares curve lies on top of FEM curve.

View large Download slide

Normal gradient, $\partial θ / \partial n$ ⁠, along interior of isothermal disk edge in upper left quadrant. The calculated shape factor is $S^{i} = 0.9986$ ⁠. Least-squares curve lies on top of FEM curve.

The singularity at each end of the isothermal section, from Eq. , is

O (r^{- 1 / 2})

⁠; and if the angle around the disk is ψ measured from the x-axis, we may set

r = R ψ

and

r = R (π / 2 - ψ)

for the two corners of the isothermal edge in upper left quadrant of the figure. The computed temperature gradient may be fitted to the analytical form suggested by Eq. , using MATLAB’s fitting algorithms. A fit that accounts for the singularity at ψ = 0 and

ψ = π / 2

⁠, each to three terms, is

\begin{matrix} \frac{\partial θ}{\partial n} = \frac{1}{R} \frac{\partial θ}{\partial ψ} ≐ 1.069 (\frac{1}{2}) [\frac{1}{ψ^{1 / 2}} + \frac{1}{{(π / 2 - ψ)}^{1 / 2}}] \\ - 0.1331 (\frac{3}{2}) [ψ^{1 / 2} + {(π / 2 - ψ)}^{1 / 2}] \\ - 0.02553 (\frac{5}{2}) [ψ^{3 / 2} + {(π / 2 - ψ)}^{3 / 2}] \end{matrix}

(74)

The root-mean-square error of the fit relative to the FEM data is 0.0041. This fit is plotted as a dotted curve in Fig. 8, but it lies directly atop the FEM result.

Recall from Sec. 5.2 that in the special case of a circular disk, the normal fluxes are equal and opposite on the interior and exterior of the disk’s isothermal edges.

The average of Eq. on

0 < ψ < π / 2

is

\bar{\frac{\partial θ}{\partial n}} = - \frac{2}{π} \int_{0}^{π / 2} \frac{\partial θ}{\partial n} d ψ

(75)

\begin{matrix} = (1.069) \frac{4}{π} \sqrt{\frac{π}{2}} - (0.1331) \frac{4}{π} {(\frac{π}{2})}^{3 / 2} - (0.02552) \frac{4}{π} {(\frac{π}{2})}^{5 / 2} \\ = (1.069) (2) \sqrt{\frac{2}{π}} - (0.1331) (2) \sqrt{\frac{π}{2}} - (0.02552) (2) {(\frac{π}{2})}^{3 / 2} \\ = 1.2718 \end{matrix}

(76)

For a shape factor $S^{i} = 1$ [1] and $Δ θ = 2$ ⁠, the theoretical average gradient is $2 (2 / π) = 4 / π = 1.2732$ ⁠, indicated by a horizontal dotted line in Fig. 8. The integrated result is lower by only 0.11%.

The computed temperature distribution along the adiabatic edge of the disk,

θ_{ae}

⁠, is shown in Fig. 9. The slope at the edges is consistent with Eq. for

β = π

and

r = R (π / 2 - ψ)

or

r = - R ψ

⁠, so that

θ (r, 0) = θ_{ae} (ψ)

⁠, for

- π / 2 < ψ < 0

⁠. With MATLAB curve-fitting routines, we again find coefficients that fit a polynomial to the FEM results (plotted as a dotted curve atop the FEM result in Fig. 9)

\begin{matrix} θ_{ae} (ψ) ≐ 1.071 [{(π / 2 - ψ)}^{1 / 2} - ψ^{1 / 2}] \\ - 0.1357 [{(π / 2 - ψ)}^{3 / 2} - ψ^{3 / 2}] \\ - 0.0244 [{(π / 2 - ψ)}^{5 / 2} - ψ^{5 / 2}] \end{matrix}

(77)

Fig. 9

View large Download slide

Temperature $θ_{ae}$ along adiabatic disk edge in lower left quadrant. Least-squares curve lies on top of FEM curve.

The root-mean-square error of the fit is $6.4266 \times 10^{- 5}$ ⁠. The corner solution fits the FEM result very well, despite the curvature of the actual boundary. Relative to the corner solution, Eqs. (66) and (67), these fits suggest that $c_{1} ≐ 1.070, c_{2} ≐ - 0.134$ ⁠, and $c_{3} ≐ - 0.025$ ⁠.

The single-layer potential is defined as $a = \partial θ / \partial n - \partial θ^{e} / \partial n$ (Eq. (53)). For the disk, from Sec. 5.2, we know that $\partial θ^{e} / \partial n = - \partial θ / \partial n$ ⁠, and so the single-layer potential $a (ψ) = 2 \partial θ / \partial n$ ⁠.⁴ Numerical evaluation of the integral in Eq. (65), in MATLAB, using the $a (ψ)$ fitted above, gives the expected constant value of 1 to high precision: the average over ψ of the integral is 1.0001 with a standard deviation of 0.00023.

7.2 Square With Isothermal and Adiabatic Sides.

Suppose that one side of a square has temperature θ = 1, that the opposing side has $θ = - 1$ ⁠, and that the two separate sides are adiabatic. The interior gradient is uniform and equal for the hot and cold sides, from the easy solution for the interior of the square. The integrals of the exterior heat flux over the hot side and over the cold side must be equal and opposite for energy conservation.

We may observe from Fig. 10 that the adiabats, and thus local heat fluxes, inside and outside the square are not the same. The exterior flux is low in the center and increases greatly as the outside corners are approached. Further, the temperature variation along the outside adiabatic boundary is not the same as that along the inside adiabatic boundary: the interior has a linear temperature distribution whereas the exterior changes steeply near the corners.

Fig. 10

View large Download slide

Heat conduction inside and outside a square with adiabatic horizontal sides, showing isotherms and adiabats [1]: (a) interior region and (b) exterior region

The exterior normal temperature gradient along an isothermal edge of the square can be computed with MATLAB, as shown in Fig. 11. Convergence of the FEM model is discussed in Appendix B. Here, the numerical square is 2 × 2 with a temperature difference of $Δ θ = 2$ ⁠. The shape factor for the exterior of a square is $S^{e} = 1$ [1]; and, with the thermal conductivity normalized to unity, the dimensionless heat transfer rate is ${\dot{Q}}_{ext} = 2$ ⁠.

Fig. 11

Exterior normal gradient, −∂θe/∂n, along hot isothermal edge of square. The interior normal gradient, ∂θ/∂n=1 for all y. Least-squares curve lies on top of FEM curve.

View large Download slide

Exterior normal gradient, $- \partial θ^{e} / \partial n$ ⁠, along hot isothermal edge of square. The interior normal gradient, $\partial θ / \partial n = 1$ for all y. Least-squares curve lies on top of FEM curve.

7.2.1 Fitting Heat Flux and Temperature With Corner Theory.

Note that

\partial θ^{e} / \partial n < 0

because the heat flows in the direction of the outward normal vector n. Corner theory, Eq. , superposed for the two corners, suggests fitting with this approximate form, where

v \in [- 1, 1]

\begin{matrix} - \partial θ^{e} / \partial n \approx c_{1} (1 / 3) [{(v + 1)}^{- 2 / 3} + {(1 - v)}^{- 2 / 3}] \\ + c_{3} (5 / 3) [{(v + 1)}^{2 / 3} + {(1 - v)}^{2 / 3}] \end{matrix}

(78)

(Note that the c₂ terms in Eq. cancel.) Using the FEM solution, we may fit the coefficients of the approximate form, yielding

\begin{matrix} - \partial θ^{e} / \partial n ≐ (0.8754) (1 / 3) [{(v + 1)}^{- 2 / 3} + {(1 - v)}^{- 2 / 3}] \\ - (0.03087) (5 / 3) [{(v + 1)}^{2 / 3} + {(1 - v)}^{2 / 3}] \end{matrix}

(79)

This approximation has a root-mean-square error of 0.0136, and is plotted as a dotted curve in Fig. 11. The result is easy to integrate

\begin{matrix} - \int_{- 1}^{+ 1} \frac{\partial θ^{e}}{\partial n} d v = (0.8754) {[{(v + 1)}^{1 / 3} - {(1 - v)}^{1 / 3}]}_{- 1}^{+ 1} \\ - (0.03087) {[{(v + 1)}^{5 / 3} - {(1 - v)}^{5 / 3}]}_{- 1}^{+ 1} \end{matrix}

(80)

\begin{matrix} = (0.8754) 2^{4 / 3} - (0.03087) 2^{8 / 3} \\ = 2.2059 - 0.1960 = 2.010 \end{matrix}

(81)

which is within 0.49% of the expected value, ${\dot{Q}}_{ext} = 2$ ⁠. The exterior normal gradient can be compared to the interior normal gradient, which is $\partial θ / \partial n = 1$ and is shown as a dashed line in Fig. 11.

The temperature on the external adiabatic side of the square,

θ_{ae}

⁠, is shown in Fig. 12. The exterior temperature is clearly different than the linear temperature distribution inside the square. A fit to

θ_{ae}

⁠, with

u \in [- 1, 1]

⁠, is

\begin{matrix} θ_{ae} ≐ \underset{c_{1}}{\underset{︸}{0.9044}} [{(u + 1)}^{1 / 3} - {(1 - u)}^{1 / 3}] + \underset{c_{2}}{\underset{︸}{0.07753}} (2 u) \\ \underset{c_{3}}{\underset{︸}{- 0.09277}} [{(u + 1)}^{5 / 3} - {(1 - u)}^{5 / 3}] \end{matrix}

(82)

Fig. 12

View large Download slide

Temperature θ along the exterior adiabatic edge of a square. Least-squares curve lies on top of FEM curve.

(Here, c₂ was constrained to make $θ_{ae} = \pm 1$ at the endpoints: $2 c_{2} = 1 - c_{1} 2^{1 / 3} - c_{3} 2^{5 / 3}$ ⁠.) This fit has a root-mean-square error of 0.0025; it is plotted as a dotted curve in Fig. 12. The value of c₁ for the temperature fit and the flux fit differ by 3.3%, likely as a result of the limited FEM resolution of the flux singularity.

We have both the interior and exterior gradient, so we are able to find the potential density,

a^{h} (v)

⁠, from its definition, Eq. , with

\partial θ / \partial n = 1

on the inside and our fit for

\partial θ^{e} / \partial n

⁠, Eq.

\begin{matrix} a^{h} (v) = \partial θ / \partial n - \partial θ^{e} / \partial n \\ = (0.8754) (1 / 3) [{(v + 1)}^{- 2 / 3} + {(1 - v)}^{- 2 / 3}] + 1 \\ - (0.03087) (5 / 3) [{(v + 1)}^{2 / 3} + {(1 - v)}^{2 / 3}] \end{matrix}

(83)

Similarly, we have b(u) along the adiabatic edges from its definition, Eq. , along with θ^e from Eq. and

θ = u

\begin{matrix} b (u) = θ^{e} - θ \\ = (0.9044) [{(u + 1)}^{1 / 3} - {(1 - u)}^{1 / 3}] \\ - (0.84494) u - (0.09277) [{(u + 1)}^{5 / 3} - {(1 - u)}^{5 / 3}] \end{matrix}

(84)

7.3 Induced Flux on a Square: Kernels and Integration.

We may now evaluate Eq. for the square. The configuration is illustrated in Fig. 13. On the hot side, parameterize the location of s with a coordinate

\hat{s}

⁠. On the cold side, parameterize

ξ

by

\hat{ξ}

⁠, and on the adiabatic sides, use

{\hat{ξ}}_{1}

and

{\hat{ξ}}_{2}

⁠. All four coordinates are

\in [- 1, 1]

⁠. Then

cos (ξ - s, n_{s}) = cos (π - α) = - cos α = - \frac{2}{| ξ - s |}

(85)

so that

\frac{\partial E (ξ | s)}{\partial n_{s}} = \frac{cos (ξ - s, n_{s})}{2 π | ξ - s |} = - \frac{1}{π | ξ - s |^{2}} = - \frac{1}{π [4 + {(\hat{s} - \hat{ξ})}^{2}]}

(86)

Fig. 13

View large Download slide

Configuration for kernels on a square

This kernel would vanish if ξ were on the hot side.

For one adiabatic side, as shown, with Eq. , letting

l = n_{ξ}

and

k = n_{s}

⁠, so that

ψ = β

and

γ = β + π / 2

\frac{\partial^{2} E (x | ξ_{2})}{\partial n_{x} \partial n_{ξ}} |_{x \to s} = \frac{\partial D_{l} (s | ξ_{2})}{\partial k} = - \frac{cos (γ + ψ)}{2 π r^{2}}

(87a)

= - \frac{cos (2 β + π / 2)}{2 π | s - ξ_{2} |^{2}} = \frac{sin (2 β)}{2 π | s - ξ_{2} |^{2}} = \frac{2 sin β cos β}{2 π | s - ξ_{2} |^{2}}

(87b)

= \frac{(1 + \hat{s}) (1 - {\hat{ξ}}_{2})}{π {[{(1 + \hat{s})}^{2} + {(1 - {\hat{ξ}}_{2})}^{2}]}^{2}}

(87c)

And, for the other adiabatic side

\frac{\partial^{2} E (x | ξ_{1})}{\partial n_{x} \partial n_{ξ}} |_{x = s} = \frac{(1 - \hat{s}) (1 - {\hat{ξ}}_{1})}{π {[{(1 - \hat{s})}^{2} + {(1 - {\hat{ξ}}_{1})}^{2}]}^{2}}

(88)

With these kernels, Eq. can be written

\begin{matrix} \int_{Γ_{s}^{h}} \int_{Γ_{ξ}^{i}} a (ξ) \frac{\partial E (ξ | s)}{\partial n_{s}} d l_{ξ} d l_{s} + \int_{Γ_{s}^{h}} \frac{\partial}{\partial n_{x}} \int_{Γ_{ξ}^{a}} b (ξ) \frac{\partial E (x | ξ)}{\partial n_{ξ}} d l_{ξ} |_{x \to s} d l_{s} \\ = - \int_{- 1}^{+ 1} \int_{- 1}^{+ 1} a^{c} (\hat{ξ}) \frac{1}{π [4 + {(\hat{s} - \hat{ξ})}^{2}]} d \hat{ξ} d \hat{s} + \int_{- 1}^{+ 1} \int_{- 1}^{+ 1} \frac{b (\hat{ξ})}{π} {\frac{(1 - \hat{s}) (1 - \hat{ξ})}{{[{(1 - \hat{s})}^{2} + {(1 - \hat{ξ})}^{2}]}^{2}} + \frac{(1 + \hat{s}) (1 - \hat{ξ})}{{[{(1 + \hat{s})}^{2} + {(1 - \hat{ξ})}^{2}]}^{2}}} d \hat{ξ} d \hat{s} \\ = - \int_{- 1}^{+ 1} \frac{a^{c} (\hat{ξ})}{2 π} [{tan}^{- 1} (\frac{1 - \hat{ξ}}{2}) + {tan}^{- 1} (\frac{1 + \hat{ξ}}{2})] d \hat{ξ} + \int_{- 1}^{+ 1} \frac{b (\hat{ξ})}{π} [\frac{1}{(1 - \hat{ξ})} - \frac{(1 - \hat{ξ})}{4 + {(1 - \hat{ξ})}^{2}}] d \hat{ξ} = 0 \end{matrix}

(89)

We may use our fits for a and b to evaluate these integrals numerically. MATLAB produces for the first integral 0.548312 and for the second integral 0.538924. The difference, 0.00938763, is 0.86% of the sum of the two values, comparable to the accuracy of the fits for a and b. We conclude that the net induced flux is zero to within the accuracy of our fitted potentials.

8 Interior and Exterior Shape Factors Are Equal

The heat conduction shape factor may be evaluated with Eqs. (4), (58), and (60).

The shape factor for the interior problem is

S^{i} = \int_{Γ_{s}^{h}} \frac{\partial θ}{\partial n_{s}} d l_{s}

(90a)

\begin{matrix} = \int_{Γ_{s}^{h}} \frac{a (s)}{2} d l_{s} + \int_{Γ_{s}^{h}} \int_{C_{ξ}} a (ξ) \frac{\partial E (ξ | s)}{\partial n_{s}} d l_{ξ} d l_{s} \\ + \int_{Γ^{h}} \lim_{x \to s^{-}} \frac{\partial}{\partial n_{x}} \int_{Γ^{a}} b (ξ) \frac{\partial E (x | ξ)}{\partial n_{ξ}} d l_{ξ} d l_{s} \end{matrix}

(90b)

= \int_{Γ_{s}^{h}} \frac{a (s)}{2} d l_{s}

(90c)

For the exterior problem, the normal vector changes sign

S^{e} = \int_{Γ_{s}^{h}} \frac{\partial θ^{e}}{\partial n_{s}^{e}} d l_{s} = - \int_{Γ_{s}^{h}} \frac{\partial θ^{e}}{\partial n_{s}} d l_{s}

(91a)

\begin{matrix} = \int_{Γ_{s}^{h}} \frac{a (s)}{2} d l_{s} - \int_{Γ_{s}^{h}} \int_{C_{ξ}} a (ξ) \frac{\partial E (ξ | s)}{\partial n_{s}} d l_{ξ} d l_{s} \\ - \int_{Γ^{h}} \lim_{x \to s^{+}} \frac{\partial}{\partial n_{x}} \int_{Γ^{a}} b (ξ) \frac{\partial E (x | ξ)}{\partial n_{ξ}} d l_{ξ} d l_{s} \end{matrix}

(91b)

= \int_{Γ_{s}^{h}} \frac{a (s)}{2} d l_{s}

(91c)

Thus

S^{e} = \int_{Γ_{s}^{h}} \frac{a (s)}{2} d l_{s} = S^{i}

(92)

9 Comments on Lienhard (2019)

In my 2019 paper [1], I showed that the present integrals Eqs. (90a) and (91a) were invariant under conformal mapping and consequently equal, as shown here using boundary integrals. The present work additionally implies that the integral of the potential density $a (ξ)$ is invariant under mapping.

The boundary integral analysis of the disk in Ref. [1] applied the Poisson integral formula with the assumption that that temperature inside and outside an adiabatic edge of the disk was the same (when comparing Eqs. (13) and (17) in Ref. [1]). The results of Sec. 5.2 in the present paper justify that assumption for the special case of a disk, while showing that the assumption cannot be used in the general case.

The Green’s function analysis in Ref. [1], for a Dirichlet boundary condition, provides the correct Green’s and influence functions for the interior and exterior of the disk. Unfortunately, the change of normal derivative between Eqs. (45) and (46) in Ref. [1] are not justified for cases other than the disk. The present work, which properly accounts for discontinuities of the normal derivative, shows that the final conclusion is nevertheless correct for arbitrary shapes (Eq. (92)). In addition, Eq. (36) in Ref. [1], while correct for the disk, does not hold in more general cases. The present work, using known kernels rather than Green’s functions, rectifies those shortcomings.

10 Summary and Conclusions

We have considered heat conduction inside and outside 2D closed curves that have two isothermal segments at different temperatures separated by two adiabatic segments. We have formulated the mixed boundary value problem using simple and double layer potentials, under the condition of zero net heat transfer to the far field. The results found are:

The outward normal temperature gradient for the region interior to an isothermal boundary is the sum of one-half the local boundary source density and integrals representing the gradient induced by sources and dipoles on other parts of the boundary (see Eq. (58)). For the region exterior to that boundary, one-half the source density is subtracted when computing that same gradient.
The induced flux has the same magnitude and orientation for the regions interior and exterior to the boundary. The flux of the local source density flows away from the boundary on each side. As a result, the combination produces unequal normal heat fluxes on the regions interior and exterior to an isothermal segment of the boundary. In addition, the temperatures are unequal for the interior and exterior of an adiabatic segment (see Eq. (59)).
In the special case of a circular disk, the induced flux is zero, with the results that: (a) the outward normal fluxes are equal and opposite on either side of an isothermal section and (b) the temperatures are equal on each side of an adiabatic segment (Sec. 4.4.2).
The integral of the induced flux over the length of either the hot or the cold boundary is zero for either the interior or the exterior region (see Eq. (60)). This behavior is analogous to the behavior of the induced charge on an ungrounded conductor in electrostatics.
Numerical solutions for the source density, $a (ξ)$ ⁠, and the dipole density, $b (ξ)$ ⁠, are given for a circular disk (Schwarz’s problem) and for a square. The numerical results accurately fit the theoretical forms for corner singularities and show an induced flux that integrates to zero (Sec. 7).
Potential theory, like conformal mapping [1], shows that the interior and exterior conduction shape factors are equal (see Eq. (92)).
An adiabatic boundary at a large but finite distance from $C$ will produce nearly the same exterior behavior as the solution for $θ \to θ^{\infty}$ as $r \to \infty$ (Appendix A).

Nomenclature

a =: simple layer (source) potential distribution, Eq. (53)
b =: double layer (dipole) potential distribution, Eq. (54)
$C$ =: boundary contour
c_m =: coefficient of corner series, Eq. (66)
$D_{l} (x | t)$ =: potential at x of dipole with orientation l at t, Eq. (15a)
dl =: differential length along contour
$E (x | t)$ =: potential at x of source at t, Eq. (6)
$e_{r}, e_{γ}$ =: unit basis vectors in polar coordinates $(r, γ)$
h(s) =: dimensionless boundary temperature
$H_{max}, N_{res}$ =: FEM parameters, see Appendix B
k =: unit vector in the k direction
$K (ψ, ψ_{0})$ =: see Eq. (64)
l =: unit vector in the l direction
m =: summation index
n =: unit normal vector
$\dot{Q}$ =: heat flow rate (W/m)
${\dot{Q}}_{ext}$ =: dimensionless heat flow rate out of square
r =: radial coordinate
R =: region inside curve
$R$ =: radius of disk
R^e =: region exterior to curve
$s, s$ =: point on contour, vector to point
S =: shape factor, Eq. (4)
t, t =: point in the plane, vector to point
T =: dimensional temperature (K)
u, v =: Cartesian coordinates, Fig. 10(a)
x, x =: point in the plane, vector to point
z =: point on $C$ in the complex plane

$α, α_{m}, β, γ, υ, ψ$ =: angles (by context)
Γ =: section of contour
$δ (x - t)$ =: Dirac delta function
ε =: a small number
η_m =: see Eq. (66)
θ =: dimensionless temperature
θ_d =: see Eq. (47a)
λ =: thermal conductivity (W/m-K)
$ξ, ξ$ =: point on $C$ ⁠, vector to point
σ =: surface in 3D
$ϕ$ =: angle in polar coordinates
$φ, Φ$ =: see Eq. (41)

a =: adiabatic
ae =: adiabatic edge
c =: cold
e =: exterior
h =: hot
i =: isothermal, or interior
$\infty$ =: at $r \to \infty$
0 =: polar coordinate for t or ξ
$\hat{}$ =: coordinate for point

FEM =: finite element method
2D =: two dimensions
3D =: three dimensions

1

Equation (25) is also referred to as a single-layer potential, and it is a special case of a Newtonian potential.

2

Note that a source on a corner point will be divided between the two sides according to the included angle, α, so that $\mp$ 1/2 becomes $- (1 - α / 2 π)$ and $+ α / 2 π$ [13].

3

As a result, we do not need to treat this integral as a principal value integral.

4

The function $a (ψ)$ is continuously differentiable for $0 < ψ < π / 2$ ⁠, and the derivative is bounded on any closed subset of that interval, e.g., $ε \leq ψ \leq π / 2 - ε$ ⁠. From the mean value theorem, $a (ψ)$ is Lipschitz continuous and, therefore, Hölder continuous on that subinterval. Near ψ = 0 and $π / 2$ ⁠, the derivative is unbounded; but we know these singularities to be integrable. We do not expect the Sokhotski–Plemelj relation to apply at ψ = 0 and $π / 2$ ⁠.

Appendix A: Far-Field Temperature, Mean Temperature on $C$ ⁠, and Adiabatic Outer Boundary

Consider first the far field (⁠

| x | = r \to \infty

⁠). Beginning with Eq.

θ^{e} (x) = θ^{\infty} + \int_{Γ^{i}} a (ξ) E (x | ξ) d l_{ξ} + \int_{Γ^{a}} b (ξ) \frac{\partial E (x | ξ)}{\partial n_{ξ}} d l_{ξ}

(A1)

we may introduce the far-field asymptotic forms of the two kernels for

r \to \infty

from Eqs. and to write

θ^{e} (x) \sim θ^{\infty} + \int_{Γ^{i}} a (ξ) [- \frac{1}{2 π} log r + \frac{1}{2 π} \frac{r_{0}}{r} cos (ϕ - ϕ_{0}) + \dots] d l_{ξ}

(A2)

+ \int_{Γ^{a}} b (ξ) [- \frac{1}{2 π r} - \frac{1}{2 π} \frac{r_{0}}{r^{2}} cos (ϕ - ϕ_{0}) - \dots] d l_{ξ}

(A3)

where

x = (r, ϕ) and ξ = (r_{0}, ϕ_{0})

⁠, and applying Eq.

= θ^{\infty} + \frac{1}{2 π r} \int_{Γ^{i}} a (ξ) r_{0} cos (ϕ - ϕ_{0}) {dl}_{ξ}

(A4)

- \frac{1}{2 π r} \int_{Γ^{a}} b (ξ) d l_{ξ} + O (\frac{1}{r^{2}})

(A5)

Because

1 / r

is not a solution of the Laplace equation in polar coordinates (cf. Eq. ), we must have

\int b (ξ) d ξ = 0

⁠. It follows that the far field of temperature is:

θ^{e} (x) = θ^{\infty} + \frac{1}{2 π r} \int_{Γ^{i}} a (ξ) r_{0} cos (ϕ - ϕ_{0}) d l_{ξ} + O (\frac{1}{r^{2}})

(A6)

The integral may be adjusted with trigonometric identities

\begin{matrix} \int_{Γ^{i}} a (ξ) r_{0} cos (ϕ - ϕ_{0}) d l_{ξ} \\ = \int_{Γ^{i}} a (ξ) r_{0} (cos ϕ cos ϕ_{0} + sin ϕ sin ϕ_{0}) d l_{ξ} \\ = cos ϕ \int_{Γ^{i}} a (ξ) r_{0} cos ϕ_{0} d l_{ξ} + sin ϕ \int_{Γ^{i}} a (ξ) r_{0} sin ϕ_{0} d l_{ξ} \\ = A cos ϕ + B sin ϕ = C cos (ϕ - α) \end{matrix}

(A7)

where A and B are defined as indicated,

C^{2} = A^{2} + B^{2}

and

tan α = B / A

⁠. Thus, the far field decays to

θ^{\infty}

plus a dipole of strength C and orientation α (cf. Eq. )

θ^{e} (x) = θ^{\infty} + \frac{C cos (ϕ - α)}{2 π r} + O (\frac{1}{r^{2}}) as r \to \infty

(A8)

Now consider the average temperature on

C

⁠, noting Eq.

\begin{matrix} \int_{C} θ^{e} (s) d l_{s} = θ^{\infty} \int_{C} d l_{s} + \int_{C_{s}} \int_{Γ_{ξ}^{i}} a (ξ) E (s | ξ) d l_{ξ} d l_{s} \\ + \int_{C_{s}} \int_{Γ_{ξ}^{a}} b (ξ) \frac{\partial E (s | ξ)}{\partial n_{ξ}} d l_{ξ} d l_{s} \end{matrix}

(A9)

The last two integrals are, in general, not equal to zero. However, in cases with sufficient symmetry, they can be zero; and then the far field temperature required for zero far-field heat transfer is simply the average temperature around the perimeter. As an example, consider the square shown in Fig. 10(a): the effect of the single integrals of a and of b on points s to the left of x = 0 is clearly opposite to their effect on points to the right of x = 0. According, the contributions of each of those integrals to the integral on s around $C$ will cancel.

And, as a counterexample, suppose that one isothermal boundary of the square is replaced by a deeply and densely V-grooved isothermal boundary at the same temperature. The corrugated boundary has a much greater length than the boundary it replaced, but much of that boundary will have low heat flux. Consequently, although the lengthy corrugated boundary will dominate the average temperature on $C$ ⁠, the total heat transfer may not be much increased. The far-field temperature for which heat flows only between the isothermal boundaries is not as strongly affected.

Finally, we note that Eq. (A8) shows that the heat flux decays as $O (1 / r^{2})$ for $r \to \infty$ ⁠. Thus, an adiabatic boundary at a finite distance from $C$ will closely approximate the result for $θ \to θ^{\infty}$ as $r \to \infty$ ⁠.

Appendix B: Finite Element Method Convergence

Convergence of the FEM model for the disk is shown in Table 2, for S (or $\partial θ / \partial n$ ⁠) approaching the theoretical value of one. Here, $H_{max}$ is the maximum mesh size of the FEM mesh. Mesh convergence to S = 1 is rapid. A radius $R_{eff} < 1$ is used for interpolation of the FEM calculation onto the curve defining the disk edge, so as to avoid NaN errors. The number of interpolated grid points, $N_{res}$ ⁠, was set to 1001 after testing grid convergence. The results in Figs. 8 and 9 are for $H_{max} = 0.001$ ⁠.

Table 2

Convergence of FEM results inside the disk

$H_{max}$	$R_{eff}$	S
0.1	0.9999	0.9069
0.05	0.9999	0.9314
0.01	0.9999	0.9692
0.005	0.9999	0.9782
0.005	0.99999	0.9788
0.001	0.9999	0.9963
0.001	0.99999	0.9986

For the square, the FEM code has some difficulty capturing the singularities at the corners. The first approach was to apply the code used for Fig. 10(b), which simulates the full perimeter of the square (see description in [1]). Table 3(a) shows the result of successive mesh refinements for that code (only $H_{max}$ shows slow grid convergence), where ${\dot{Q}}_{ext}$ is the dimensionless heat flow out of one isothermal edge of the square. The mesh used was not adaptive, and the simulations were slow to converge. (The 0.005 mesh required a 66 million node FEM computation, taking about 4 h on a laptop.)

Table 3

Convergence of the uncorrected FEM results outside a square toward the theoretical value of ${\dot{Q}}_{ext} = 2$ ⁠: (a) full-square simulation and (b) quarter-square simulation

$H_{max}$	$N_{res}$	${\dot{Q}}_{ext}$
(a)
0.4	2001	1.1346
0.2	2001	1.3052
0.1	2001	1.4446
0.05	2001	1.5588
0.025	1001	1.6512
0.025	2001	1.6513
0.0125	1001	1.7262
0.0125	8001	1.7275
0.005	2001	1.8047
$H_{max}$	$R$	${\dot{Q}}_{ext}$
(b)
0.0500	6	1.650
0.0500	8	1.591
0.0500	9	1.576
0.0500	10	1.565
0.0500	11	1.557
0.0500	12	1.551
0.0050	9	1.8252
0.0025	9	1.8712

$H_{max}$	$N_{res}$	${\dot{Q}}_{ext}$
(a)
0.4	2001	1.1346
0.2	2001	1.3052
0.1	2001	1.4446
0.05	2001	1.5588
0.025	1001	1.6512
0.025	2001	1.6513
0.0125	1001	1.7262
0.0125	8001	1.7275
0.005	2001	1.8047
$H_{max}$	$R$	${\dot{Q}}_{ext}$
(b)
0.0500	6	1.650
0.0500	8	1.591
0.0500	9	1.576
0.0500	10	1.565
0.0500	11	1.557
0.0500	12	1.551
0.0050	9	1.8252
0.0025	9	1.8712

The code can be improved by switching to a quarter-square simulation, with a circular outer boundary of radius $R$ and $N_{res} = 2001$ (Table 3(b)). The region simulated is substantially smaller, allowing a finer mesh to be used without exceeding MATLAB’s maximum matrix size of 32 GB. At large radius, the solution is essentially isothermal, and value of $R = 9$ seems adequate to maintain the outer boundary condition accurately.

The results in Figs. 11 and 12 are for $H_{max} = 0.005, N_{res} = 1001$ ⁠, and $R = 9$ on a half-square model.

The gradient computed by the FEM code at the last few grid points near the corner goes to an unphysical value of zero: the actual gradient vector is undefined at the corner, but its magnitude tends to infinity as the corner is approached along an isothermal edge. A simple correction to the integral for ${\dot{Q}}_{ext}$ ⁠, addressing the lost area near the corner, can be made as follows. We integrate the fit Eq. (79) from –1 to a position slightly away from the singularity, $- 1 + δ y$ ⁠. That area is $\approx (0.8754) {(δ y)}^{1 / 3}$ ⁠. Then, the numerical integration of the FEM result may be confined to $| y | \leq 1 - δ y$ ⁠; and when $2 (0.8754) {(δ y)}^{1 / 3}$ is added to the numerical integral, we obtain ${\dot{Q}}_{ext} = 2.0202$ ⁠, with the following matlab inputs: $H_{max} = 0.005, N_{res} = 1001, δ y = 0.01$ ⁠, and $R$ = 9. The uncorrected value is ${\dot{Q}}_{ext} = 1.8213$ ⁠, and the theoretically known value is ${\dot{Q}}_{ext} = S Δ θ = 2$ ⁠. The distance δy must be large enough to avoid the unphysical gradient, 0, that the FEM code returns at $| y | = 1$ ⁠.

References

1.

Lienhard

,

J. H.

, V ,

2019

, “

Exterior Shape Factors From Interior Shape Factors

,”

ASME J. Heat Transfer-Trans. ASME

,

141

(

6

), p.

061301

. 10.1115/1.4042912

Google Scholar

Crossref

2.

Stakgold

,

I.

,

1967

,

Boundary Value Problems of Mathematical Physics

, Vol.

2

,

Macmillan

,

New York

.

3.

Nehari

,

Z.

,

1952

,

Conformal Mapping

,

McGraw-Hill Book Company

,

New York

.

4.

Kellogg

,

O. D.

,

1929

,

Foundations of Potential Theory

,

Verlag Von Julius Springer

,

Berlin, Germany

.

Google Scholar

Crossref

5.

MacMillan

,

W. D.

,

1930

,

The Theory of the Potential

,

McGraw-Hill Book Company

,

New York

.

6.

Morse

,

P. M.

, and

Feshbach

,

H.

,

1953

,

Methods of Theoretical Physics

,

McGraw-Hill Book Company

,

New York

.

7.

Martinsson

,

G.

,

2014

, CBMS Conference on Fast Direct Solvers, Dartmouth College, Lectures 1 and 8, Hanover, NH, accessed Oct. 17, 2022, https://math.dartmouth.edu/~fastdirect/

8.

Atkinson

,

K. E.

,

1997

,

Boundary Integral Equations on a Piecewise Smooth Planar Boundary

(Chap. 8, Cambridge Monographs on Applied and Computational Mathematics),

Cambridge University Press

,

Cambridge

, pp.

384

–

426

.

9.

Atkinson

,

K. E.

,

2012

, “Boundary Integral Codes for the Planar Laplace Equation on Interior or Exterior Regions Having Dirichlet or Neumann Boundary Conditions,” MATLAB Codes for Planar Regions, Software package,

University of Iowa

,

Ames, IA

, accessed Oct. 24, 2022, https://homepage.divms.uiowa.edu/~atkinson/laplace.html

10.

Green

,

G.

,

1828

,

An Essay on the Application of Mathematical Analysis to the Theories of Electricity and Magnetism

,

Author, Nottingham

,

UK

.

11.

Poincaré

,

H.

,

Le Roy

,

É.

, and

Vincent

,

G.

,

1899

,

Théorie du Potentiel Newtonien: Leçons Professés à la Sorbonne Pendant le Premier Semestre 1894-1895

,

Georges Carré et C. Naud, Paris

,

France

.

12.

Xu

,

Y.

,

Chen

,

H.-L.

, and

Zou

,

Q.

,

2004

, “

Limit Values of Derivatives of the Cauchy Integrals and Computation of the Logarithmic Potentials

,”

Computing

,

73

(

4

), pp.

295

–

327

.10.1007/s00607-003-0072-4

Google Scholar

Crossref

13.

Muskhelishvili

,

N. I.

,

1946

,

Singular Integral Equations: Boundary Problems of Function Theory and Their Application to Mathematical Physics

,

P. Noordhoff N.V., Groningen

,

Holland

.

14.

Jackson

,

J. D.

,

1999

,

Classical Electrodynamics

, 3rd ed.,

Wiley

,

New York

.

15.

Günter

,

N. M.

,

1967

,

Potential Theory and Its Applications to Basic Problems of Mathematical Physics

,

Frederick Ungar Publishing Co

,

New York

(Originally published in French in 1934. This edition is an English translation of a later Russian edition).

16.

Schmidt

,

E.

,

1914

, “

Bemerkung Zur Potentialtheorie

,”

Mathematische Abhandlungen H. A. Schwarz Gewidmet

,

J. Springer

,

Berlin

, pp.

365

–

383

(Mathematical Treatises dedicated to Hermann Amandus Schwarz by friends and students on August 6, 1914 to mark the 50th anniversary of his doctorate).

Google Scholar

Crossref

17.

Hinrichsen

,

J. J. L.

,

1935

, “

Note on Potential Theory in n-Space

,”

Am. J. Math.

,

57

(

4

), pp.

921

–

927

.10.2307/2371029

Google Scholar

Crossref

18.

Schwarz

,

H. A.

,

1869

, “

Über Einige Abbildungsaufgaben

,”

J. Die Reine Angew. Math.

,

70

, pp.

105

–

120

.10.1016/j.ymssp.2016.09.007

19.

Igarashi

,

H.

, and

Honma

,

T.

,

1996

, “

A Boundary Element Method for Potential Fields With Corner Singularities

,”

Appl. Math. Modell.

,

20

(

11

), pp.

847

–

852

.10.1016/S0307-904X(96)00091-1

Google Scholar

Crossref

20.

Atkinson

,

K.

, and

De Hoog

,

F.

,

1984

, “

The Numerical Solution of Laplace’s Equation on a Wedge

,”

IMA J. Numer. Anal.

,

4

(

1

), pp.

19

–

41

.10.1093/imanum/4.1.19

Google Scholar

Crossref

21.

Kobayashi

,

S.

,

Nishimura

,

N.

, and

Kawakami

,

T.

,

1984

, “

Simple Layer Potential Method for Domains Having External Corners

,”

Appl. Math. Modell.

,

8

(

1

), pp.

61

–

65

.10.1016/0307-904X(84)90180-X

Google Scholar

Crossref

Meaning	Primary symbol	Full expression	Equivalent symbol
Potential at s caused by unit source density at ξ	$E (s \| ξ)$	$\frac{1}{2 π} log \frac{1}{\| s - ξ \|}$	$E (ξ \| s)$
Outward normal gradient at s caused by unit source density at ξ	$\frac{\partial E (ξ \| s)}{\partial n_{s}}$	$\frac{cos (ξ - s, n_{s})}{2 π \| ξ - s \|}$	$D_{n_{s}} (ξ \| s)$
Potential at s caused by unit outward-normal dipole density at ξ	$\frac{\partial E (s \| ξ)}{\partial n_{ξ}}$	$\frac{cos (s - ξ, n_{ξ})}{2 π \| s - ξ \|}$	$D_{n_{ξ}} (s \| ξ)$

Meaning	Primary symbol	Full expression	Equivalent symbol
Potential at s caused by unit source density at ξ	$E (s \| ξ)$	$\frac{1}{2 π} log \frac{1}{\| s - ξ \|}$	$E (ξ \| s)$
Outward normal gradient at s caused by unit source density at ξ	$\frac{\partial E (ξ \| s)}{\partial n_{s}}$	$\frac{cos (ξ - s, n_{s})}{2 π \| ξ - s \|}$	$D_{n_{s}} (ξ \| s)$
Potential at s caused by unit outward-normal dipole density at ξ	$\frac{\partial E (s \| ξ)}{\partial n_{ξ}}$	$\frac{cos (s - ξ, n_{ξ})}{2 π \| s - ξ \|}$	$D_{n_{ξ}} (s \| ξ)$

Steady Two-Dimensional Conduction: Simple and Double Layer Potentials, Corner Singularities, and Induced Heat Flux

Abstract

1 Introduction

2 Preliminaries: Sources, Dipoles, and Gradients

2.1 Fundamental Solution: E(x|t)⁠.

2.2 Dipoles.

2.3 Dipoles as Induced Gradients in Two Dimensions.

2.4 Gradient of a Dipole Potential in Two Dimensions.

2.5 Some Properties of Harmonic Functions.

3 Solution for the Interior and Exterior Dirichlet Problems by a Simple-Layer Potential

3.1 Modifications of Equations (25) and (26) for Two Dimensions.

3.2 Two-Dimensional Contours With Zero Net Heat Transfer.

3.3 Example: Simple-Layer Potential on a Unit Disk.

4 Heat Flux at the Boundaries

4.1 The Simple-Layer Normal Derivative as a Cauchy Integral.

4.2 Sokhotski–Plemelj Theorem.

4.3 Plemelj Formulæ for the Normal Derivative.

4.3.1 More About Induced Heat Flux.

4.4 Special Cases of the Normal Derivative.

4.4.1 Straight Lines and Parallel Plates.

4.4.2 Circular Disk.

4.5 Plemelj Formulæ for the Integral of the Normal Derivative (Double-Layer Potential).

5 Mixed Boundary Conditions as Simple and Double Layer Potentials

5.1 Induced Heat Flux.

5.2 Mixed Boundary Conditions on a Disk.

5.2.1 Fredholm Equation for Schwarz’s Problem.

6 Singularities at Corners and Dirichlet-to-Neumann Transitions

6.1 Corner Singularities.

6.2 Change From Dirichlet to Neumann Boundary Condition on a Straight Section.

7 Finite Element Method Studies of Potential and Induced Flux

7.1 Disk With Isothermal and Adiabatic Edges.

7.2 Square With Isothermal and Adiabatic Sides.

7.2.1 Fitting Heat Flux and Temperature With Corner Theory.

7.3 Induced Flux on a Square: Kernels and Integration.

8 Interior and Exterior Shape Factors Are Equal

9 Comments on Lienhard (2019)

10 Summary and Conclusions

Nomenclature

Appendix A: Far-Field Temperature, Mean Temperature on C⁠, and Adiabatic Outer Boundary

Appendix B: Finite Element Method Convergence

References

Get Email Alerts

Cited By

Related Articles

Related Proceedings Papers

Related Chapters

ASME Journals

ASME Conference Proceedings

ASME eBooks

Resources

Opportunities

This Feature Is Available To Subscribers Only

2.1 Fundamental Solution: $E (x | t)$ ⁠.

Appendix A: Far-Field Temperature, Mean Temperature on $C$ ⁠, and Adiabatic Outer Boundary