## 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.

*R*bounded by a simple closed curve, $C$, that contains two isothermal segments, $\Gamma h$ and $\Gamma c$, separated by two adiabatic segments, $\Gamma 1a$ and $\Gamma 2a$. Let

*s*be a point on $C$. Steady heat conduction in

*R*is governed by the Laplace equation for temperature,

*θ*, scaled to $\u22121\u2264\theta \u22641$

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

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, $ne$, 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, $\Gamma h$, to the cold segment of $\Gamma c$. In this case, the far field temperature $\theta \u221e$ 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].

*S*, in 2D provides the heat transfer rate, $Q\u02d9$ (W/m), between the hot and cold boundaries in terms of the dimensional temperature difference, $\Delta T=Th\u2212Tc$, and thermal conductivity,

*λ*, as $Q\u02d9=S\lambda \Delta T$. For either the interior or the exterior of $C$, the shape factor may be written as an integral on $\Gamma h$ of the appropriate outward normal derivative of

*θ*[1]

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)$.

*x*caused by a unit point source at a point

*t*[2]

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

*R*using the divergence theorem

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 $\u2212\u2207x2T=(P/\lambda )\delta (x\u2212t)$ and the leaving heat flux is $\u2212\lambda \u2202T/\u2202nx$. The net heat flow out of *R* through $C$ is *P* if $t\u2208R$ 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.

*multipole expansion*)

where $r>=max(r,r0)$ and $r<=min(r,r0)$.

### 2.2 Dipoles.

*ε*, the result of letting $\epsilon \u21920$ 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+\epsilon l$ approaching a sink at

*t*(Fig. 3). Call the potential of this dipole $Dl(x|t)$. In 2D

where *γ* is the angle between the vectors $x\u2212t$ and **l**. Note that $Dl(x|t)=\u2212Dl(t|x)$.

### 2.3 Dipoles as Induced Gradients in Two Dimensions.

*t*produced by a source at

*x*is simply $\u2207tE(t|x)=\u2207tE(x|t)$, and this gradient is collinear with $x\u2212t$ (pointing from

*t*to

*x*). The derivative of $E(x|t)$ at

*t*in the direction

**l**is

Thus, we may consider $Dl(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*, $nt$. If both *x* and *t* are on this surface and $l=nt,\u2009\u2009cos\u2009\gamma \u21920$ as $|x\u2212t|\u21920$. Analysis shows that $Dnt(x|t)$ is continuous at *x *=* t* when *x* approaches *t* along the surface (as opposed to approaching from another direction).

*x*produced by a source at

*t*is

where $cos(t\u2212x,nx)$ denotes the cosine of the angle between vectors $t\u2212x$ and $nx$. This expression may be integrated over the surface that contains the source distribution, *a*(*t*), holding the vector $nx$ 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.

**l**at

*t*and that we wish to know the potential’s gradient in the direction

**k**at point

*x*(Fig. 4). We find

*γ*is the angle between the vectors $x\u2212t$ and

**l**. For convenience, we may use polar coordinates $(r,\gamma )$ with origin at

*t*and $r=|x\u2212t|$. Then

**k**makes an included angle

*ψ*with respect to $x\u2212t$: $k=cos\u2009\psi \u2009er\u2212sin\u2009\psi \u2009e\gamma $. We find

### 2.5 Some Properties of Harmonic Functions.

*r*>

*0, the separation of variables in polar coordinates gives*

To have a bounded solution as $r\u2192\u221e$, we require *B *=* *0 and $cm+=0$. As consequences, $\theta \u2192A=\theta \u221e$ as $r\u2192\u221e$; and for every circle $r=constant$, the average of *θ* around the circumference is $\theta \u221e$.

This gradient $\u21920$ as $r\u2192\u221e$.

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 $\theta \u221e$, and has zero average heat flux through it.

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

*R*) and outside (

*R*) a 3D surface

^{e}*σ*may be written as a simple-layer potential [10]

^{1}

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 function $a(\xi )$ may be considered to be a heat source distribution on *σ*, with heat flowing out to each side.

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

*R*, and suppose it to be bounded by the contour $C$ and a very distant circular contour of radius

^{e}*r*

_{0}, on which $\xi =(r0,\varphi 0)$. We treat $E(x|\xi )$ and $\theta e(\xi )$ as functions of

*ξ*and apply Green’s second identity to these functions

*ξ*indicates that differentiation or integration is with respect to the

*ξ*coordinate. The left-hand side simplifies upon substituting Eqs. (1) and (5)

These results show that $\theta \u221e$ is imposed by the far field in 2D.

so that $\theta e(x)\u2192\theta \u221e$ as expected.

### 3.2 Two-Dimensional Contours With Zero Net Heat Transfer.

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

*θ*and

*θ*just as would the separation of variables solution, Eq. (23)

^{e}## 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.

### 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\u2208C$, and consider the fixed point $x\u2209C$, near *s*. We can differentiate at under the integral with respect to *x* in the direction $ns$. The angle between $ns$ and the real axis is $\varphi 0$ (Fig. 5). Here, *x*, *s*, and *z* are points in the complex plane.

where the angle for the cosine in Eq. (38*a*) is $\varphi 0\u2212arg(z\u2212x)$.

*d*(

*x*) defined as

where *a*(*z*) is real valued.

### 4.2 Sokhotski–Plemelj Theorem.

where $\phi (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\u2209C$; further, $\Phi =O(|x|\u22121)$ as $|x|\u2192\u221e$. When $x=s\u2208C$, the integral exists only in the Cauchy principal value sense.

*s*on $C$ to the principal value of the integral (/$\u222b$)

*s*

^{+}indicates approaches from outside the contour and

*s*

^{–}indicates approach from inside. The analysis requires that $\phi $ be Hölder continuous

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

### 4.3 Plemelj Formulæ for the Normal Derivative.

**s**on $C$, is

As noted in Sec. 2.3, this kernel does not have a singularity at $\xi =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. (45*b*), includes a source term and an integral representing gradients in the normal direction at $s$ induced by sources at other points $\xi $ 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.^{2} The induced gradient term depends on the distance to the source point, $|\xi \u2212s|$, and the orientation of $\xi \u2212s$ relative to the fixed point’s outward normal vector, $ns$. 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\u2192s\xb1$. 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(\xi )=0$ (the boundary contains no heat source distribution), and so Eq. (45*b*) 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(\xi )=0$ on the disk.

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. (45*b*) 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(\xi \u2212s,n)=0$ if $\xi $ 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(\xi \u2212s,ns)/|\xi \u2212s|=\u22121/2R$, and with Eq. (35), the integral in Eq. (45*b*) will vanish at every *s*, leaving $\u2202\theta /\u2202ns=\u2212\u2202\theta e/\u2202ns$.

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).

*z*(Fig. 5). We proceed as in the simple-layer analysis

*b*(

*z*), $\u2212ei\pi /2=1/i$, and

*g*(

*x*) defined as

*s*, and with

*b*(

*s*) real valued, Eq. (42) gives

for $n\xi $ the outward normal at the moving point *ξ*.

The double-layer potential *θ _{d}* is discontinuous at $C$, but the normal derivative $\u2202\theta d/\u2202nx$ 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.

## 5 Mixed Boundary Conditions as Simple and Double Layer Potentials

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

*b*), (16), and (17)

The double layer potential has the same normal derivative for either $x\u2192s\xb1$ (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 $\Gamma i$ has a logarithmic singularity, while the second integral has a nonsingular kernel (Sec. 2.3).^{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 normal derivative of the second integral has no singularity for $s\u2208\Gamma i$, in contrast to $s\u2208\Gamma 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).

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.

Meaning | Primary symbol | Full expression | Equivalent symbol |
---|---|---|---|

Potential at s caused by unit source density at ξ | $E(s|\xi )$ | $12\pi log\u20091|s\u2212\xi |$ | $E(\xi |s)$ |

Outward normal gradient at s caused by unit source density at ξ | $\u2202E(\xi |s)\u2202ns$ | $\u2009cos(\xi \u2212s,ns)2\pi |\xi \u2212s|$ | $Dns(\xi |s)$ |

Potential at s caused by unit outward-normal dipole density at ξ | $\u2202E(s|\xi )\u2202n\xi $ | $\u2009cos(s\u2212\xi ,n\xi )2\pi |s\u2212\xi |$ | $Dn\xi (s|\xi )$ |

Meaning | Primary symbol | Full expression | Equivalent symbol |
---|---|---|---|

Potential at s caused by unit source density at ξ | $E(s|\xi )$ | $12\pi log\u20091|s\u2212\xi |$ | $E(\xi |s)$ |

Outward normal gradient at s caused by unit source density at ξ | $\u2202E(\xi |s)\u2202ns$ | $\u2009cos(\xi \u2212s,ns)2\pi |\xi \u2212s|$ | $Dns(\xi |s)$ |

Potential at s caused by unit outward-normal dipole density at ξ | $\u2202E(s|\xi )\u2202n\xi $ | $\u2009cos(s\u2212\xi ,n\xi )2\pi |s\u2212\xi |$ | $Dn\xi (s|\xi )$ |

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

### 5.1 Induced Heat Flux.

*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

The limit should be continuous, as the integral is nonsingular for $x\u2192s\u2208\Gamma h$ and $\xi \u2208\Gamma a$, except at the point where $\Gamma h$ and $\Gamma a$ intersect. The same conclusion applies to the net induced flux on the cold isothermal boundary, $\Gamma c$.

As noted in Sec. 4.5, the normal derivative of second integral has the same limit for $x\u2192s\xb1$, but is undefined for $x=s\u2208\Gamma 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 $\theta (s)=1$, the boundary in the third quadrant is isothermal at $\theta (s)=\u22121$, and the other two boundaries are adiabatic. This problem was solved by Schwarz in 1869 using conformal mapping [18].

For a circular disk, as in Sec. 4.4.2, the kernel of the integral on $a(\xi )$ in Eq. (57) reduces to a constant and the integral vanishes. The integral in $b(\xi )$ retains an *x* dependence, but the heat flux boundary condition on the adiabatic sides is satisfied by $b(\xi )=0$, and Eq. (59) then indicates equal interior and exterior temperatures on the adiabatic edges. Further, with $b(\xi )=0$, Eq. (56) reduces to a Fredholm problem for the potential on the isothermal sections of the boundary, with $a(\xi )=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 $\u2213a(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.

*ψ*is the angle around the perimeter of the a disk measured from the

*x*-axis. We know that $b(\psi )=0$ and that the heat flux boundary conditions on the adiabatic edges are met. So, Eq. (59) reduces to the following equation on the isothermal edges of the unit disk:

Note that the kernel is weakly singular at $\psi =\psi 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.

*β*, and we have one edge that is isothermal and one that is adiabatic, the solutions of the Laplace equation, from Eq. (23), will be [19]

*θ*normal to a line of constant $\varphi $ is

When $\eta 1<1$ ($\beta >\pi /2$), the derivative contains an integrable singularity as $r\u21920$. For smaller *β*, the derivative is nonsingular.

where the final equality ($c1=1$, *c _{m}* = 0 for

*m*>

*1) is from my knowledge of the solution: $\theta (r,\varphi )=1\u2212r\u2009cos\u2009\varphi $ (note that*

*n*is outward normal here).

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

*r*=

*0 (e.g., locally on the rim of a disk), $\beta =\pi $ and $\eta m=(m\u22121/2)$, so that*

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 ($\varphi =0$) the temperature distribution inside and outside will be the same *locally*.

## 7 Finite Element Method Studies of Potential and Induced Flux

### 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.

*b*), is $O(r\u22121/2)$; and if the angle around the disk is

*ψ*measured from the

*x*-axis, we may set $r=R\psi $ and $r=R(\pi /2\u2212\psi )$ 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. (73

*b*), using MATLAB’s fitting algorithms. A fit that accounts for the singularity at

*ψ*= 0 and $\psi =\pi /2$, each to three terms, is

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.

For a shape factor $Si=1$ [1] and $\Delta \theta =2$, the theoretical average gradient is $2(2/\pi )=4/\pi =1.2732$, indicated by a horizontal dotted line in Fig. 8. The integrated result is lower by only 0.11%.

The root-mean-square error of the fit is $6.4266\xd710\u22125$. 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 $c1\u22501.070,\u2009c2\u2250\u22120.134$, and $c3\u2250\u22120.025$.

The single-layer potential is defined as $a=\u2202\theta /\u2202n\u2212\u2202\theta e/\u2202n$ (Eq. (53)). For the disk, from Sec. 5.2, we know that $\u2202\theta e/\u2202n=\u2212\u2202\theta /\u2202n$, and so the single-layer potential $a(\psi )=2\u2202\theta /\u2202n$.^{4} Numerical evaluation of the integral in Eq. (65), in MATLAB, using the $a(\psi )$ 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 $\theta =\u22121$, 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.

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 $\Delta \theta =2$. The shape factor for the exterior of a square is $Se=1$ [1]; and, with the thermal conductivity normalized to unity, the dimensionless heat transfer rate is $Q\u02d9ext=2$.

#### 7.2.1 Fitting Heat Flux and Temperature With Corner Theory.

**n**. Corner theory, Eq. (70

*b*), superposed for the two corners, suggests fitting with this approximate form, where $v\u2208[\u22121,1]$

*c*

_{2}terms in Eq. (70

*b*) cancel.) Using the FEM solution, we may fit the coefficients of the approximate form, yielding

which is within 0.49% of the expected value, $Q\u02d9ext=2$. The exterior normal gradient can be compared to the interior normal gradient, which is $\u2202\theta /\u2202n=1$ and is shown as a dashed line in Fig. 11.

(Here, *c*_{2} was constrained to make $\theta ae=\xb11$ at the endpoints: $2c2=1\u2212c121/3\u2212c325/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*_{1} 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.

### 7.3 Induced Flux on a Square: Kernels and Integration.

**s**with a coordinate $s\u0302$. On the cold side, parameterize $\xi $ by $\xi \u0302$, and on the adiabatic sides, use $\xi \u03021$ and $\xi \u03022$. All four coordinates are $\u2208[\u22121,1]$. Then

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

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

## 9 Comments on Lienhard (2019)

In my 2019 paper [1], I showed that the present integrals Eqs. (90*a*) and (91*a*) 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(\xi )$ 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(\xi )$, and the dipole density, $b(\xi )$, 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 $\theta \u2192\theta \u221e$ as $r\u2192\u221e$ (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)

- $Dl(x|t)$ =
potential at

*x*of dipole with orientation*l*at*t*, Eq. (15*a*)*dl*=differential length along contour

- $E(x|t)$ =
potential at

*x*of source at*t*, Eq. (6)- $er,e\gamma $ =
unit basis vectors in polar coordinates $(r,\gamma )$

*h*(*s*) =dimensionless boundary temperature

- $Hmax,Nres$ =
FEM parameters, see Appendix B

**k**=unit vector in the

*k*direction- $K(\psi ,\psi 0)$ =
see Eq. (64)

**l**=unit vector in the

*l*direction*m*=summation index

**n**=unit normal vector

- $Q\u02d9$ =
heat flow rate (W/m)

- $Q\u02d9ext$ =
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

- $\alpha ,\alpha m,\beta ,\gamma ,\upsilon ,\psi $ =
angles (by context)

- Γ =
section of contour

- $\delta (x\u2212t)$ =
Dirac delta function

*ε*=a small number

*η*=_{m}see Eq. (66)

*θ*=dimensionless temperature

*θ*=_{d}see Eq. (47

*a*)*λ*=thermal conductivity (W/m-K)

- $\xi ,\xi $ =
point on $C$, vector to point

*σ*=surface in 3D

- $\varphi $ =
angle in polar coordinates

- $\phi ,\u2009\Phi $ =
see Eq. (41)

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

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

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

The function $a(\psi )$ is continuously differentiable for $0<\psi <\pi /2$, and the derivative is bounded on any closed subset of that interval, e.g., $\epsilon \u2264\psi \u2264\pi /2\u2212\epsilon $. From the mean value theorem, $a(\psi )$ is Lipschitz continuous and, therefore, Hölder continuous on that subinterval. Near *ψ* = 0 and $\pi /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 $\pi /2$.

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

*A*and

*B*are defined as indicated, $C2=A2+B2$ and $tan\u2009\alpha =B/A$. Thus, the far field decays to $\theta \u221e$ plus a dipole of strength

*C*and orientation

*α*(cf. Eq. (15

*a*))

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/r2)$ for $r\u2192\u221e$. Thus, an adiabatic boundary at a finite distance from $C$ will closely approximate the result for $\theta \u2192\theta \u221e$ as $r\u2192\u221e$.

### Appendix B: Finite Element Method Convergence

Convergence of the FEM model for the disk is shown in Table 2, for *S* (or $\u2202\theta /\u2202n$) approaching the theoretical value of one. Here, $Hmax$ is the maximum mesh size of the FEM mesh. Mesh convergence to *S *=* *1 is rapid. A radius $Reff<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, $Nres$, was set to 1001 after testing grid convergence. The results in Figs. 8 and 9 are for $Hmax=0.001$.

$Hmax$ | $Reff$ | 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 |

$Hmax$ | $Reff$ | 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 $Hmax$ shows slow grid convergence), where $Q\u02d9ext$ 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.)

$Hmax$ | $Nres$ | $Q\u02d9ext$ |
---|---|---|

(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 |

$Hmax$ | $R$ | $Q\u02d9ext$ |

(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 |

$Hmax$ | $Nres$ | $Q\u02d9ext$ |
---|---|---|

(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 |

$Hmax$ | $R$ | $Q\u02d9ext$ |

(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 $Nres=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 $Hmax=0.005,\u2009Nres=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 $Q\u02d9ext$, 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, $\u22121+\delta y$. That area is $\u2248(0.8754)(\delta y)1/3$. Then, the numerical integration of the FEM result may be confined to $|y|\u22641\u2212\delta y$; and when $2(0.8754)(\delta y)1/3$ is added to the numerical integral, we obtain $Q\u02d9ext=2.0202$, with the following matlab inputs: $Hmax=0.005,\u2009Nres=1001,\u2009\delta y=0.01$, and $R$* *=* *9. The uncorrected value is $Q\u02d9ext=1.8213$, and the theoretically known value is $Q\u02d9ext=S\Delta \theta =2$. The distance *δy* must be large enough to avoid the unphysical gradient, 0, that the FEM code returns at $|y|=1$.

## References

*n*-Space