1. Introduction
When a passive scalar field
$c(\boldsymbol{x},t)$
, such as temperature or a solute concentration, is transported within a fluid flow
$\boldsymbol{u}(\boldsymbol{x},t)$
, it is advected and diffused so that its evolution is governed by
\begin{equation} \frac {\partial c}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }c = \alpha {\nabla} ^2 c, \end{equation}
where
$\alpha$
is the diffusivity of
$c$
. Equation (1.1) must be supplemented by suitable boundary conditions that depend on the particular problem under consideration. If the system has developed into a steady-state, (1.1) simplifies to
\begin{equation} \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }c = \alpha {\nabla} ^2 c. \end{equation}
This is simply a linear elliptic equation for
$c$
, but due to the non-constant nature of the velocity field, obtaining a numerical solution generally requires discretising the entire fluid. To produce solutions with high accuracy, fine meshes are needed and the associated linear algebra problem becomes costly. Standard boundary-integral techniques are also difficult to apply directly to (1.2), since closed forms of the Green’s function are not available for generic velocity fields.
In the present paper, we restrict our attention to the case where
$\boldsymbol{u}$
is a steady two-dimensional (2-D) incompressible potential flow with potential
$\phi$
, so that
$\boldsymbol{u}=\boldsymbol{\nabla }\phi (x,y)$
and
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u} = 0$
, past a collection of stationary impermeable obstacles, where each obstacle maintains an arbitrary time-independent boundary temperature profile. Potential flows serve as accurate models of flow in Hele-Shaw geometries, flows through porous media and high-Reynolds-number flows in situations where viscous effects are confined to boundary layers attached to obstacle boundaries. Such flows may be represented in the complex plane in terms of a complex potential
$W(z)=\phi (z)+\mathrm{i}\psi (z)$
, where
$W(z)$
is a complex analytic function of
$z=x+\mathrm{i}y$
in the fluid and
$\psi (z)$
is the harmonic conjugate of
$\phi (z)$
. The complex flow velocity may then be expressed as
$u=\overline {{\rm d}W/{\rm d}z}$
.
In addition to describing the flow velocity, the complex potential
$W(z)$
can be interpreted as a multiply connected conformal map from the physical domain (
$z$
-domain) to the mapped domain (
$w$
-domain), also referred to as streamline coordinates. Under this transformation, both the governing equation (1.2) and the problem geometry are simplified greatly. Each impermeable obstacle is transformed into a horizontal slit, since the streamfunction
$\psi$
must be constant on the boundary of such an obstacle, and the fluid region is transformed to the entire complex plane exterior to these slits. This transformation dates back to Boussinesq (Reference Boussinesq1902). By the fortuitous cancellation of Jacobian terms, the governing equation (1.2) is transformed into a constant-coefficient equation in the mapped domain, as we will detail in § 2.2. This allows the partial differential equation (PDE) to be efficiently solved in the mapped domain either analytically (in some cases) or with boundary-integral techniques.
Historically, analytical solutions for advection–diffusion have been sought using this approach for potential flows past canonical isolated obstacles such as flat plates, circles and elliptic cylinders (Boussinesq (Reference Boussinesq1905), Hsu (Reference Hsu1965)); see Galante & Churchill (Reference Galante and Churchill1990) for a general review. Choi et al. (Reference Choi, Margetis, Squires and Bazant2005) studied advection–diffusion past an isothermal simply connected obstacle of arbitrary shape within a uniform potential flow. Using a single-layer boundary integral formulation in the mapped domain, they derived asymptotic expansions in the Péclet number for the far-field temperature profile and net heat flux. For example, for potential flow past a disk of radius
$R$
, they found that as the Péclet number increases, the Nusselt number scales asymptotically as
\begin{equation} {\textit{Nu}} \underset {{\textit{Pe}} \rightarrow \infty }{\sim } 8 \sqrt {\frac {{\textit{Pe}}}{\pi }}. \end{equation}
The Nusselt number is defined as the ratio of the total heat flux out of the obstacle to the quantity
$\alpha \Delta T$
, where
$\Delta T$
is the obstacle temperature relative to infinity. The Péclet number is given by
${\textit{Pe}} = U R / \alpha$
, where
$U$
is the far-field (free stream) velocity. In the limit of small Péclet number, the study found that
\begin{equation} {\textit{Nu}} \underset {{\textit{Pe}} \rightarrow 0}{\sim } -\frac {2\pi }{\gamma +\log {\left (\dfrac {{\textit{Pe}}}{4}\right )}}, \end{equation}
where
$\gamma = 0.57721566 \ldots$
is the Euler–Mascheroni constant. Related investigations have studied advection-diffusion limited aggregation (ADLA) (Bazant et al. Reference Bazant, Choi and Davidovitch2003; Bazant & Moffatt Reference Bazant and Moffatt2005; Davidovitch, Choi & Bazant Reference Davidovitch, Choi and Bazant2005; Rycroft & Bazant Reference Rycroft and Bazant2016) and artificial freezing (Alimov, Kornev & Mukhamadullina Reference Alimov, Kornev and Mukhamadullina1994; Mukhamadullina, Kornev & Alimov Reference Mukhamadullina, Kornev and Alimov1998) around an isothermal, simply connected obstacle. More recently, Goyette et al. (Reference Goyette, Boulais, Normandeau, Laberge, Juncker and Gervais2019) studied advection–diffusion in a flow consisting of source/sink singularities, in the absence of a free stream or any flow obstacles.
In the present paper, we consider the more general problem of advection–diffusion in multiply connected domains. Importantly, we allow the boundary of each obstacle to possess a spatially varying temperature distribution, rather than limiting ourselves to the isothermal case. First, we develop an accurate numerical procedure for computing the potential flow around a generic collection of smooth objects. This method is enabled by the AAA-LS algorithm (Costa & Trefethen Reference Costa and Trefethen2023), a flexible tool for solving Laplace problems using rational function approximation. In the mapped domain, we then solve for the temperature field using a mixed single-plus-double-layer boundary integral method that exploits the homogeneous Green’s function of the transformed PDE.
Using these tools, we investigate both similarities and key differences between the behaviour of our solutions and the isothermal, simply connected problem that was examined by Choi et al. (Reference Choi, Margetis, Squires and Bazant2005). We first compare the behaviour of the simply connected variable-temperature problem with the results of Choi et al. (Reference Choi, Margetis, Squires and Bazant2005) in the high- and low-Péclet-number limits. Notably, we demonstrate new scalings for the Nusselt number, both numerically and analytically, in the low-Péclet-number limit, based on the structure of the imposed temperature boundary conditions. We then investigate the high-Péclet-number limit and multiply connected problems.
The methods developed herein, being highly general and accurate to many digits of precision, should prove useful for researchers looking to solve advection–diffusion problems in multiply connected domains. Prior work considering non-trivial multiply connected geometries was highly limited by the lack of flexible and efficient numerical algorithms for these cases. Our new algorithms for multiply connected conformal maps using AAA-LS, and our efficient boundary integral solvers, are capable of solving many such problems with ease. Future applications of this line of work might include investigations of multiply connected ADLA and artificial growth problems, as well as extensions of these techniques to quasi-static advection–diffusion processes.
The remainder of this paper is organised as follows. In § 2, we mathematically describe the multiply connected conformal map corresponding to the potential flow solution and demonstrate how the governing PDE is simplified under this transformation. In § 3, we formulate our conformal mapping computations as well as the boundary integral technique used to obtain the temperature field in the mapped domain. We also derive and verify expressions for the flux from each obstacle. In § 4, we compare our variable temperature results, in both simply and multiply-connected domains, with existing asymptotic results pertaining to the special case of simply connected isothermal obstacles; in doing so, we demonstrate similarities and key differences in our generalised setting. In § 4.3, we derive asymptotic scalings for the Nusselt number, under variable temperature data, in the low-Péclet-number limit. In § 4.4, we showcase the utility of our method by studying advection–diffusion past multiple obstacles over a wide range of Péclet numbers. We conclude in § 5 by noting potential future extensions and applications of our results.
2. Mathematical formulation
Here, we formulate the problem of finding the incompressible potential flow comprising a uniform free stream past a collection of impermeable obstacles. We then interpret the complex potential of this flow as a conformal mapping from the physical domain to streamline coordinates, and demonstrate how the governing PDE for steady advection-diffusion is greatly simplified in the mapped domain.
2.1. Complex potential,
$W(z)$
, for flow past multiple obstacles
We now formulate the boundary value problem describing the incompressible potential flow around impermeable obstacles with a prescribed free stream (far-field) velocity. We consider flow past a collection of
$N^b$
many obstacles with boundaries given by smooth closed curves
$\partial B_{\!j}$
for
$j \in \{1,\ldots ,N^b\}$
. We write the free stream velocity as
$\boldsymbol{u}_\infty = (U \cos \varTheta , U \sin \varTheta )$
, where
$U\gt 0$
is the free stream velocity and
$\varTheta$
is the angle of the stream with respect to the
$x$
-axis. We non-dimensionalise lengths by
$L$
, a characteristic scale of the obstacles, and time by
$L/U$
. The non-dimensional potential flow problem is then to find a velocity
$\boldsymbol{u} = \boldsymbol{\nabla }\phi (x,y)$
such that
$\boldsymbol{\nabla }\boldsymbol{\cdot }\boldsymbol{u}=0$
,
$\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{n}_{\!j} = 0$
, where
$\boldsymbol{n}_{\!j}$
is the normal vector along the boundary of each obstacle and
$\boldsymbol{u}$
tends to
$(\cos \varTheta , \sin \varTheta )$
, with unit magnitude, at infinity.
The incompressibility condition implies
$\boldsymbol{\nabla} ^2 \phi = 0$
, so we proceed by considering the complexified problem with
$z = x + \mathrm{i} y$
. The flow potential can then be viewed as the real part of a complex potential
$W(z) = \phi (z) + \mathrm{i} \psi (z)$
, which is analytic in the fluid
$\mathbb{C}\backslash (\cup _{k=1}^{N^b}B_k )$
, and where the streamfunction
$\psi (z)$
is the harmonic conjugate of
$\phi (z)$
. The complex velocity can then be expressed as
$u = \overline {{\rm d}W/{\rm d}z}$
. The impermeability conditions can be expressed by requiring the streamfunction to be constant on the boundary of each obstacle. We thus seek a single-valued complex function
$W(z)$
that is analytic in the fluid and has the following properties:
\begin{eqnarray} \begin{gathered} W(z) \sim e^{-\mathrm{i}\varTheta } z,\; |z|\rightarrow \infty, \\ \mathrm{Im}\left \{W(z)\right \}=\psi _{\!j},\;z\in \partial B_{\!j}, \end{gathered} \end{eqnarray}
for some a priori unknown constants,
$\left \{\psi _1,\ldots , \psi _{N^b} \right \}$
. The flow is assumed to be free of singularities (e.g. source/sinks and point vortices) and further possesses zero-circulation around each obstacle. Writing
$W(z) = e^{-\mathrm{i} \varTheta } z + \delta W(z)$
, the free stream component can be subtracted and we can equivalently seek a complex analytic function
$\delta W(z)$
satisfying
\begin{eqnarray} \begin{gathered} \delta W(z) \rightarrow 0, \; |z| \rightarrow \infty, \\ \mathrm{Im}\left \{\delta W(z)\right \} = \psi _{\!j} - \mathrm{Im}\left \{e^{-\mathrm{i} \varTheta } z\right \}, \; z \in \partial B_{\!j}. \end{gathered} \end{eqnarray}
For non-trivial geometries,
$W(z)$
typically cannot be attained in a simple closed form. However, the problem for
$\delta W(z)$
can be solved numerically using a least-squares approach over sets of functions that are holomorphic in the fluid and decay at infinity. For instance, solutions can be found using truncated Laurent series centred within each obstacle. However, this approach can converge very slowly in some cases, for instance, when the boundaries are non-smooth, slender or concave. In these cases, the Laurent series ansatz can be supplemented with additional poles in the obstacle interiors; these poles can greatly enhance the solution accuracy, provided that their locations are chosen appropriately. We take the approach of using the AAA-LS algorithm (Costa & Trefethen Reference Costa and Trefethen2023) to automatically determine efficient pole placements for generic smooth obstacles. The numerical details of this approach will be discussed in § 3.1.1.
Two examples of potential flow solutions around obstacles and their interpretation as a conformal maps to streamfunction coordinates. (a) Free stream flow with
$\boldsymbol{u}_{\infty } = (1,0)$
around a circle. The complex potential is given by the classical Joukowsky map, bringing the unit circle to a slit along the real axis. (b) Free stream flow with
$\boldsymbol{u}_{\infty } = (1,0)$
around a generic smooth obstacle parametrised as
$z_{\!j}(\theta )$
. The complex potential is found numerically, but also maps the obstacle boundary to a horizontal slit in the mapped domain. A parametrisation of the obstacle, which we term the Joukowsky parametrisation
$\tilde {z}(\eta )$
, can be found whose image traverses the slit in a fashion like
$\cos \eta$
, analogous to the case of the circle under the classical Joukowsky map.
The complex potential governing the fluid flow,
$W(z)$
, may alternatively be interpreted as a conformal map between the region exterior to the obstacles in the physical domain and the entire complex plane exterior to a collection of slits, with the slits representing the images of
$\partial B_{\!j}$
under
$W(z)$
. To illustrate the geometry of the mapping, consider the flow past a unit disk with
$\varTheta =0$
, as is depicted in figure 1(a). Here, the mapping is given by the classical Joukowsky map
$W(z) = z + 1/z$
and the image of the fluid region is the entire complex plane exterior to the slit
$[-2,2]$
on the real axis, as is depicted in the top of figure 1(a). Points on the circle
$z=e^{{i}\eta }$
are mapped to a double-cover of the slit as
$w=2\cos \eta$
. The images of points just outside the circle
$z=(1+\epsilon ) e^{{i}\eta }$
traverse the slit in a positively oriented fashion starting from the trailing edge (
$w=2$
).
For a generic non-circular obstacle, the interpretation of
$W(z)$
as a conformal map is retained. The obstacle boundary, say given by some
$2 \pi$
-periodic parametrised curve
$z_{\!j}(\theta )$
, is mapped to a horizontal slit defined by the interval
$S_{\!j} = [\phi _{\!j}^{{min}}, \phi _{\!j}^{{max}}] + \mathrm{i} \psi _{\!j}$
. The length of the slit is
$2s_{\!j} \equiv \phi _{\!j}^{{max}}-\phi _{\!j}^{{min}}$
, and it is centred at
$\phi _{\!j} \equiv (\phi _{\!j}^{{max}}+\phi _{\!j}^{{min}})/2$
. A new parametrisation of the obstacle boundary
$\partial B_{\!j}$
, which we term the Joukowsky parametrisation
$\tilde {z}_{\!j}(\eta )$
, can be found such that its image under the mapping is
$W(\tilde {z}_{\!j}(\eta )) = \phi _{\!j} + s_{\!j} \cos \eta + \mathrm{i} \psi _{\!j} \equiv \phi _{\!j}(\eta ) + \mathrm{i} \psi _{\!j}$
. This is illustrated in figure 1(b). In multiply connected geometries, the same interpretation is retained and each obstacle is mapped to a separate horizontal slit, generally of different lengths and at different streamfunction values. A triply connected example is illustrated in figure 3.
2.2. Governing equation in streamline coordinates
Under the conformal mapping provided by the complex potential, the governing equation for steady advection–diffusion simplifies dramatically. Non-dimensionalising lengths by
$L$
, time by
$L/U$
and temperature by
$\Delta T$
, a characteristic scale for the variation of the obstacle temperatures from the far-field, (1.2) becomes
\begin{equation} {\textit{Pe}} \, \boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla }c = {\nabla} ^2 c, \end{equation}
where the Péclet number
${\textit{Pe}}=UL/\alpha$
and all quantities are non-dimensional. This equation can be interpreted as representing the advection and diffusion of
$c$
with unit diffusivity by a potential flow with a free stream velocity of amplitude
${\textit{Pe}}$
, in which case, it can be incorporated into the conformal map and hence the obstacle slits in the mapped domain will scale with
${\textit{Pe}}$
. Alternatively, the equation can be interpreted as representing the advection and diffusion of
$c$
with diffusivity
${\textit{Pe}}^{-1}$
by a potential flow with unit free stream velocity. In this case, the conformal map depends only on
$\varTheta$
, and the Péclet dependence remains present in the PDE in the mapped domain. We proceed with the latter approach for computational convenience.
Under the conformal mapping defined by
$W(z)$
, (2.3) can be rewritten in streamfunction coordinates as
\begin{equation} {\textit{Pe}} J c_\phi = J(c_{\phi \phi } + c_{\psi \psi }), \end{equation}
where
$J = |W'(z)|^2$
is the determinant of the Jacobian of the conformal map. Only when
$\boldsymbol{u}$
is a potential flow does the Jacobian factor from the advective term exactly match the Jacobian factor in the Laplacian under this change of coordinates (Bazant Reference Bazant2004). Since this Jacobian is non-zero in the fluid, where
$W(z)$
is holomorphic, the Jacobian factors can be cancelled and we are simply left with a constant-coefficient equation in the mapped domain:
\begin{equation} {\textit{Pe}}\, c_{\phi } = c_{\psi \psi }+c_{\phi \phi }. \end{equation}
The boundary conditions for
$c(\phi ,\psi )$
are now specified on the boundaries of each obstacle, which are horizontal slits in the mapped domain.
The constant-coefficient nature of (2.5) dramatically simplifies the PDE’s solution, since the free-space Green’s function of this operator, and its normal derivative, can be written explicitly as
\begin{equation} G(\phi , \psi ) =\frac {1}{2\pi } \exp \left(\frac {{\textit{Pe}}}{2} \phi \right) K_0\left(\frac {{\textit{Pe}}}{2}\sqrt {\phi ^2 + \psi ^2}\right) \end{equation}
and
\begin{equation} G_\psi (\phi , \psi ) = - \frac {{\textit{Pe}}}{4\pi } \exp \left(\frac {{\textit{Pe}}}{2} \phi \right) K_1\left(\frac {{\textit{Pe}}}{2} \sqrt {\phi ^2 + \psi ^2}\right) \frac {\psi }{\sqrt {\phi ^2 + \psi ^2}}, \end{equation}
where
$K_0$
and
$K_1$
are modified Bessel functions of the second kind (see, for example, Aichi Reference Aichi1920). Having simple known Green’s functions enables the efficient solution of the PDE using boundary integral methods. Importantly, we note that
$G$
is continuous across each slit (at
$\psi =0$
), while
$G_\psi$
is discontinuous. Thus, to solve problems with arbitrary boundary data, including different Dirichlet values on the top and bottom of each slit, we will adopt a combined single-plus-double-layer boundary integral formulation, where distributions of both (2.6) and (2.7) are included. For completeness, we note that if the Péclet scaling is instead incorporated into the conformal map, as described previously, the PDE and Green’s functions in the mapped domain become independent of Péclet number.
2.3. Variable-temperature boundary conditions
We consider advection–diffusion between obstacles with arbitrary Dirichlet temperature profiles on their boundaries. In the physical domain, the temperature conditions are specified on each obstacle boundary
$\partial B_{\!j}$
as
\begin{equation} c(z_{\!j}(\theta )) = c_{\!j}(\theta ), \end{equation}
or, alternatively, using the Joukowski parametrisations of each obstacle, as
\begin{equation} c(\tilde {z}_{\!j}(\eta )) = \tilde {c}_{\!j}(\eta ). \end{equation}
In the mapped domain, the boundary condition specification is slightly more subtle. Each boundary
$\partial B_{\!j}$
is mapped to a horizontal slit, with the top and bottom of the slit being the images of distinct points on the obstacle boundary in the physical domain. Since the Dirichlet data need not be constant in the physical domain, boundary values on each slit in the mapped domain can thus be different when approached from above or below. These top and bottom boundary conditions can be written simply using the Joukowsky parametrisation as
\begin{align} c(\phi _{\!j}(\eta ), \psi _{\!j} + \epsilon ) &= \tilde {c}_{\!j}(\eta ), \quad 0 \leqslant \eta \lt \pi, \end{align}
\begin{align} c(\phi _{\!j}(\eta ), \psi _{\!j} - \epsilon ) &= \tilde {c}_{\!j}(\eta ), \quad \pi \leqslant \eta \lt 2 \pi, \end{align}
for
$j \in \{1,\ldots ,N^b\}$
and where
$\epsilon \rightarrow 0^+$
. Because the temperature is generally discontinuous across each slit in the mapped domain, it is essential that we include a discontinuous double-layer contribution in our representation of
$c(\phi , \psi )$
, as will be developed in § 3.3.
3. Numerical formulation
We have seen that by moving to streamline coordinates, the equations governing conjugate heat transfer in a potential flow can be greatly simplified. Namely, the variable coefficient elliptic PDE (1.2) is transformed into the constant-coefficient PDE (2.5). In simply connected geometries, conformal maps (and thus potential flow solutions) are analytically known in a variety of geometries (Galante & Churchill (Reference Galante and Churchill1990)). In multiply connected problems, however, the set of analytically known conformal maps is much more limited. In the present section, we describe our numerical formulation and solution of multiply connected advection–diffusion problems. First, we compute the necessary conformal transformations by leveraging the AAA algorithm to efficiently compute the complex potential,
$W(z)$
, for flows around arbitrary numbers of obstacles. We then formulate a boundary integral solution to the advection–diffusion problem in the mapped domain. Finally, the solution is mapped backed to the physical domain via the inverse conformal mapping to
$W(z)$
.
3.1. Computation of conformal maps
The class of potential flow geometries for which the complex potential,
$W(z)$
, can be found in an analytic closed form is somewhat limited. In simply connected geometries, there exist a wide variety of known conformal maps from the unit circle to geometries of interest, such as ellipses and polygons (Nehari Reference Nehari2012), which thus allow for a closed-form representation of
$W(z)$
for those domains. However, in multiply connected geometries, the class of known conformal maps remains limited, despite recent developments that have unlocked a variety of mappings using the Prime function framework of Crowdy (Reference Crowdy2020). To address this issue, here we develop a general numerical method for deducing
$W(z)$
for an arbitrary collection of
$N^b$
obstacles via series solutions supplemented by specifically placed simple poles. We restrict our consideration to smooth obstacles that are well separated. In this regime, the potential flow can be efficiently solved to high precision using the AAA-LS algorithm as is described later. However, this method does not expedite convergence rates in close-contact problems. Extensions of the AAA-LS method to handle such problems is an open research problem and we do not attempt to address this issue in the present work. We thus restrict our attention to well-separated obstacles, specifically ones that are separated by at least their local radii of curvature.
3.1.1. Laurent series representation
The most straightforward representation of
$\delta W(z)$
is as a collection of truncated Laurent series centred within each obstacle and containing only negative powers to ensure that
$\delta W(z)$
is analytic in the fluid and decays at infinity. The form of
$W(z)$
is thus given by
\begin{equation} W(z)\approx \, e^{-\mathrm{i}\varTheta }z+\sum _{j=1}^{N^b}\sum _{k=1}^{N^L_{\!j}}\frac {a_{j,k}}{\left (z-z_{\!j}\right )^k}, \end{equation}
where
$z_{\!j}$
is a point centred within
$B_{\!j}$
and
$N^L_{\!j}$
is the degree of truncation of the Laurent expansion centred within the
$j{\mathrm{th}}$
obstacle. This representation can be used to find an approximation of
$W(z)$
by solving a least-squares problem for the series coefficients
$\{a_{j,k}\}$
by enforcing the boundary conditions
$\mathrm{Im}\{W(z)\}=\psi _{\!j}$
at
$N^s_{\!j}$
-many sample points
$z^s_{j,k}$
along each obstacle boundary
$\partial B_{\!j}$
. In practice, this least-squares solve converges best when
$N^s_{\!j} \sim 6 N^L_{\!j}$
and the system is solved using Arnoldi orthogonalisation (Brubeck, Nakatsukasa & Trefethen Reference Brubeck, Nakatsukasa and Trefethen2021).
In theory, when the full Laurent series is retained (
$N_{\!j}^L$
tends to infinity), (3.1) holds exactly in the fluid domain for any obstacle shapes. Unfortunately, when obstacles have corners, regions of high curvature or include concavities, such Laurent series representations often converge extremely slowly (Trefethen Reference Trefethen2020, Reference Trefethen2024). For instance, to attain a mapping
$W(z)$
to an accuracy of
$10^{-12}$
for the irregular shape in figure 1, the number of Laurent series coefficients that must be retained grows considerably to
$N_{\!j}^L=340$
. Such large
$N_{\!j}^L$
invariably makes solving for the series coefficients, as well as subsequent evaluations of
$W(z)$
, slow, so alternative strategies are needed for all but the simplest geometries.
3.1.2. Supplementing Laurent series with simple poles
For non-smooth or non-convex obstacles, supplementing the Laurent series with a set of strategically placed poles – a form of the method of fundamental solutions (MFS) – substantially improves the convergence of rational approximation of exterior harmonic functions. The key challenge is in finding a set of suitable pole locations. When obstacles have sharp corners, Gopal & Trefethen (Reference Gopal and Trefethen2019) showed that an exponentially clustered set of poles near corners resolves corner singularities with root-exponential convergence; this result is deeply connected to the work of Newman (Reference Newman1964). The AAA-LS algorithm (Costa & Trefethen Reference Costa and Trefethen2023) has recently emerged as a robust and automatic way to determine nearly ideal MFS poles for Laplace problems in both smooth and non-smooth domains. The AAA algorithm computes a rational barycentric approximation of a function defined on a contour in the complex plane (Nakatsukasa, Sète & Trefethen Reference Nakatsukasa, Sète and Trefethen2018). The poles of this rational function depend on the function being approximated as well as the contour on which the function is sampled. The AAA-LS algorithm computes the AAA approximation to some boundary data on
$\partial B_{\!j}$
and uses the interior poles of the resulting rational function as MFS support points. Here, we apply this algorithm to the Schwarz function,
$S(z) = \overline {z}$
for
$z \in \partial B_{\!j}$
, on obstacle boundaries (for more on approximation of the Schwarz function, see Trefethen (Reference Trefethen2025)). In practice, this method seems to automatically produce efficient pole placement locations for solving Laplace problems in many geometries.
To efficiently produce multiply connected potential flow solutions, we supplement low-order Laurent series within each obstacle with simple poles produced via the AAA-LS procedure. Although this method is able to accurately solve for potential flows around generic obstacles, here we only consider smooth obstacles to simplify our later boundary integral approach to the advection–diffusion problem. We determine the poles by applying AAA to the Schwarz function of each obstacle independently and retaining the
$N_{\!j}^p$
-many interior poles; we denote the collection of all such poles by
$\{z^p_k\}$
with
$k=1,\ldots ,N^p$
. Such poles are depicted in figures 2 and 3.
Steady-state temperature under advection–diffusion by a potential flow with
${\textit{Pe}}=1$
past a single obstacle with a varying prescribed surface temperature. (a) Solution in the
$w$
-domain, where the obstacle takes the form of a horizontal slit. A boundary-integral representation in this domain is used to solve the PDE. (b) Solution in the physical
$z$
-domain, produced by composing the conformal map with the boundary integral solution in the
$w$
-domain. White lines indicate the streamlines of the potential flow and the coloured line on the border of the obstacle indicates the prescribed temperature boundary condition. The red dots are the poles used to solve for the potential flow via the AAA-LS algorithm.
Steady-state temperature under advection–diffusion by a potential flow with
${\textit{Pe}}=1$
past three obstacles with various prescribed surface temperatures. (a) Solution in the
$w$
-domain, where each obstacle takes the form of a horizontal slit. A boundary-integral representation in this domain is used to solve the PDE. (b) Solution in the physical
$z$
-domain, produced by composing the conformal map with the boundary integral solution in the
$w$
-domain. White lines indicate the streamlines of the potential flow and the coloured lines on the borders of the obstacles indicate the prescribed temperature boundary conditions. The red dots are the poles used to solve for the potential flow via the AAA-LS algorithm.
The new AAA-LS representation of
$W(z)$
then takes the form
\begin{equation} W(z) \approx e^{-\mathrm{i}\varTheta }z+\sum _{j=1}^{N^b}\sum _{k=1}^{N^L_{\!j}}\frac {a_{j,k}}{\left (z-z_{\!j}\right )^k}+\sum _{k=1}^{N^p}\frac {b_k}{z-z^p_k}, \end{equation}
where now the coefficients
$\{a_{j,k}\}$
and
$\{b_k\}$
are simultaneously found through the solution of the associated least-squares problem. For example, to compute a mapping
$W(z)$
to an accuracy of
$10^{-12}$
for the irregular shape in figure 1, only
$N^L = 20$
Laurent series terms are needed in addition to the
$N^p=42$
simple poles with locations furnished by AAA. Without the AAA poles, as was mentioned earlier,
$N^L=340$
Laurent series terms must be retained to achieve the same accuracy.
While the poles found via AAA can efficiently represent data on each independent non-circular boundary, we have found that AAA poles do not help to resolve close-contact problems. Additionally, solving the AAA problem globally tends to produce lines of poles in the fluid region between nearby obstacles, but these poles are discarded and therefore still cannot help in the least-squares step. Future modifications of this approach, for instance, incorporating points from alternative fundamental-solution techniques (e.g. Cheng & Greengard Reference Cheng and Greengard1998) or placing higher-order poles at the AAA-LS simple pole locations, may extend its capabilities for close-contact problems, but we leave these explorations to future work.
3.1.3. Solving the least squares problem
To solve for the coefficients
$\{a_{j,k}\}$
and
$\{b_k\}$
in the representation of
$W(z)$
defined in (3.2), we substitute (3.2) into (2.1). Sampling the equation at points along the domain boundary defines a linear least squares problem for the coefficients. We denote the boundary sample points on the
$j{\mathrm{th}}$
obstacle boundary
$\partial B_{\!j}$
as
$\{z^s_{j,k}\}$
for
$k\in \{1,\ldots ,N^s_{\!j}\}$
, where
$N^s_{\!j}$
is the total number of sample points on
$\partial B_{\!j}$
. The constraint associated with the
$k{\mathrm{th}}$
sample point on the
$j{\mathrm{th}}$
obstacle may be expressed simply as
\begin{equation} \mathrm{Im}\left \{W(z^s_{j,k})\right \}\approx \psi _{\!j}. \end{equation}
These equations are imposed for each sample point on each obstacle, forming a linear least squares problem for the representation coefficients,
$\{a_{j,k}\}$
and
$\{b_k\}$
, along with the a priori unknown streamfunction values
$\{\psi _{\!j}\}$
. Performing Vandermonde with Arnoldi orthogonalisation on the Laurent series portion of the matrix prevents ill-conditioning of the associated Vandermonde blocks (Brubeck et al. Reference Brubeck, Nakatsukasa and Trefethen2021). For more details of the associated linear algebra problem, see Baddoo (Reference Baddoo2020). Therein, Baddoo places poles manually, whereas we use the AAA algorithm to determine the pole locations; however, the linear algebra problem remains the same. With the proper Arnoldi orthogonalisation, the system can be solved using standard linear least squares techniques. The error in the resulting map can be assessed simply by evaluating the residual
$|\mathrm{Im}\{W(z)\}-\psi _{\!j}|$
at additional points along each boundary
$z \in \partial B_{\!j}$
, and
$N_{\!j}^L$
can be increased until the convergence is suitable. In practice, we find using
$N^s_{\!j} \sim 6 (N^L_{\!j} + N^p_{\!j})$
works well.
3.2. Joukowsky parametrisation
The numerical conformal map
$W(z)$
solved previously maps each obstacle
$\partial B_{\!j}$
in the
$z$
-domain to a horizontal slit
$S_{\!j}$
in the
$w$
-domain, where the PDE has constant coefficients and can be efficiently solved using boundary integral techniques. As a first step in solving the boundary integral equation in the
$w$
-domain, it is useful to reparametrise the obstacles after computing the flow map
$W(z)$
. The obstacles are initially prescribed with some positively oriented shape and temperature parametrisations as
$z_{\!j}(\theta )$
and
$c_{\!j}(\theta )$
, which are generic parametrisations with no restrictions, other than that they are smooth
$2\pi$
-periodic functions of
$\theta$
.
The flow map
$W(z)$
maps the boundary of each object to a horizontal slit in the
$w$
-domain described by
$S_{\!j} = \{\phi + \mathrm{i} \psi _{\!j} : \phi ^{{min}}_{\!j} \leqslant \phi \leqslant \phi ^{{max}}_{\!j}\}$
. Motivated by the Joukowsky map of a circle to the slit, we seek a reparametrisation
$\theta _{\!j}(\eta )$
for each obstacle such that
\begin{equation} W(z_{\!j}(\theta _{\!j}(\eta ))) = \underbrace {\underbrace {\frac {\phi _{\!j}^{{min}}+\phi _{\!j}^{{max}}}{2}}_{\phi _{\!j}} +\underbrace {\frac {\phi _{\!j}^{{max}}-\phi _{\!j}^{{min}}}{2}}_{s_{\!j}} \cos (\eta )}_{\phi _{\!j}(\eta )} + \mathrm{i} \psi _{\!j}. \end{equation}
This takes the form of the image of the Joukowsky map of the circle parametrised as
$e^{i\eta }$
, stretched by
$s_{\!j}$
and shifted by
$\phi _{\!j}(0) + \mathrm{i} \psi _{\!j}$
. Importantly, it provides a parametrisation for the obstacle in the
$z$
-domain that conforms to the values of
$\phi$
on the boundary under the conformal map.
The parametrisation is unique and maintains the orientation of the original parametrisation if we specify
\begin{equation} \theta _{\!j}(\eta ) = \eta + \delta \theta _{\!j}(\eta ), \end{equation}
where
$\delta \theta _{\!j}(\eta )$
is a
$2\pi$
-periodic function of
$\eta$
, and take
$\theta _{\!j}(0) = \mathrm{argmax}_{\theta \in [0, 2\pi )}{} \mathrm{Re}\left \{ W(z_{\!j}(\theta ))\right \}$
. A high-precision Fourier approximation of
$\delta \theta _{\!j}(\eta )$
is found using ApproxFun.jl (https://github.com/JuliaApproximation/ApproxFun.jl) and rootfinding to determine
$\delta \theta _{\!j}(\eta )$
for various
$\eta$
. We use this function to compute the Joukowsky parametrisations,
$\tilde {z}_{\!j}(\eta ) = z_{\!j}(\theta _{\!j}(\eta ))$
and
$\tilde {c}_{\!j}(\eta ) = c_{\!j}(\theta _{\!j}(\eta ))$
, for each obstacle. Numerically, these are each represented with high-precision Fourier approximations using ApproxFun.jl.
3.3. Boundary integral formulation
With the flow-induced Joukowsky parametrisations at hand, we formulate the boundary integral representation for the temperature in the
$w$
-domain using a mixed single/double-layer representation:
\begin{equation} c(\phi , \psi ) = \sum _{j=1}^{N^b} S[f_{\!j},s_{\!j}](\phi -\phi _{\!j},\psi -\psi _{\!j}) + D[g_{\!j},s_{\!j}](\phi -\phi _{\!j},\psi -\psi _{\!j}), \end{equation}
where the single/double-layer operators are given by
\begin{align} S[f,s](\phi ,\psi ) &= \int _{-s}^s f(s') G(\phi -s', \psi ) \,{\rm d}s', \end{align}
\begin{align} D[g,s](\phi ,\psi ) &= \int _{-s}^s g(s') G_\psi (\phi -s', \psi ) \,{\rm d}s', \end{align}
and the Green’s function and its vertical/normal derivative are given by (2.6) and (2.7). We note that the single-layer operator
$S$
is continuous across
$\psi =0$
, while the double-layer operator
$D$
has limits
$D[g,s](\phi , 0^\pm ) = \mp (1/2) g(\phi )$
for
$|\phi | \lt s$
. We include both terms in our representation so that generic variable temperature conditions, corresponding to different limits approaching each slit from above and below, can be applied on the obstacles.
The single and double layer densities,
$f_{\!j}$
and
$g_{\!j}$
, are determined by requiring the temperature to approach the prescribed temperature boundary conditions
$c_{\!j}$
on the top and bottom of each slit:
\begin{align}& \lim _{\epsilon \rightarrow 0^+} c(\phi _{\!j} + s_{\!j} \cos (\eta ), \psi _{\!j} + \epsilon ) = \tilde {c}_{\!j}(\eta ), \quad 0 \leqslant \eta \lt \pi, \end{align}
\begin{align}& \lim _{\epsilon \rightarrow 0^+} c(\phi _{\!j} + s_{\!j} \cos (\eta ), \psi _{\!j} - \epsilon ) = \tilde {c}_{\!j}(\eta ), \quad -\pi \leqslant \eta \lt 0 . \end{align}
Subtracting these equations provides local equations determining the double-layer density on each slit, since the single-layer contributions (and double-layer contributions from other slits) are continuous across each slit:
\begin{align} \tilde {c}_{\!j}(\eta ) - \tilde {c}_{\!j}(-\eta ) &= \lim _{\epsilon \rightarrow 0^+} c(\phi _{\!j}(\eta ), \psi _{\!j} + \epsilon ) - c(\phi _{\!j}(\eta ), \psi _{\!j} - \epsilon ) \end{align}
\begin{align} &= D[g_{\!j},s_{\!j}](s_{\!j} \cos (\eta ), 0^+) - D[g_{\!j},s_{\!j}](s_{\!j} \cos (\eta ), 0^-) \end{align}
\begin{align} &= - g_{\!j}(s_{\!j} \cos (\eta )). \end{align}
That is, the double-layer densities
$g_{\!j}$
can be simply and individually determined by the odd part of
$\tilde {c}_{\!j}(\eta )$
for each obstacle, corresponding to the sine part of the function’s Fourier series. Adding the top and bottom conditions on each slit then produces coupled boundary integral equations for the single-layer densities:
\begin{align} \tilde {c}_{\!j}(\eta ) + \tilde {c}_{\!j}(-\eta ) &= \lim _{\epsilon \rightarrow 0^+} c(\phi _{\!j}(\eta ), \psi _{\!j} + \epsilon ) + c(\phi _{\!j}(\eta ), \psi _{\!j} - \epsilon ) \end{align}
\begin{align} &= 2 \sum _{j'} S[f_{j'},s_{j'}](\phi _{\!j}(\eta )-\phi _{j'}, \psi _{\!j}-\psi _{j'}) \end{align}
\begin{align} &\quad + 2 \sum _{j' \neq j} D[g_{j'},s_{j'}](\phi _{\!j}(\eta )-\phi _{j'}, \psi _{\!j}-\psi _{j'}). \end{align}
These are
$N^b$
-many coupled first-kind Fredholm integral equations for
$f_{\!j}$
, where the inhomogeneous term consists of the even (cosine) parts of
$\tilde {c}_{\!j}(\eta )$
minus the double-layer contributions from the other slits.
3.4. Boundary integral discretisation
From its identification with the sine part of
$\tilde {c}_{\!j}(\eta )$
, the double-layer densities can be written as
$g_{\!j}(s) = \tilde {g}_{\!j}(s) \sqrt {s_{\!j}^2 - s^2}$
, where
$\tilde {g}_{\!j}(s)$
is a smooth function on the interval
$[-s_{\!j},s_{\!j}]$
. From the theory of first-kind log kernel integral equations on open arcs (Atkinson & Sloan Reference Atkinson and Sloan1991), the single-layer densities can be written as
$f_{\!j}(s) = \tilde {f}_{\!j}(s) / \sqrt {s_{\!j}^2 - s^2}$
, where
$\tilde {f}_{\!j}(s)$
is a smooth function on the interval
$[-s_{\!j},s_{\!j}]$
. We therefore numerically represent each
$\tilde {f}_{\!j}$
and
$\tilde {g}_{\!j}$
as Chebyshev series over
$[-s_{\!j}, s_{\!j}]$
, and solve for them (via (3.11) and (3.12), respectively) using Chebyshev collocation. The single-layer integrals are discretised using a Lagrange polynomial expansion for
$f_{\!j}$
over the Chebyshev nodes, and the integrals are evaluated to high precision using QuadGK.jl (https://github.com/JuliaMath/QuadGK.jl). The resulting system is solved to determine the Chebyshev representation of
$f_{\!j}$
.
Once
$g_{\!j}$
and
$f_{\!j}$
have been determined, the temperature at any point in the
$w$
-domain can be computed using (3.6). The integrals are again efficiently evaluated using QuadGK.jl, which adaptively splits the integrand to resolve the singularities introduced by the kernel. The temperature at any point in the
$z$
-domain can be computed simply by applying the conformal map
$W(z)$
and then evaluating
$c$
at the corresponding position in the
$w$
-domain.
We demonstrate this solution procedure in figure 2, where we solve for advection–diffusion around a single obstacle with a varying surface temperature. The figure shows the solution in both the
$w$
-domain and the
$z$
-domain. Figure 3 shows the solution for potential flow around three obstacles, again in both the
$w$
and
$z$
domains.
3.5. Nusselt number computation
The net flux out of a given obstacle in the physical domain can be computed by integrating the normal component of the flux density,
$\boldsymbol{q}=-\alpha \boldsymbol{\nabla }c + c \boldsymbol{u}$
, around the boundary. Note the advective contribution is zero here since we are taking impermeable boundary conditions (
$\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{n} = 0$
) on each obstacle. Written in the physical domain, the flux out of the
$j{\mathrm{th}}$
obstacle is given by
\begin{equation} \mathcal{F}_{\!j} = \int _{\partial B_{\!j}} \boldsymbol{q} \boldsymbol{\cdot }\boldsymbol{n} \,{\rm d}l = - \alpha \int _{\partial B_{\!j}} \frac {\partial c}{\partial n} \,{\rm d}l, \end{equation}
where
$\boldsymbol{n}$
is the unit normal vector and
$n$
is the normal distance to the boundary, so that
$\partial c/\partial n=\boldsymbol{\nabla }c\boldsymbol{\cdot }\boldsymbol{n}$
. In streamline coordinates, because each boundary corresponds to a surface of constant
$\psi$
, we have for the top of the slit that
$ (\partial c / \partial n ) = |\partial \psi /\partial n|(\partial c/\partial \psi )=|J|^{(1/2)}(\partial c/\partial \psi )$
, whereas for the bottom of the slit, we have
$ (\partial c / \partial n ) =-|J|^{(1/2)}(\partial c/\partial \psi )$
. The distance along the boundary in the mapped domain transforms as
${\rm d}l=|J|^{-(1/2)}{\rm d}\phi$
. The Nusselt number associated with the
$j{\mathrm{th}}$
obstacle can then be written in terms of contributions from the top and bottom of the slit in the
$w$
-domain as
\begin{equation} {\textit{Nu}}_{\!j}= - \int _{\phi _{\!j}^{{min}}}^{\phi _{\!j}^{{max}}}\left (\frac {\partial c}{\partial \psi }\Big |_{\psi =\psi _{\!j}+\epsilon }-\frac {\partial c}{\partial \psi }\Big |_{\psi =\psi _{\!j}-\epsilon }\right )\,{\rm d}\phi . \end{equation}
Using the representation in (3.6), the double-layer contributions to the flux cancel and (3.14) further simplifies to
\begin{equation} {\textit{Nu}}_{\!j} = \int _{\phi _{\!j}^{{min}}}^{\phi _{\!j}^{{max}}} f_{\!j}(\phi -\phi _{\!j}) \,{\rm d}\phi = \int _{-s_{\!j}}^{s_{\!j}} f_{\!j}(s) \,{\rm d}s. \end{equation}
To verify the validity of (3.15), we consider the test case of two circles of equal radii,
$R$
, centred on the coordinates
$z_{\!j} = (-1)^j \mathrm{i} d$
and with temperatures
$c_{\!j}= (-1)^j \Delta T/2$
for
$j=\{1,2\}$
. In the absence of fluid flow, the net heat flux between the circles is given by
$\mathcal{F}=2\pi \alpha \Delta T/(\log { (1/\rho )})$
, where
$\rho =(d-\sqrt {d^2-R^2})^2/R^2$
. Thus, in the limit of
${\textit{Pe}} \rightarrow 0$
, the flux in and out of the top and bottom circles, respectively, should approach the conductive heat flux. We solve this problem numerically with
$d=2R$
and find precise agreement with error
$(|{\textit{Nu}}_2-\mathcal{F}|)\approx 10 \, {\textit{Pe}}^2$
for
$10^{-6} \leqslant {\textit{Pe}} \leqslant 10^{-2}$
.
4. Results
Here, we use the combined conformal mapping and boundary integral solver to produce solutions to a variety of variable-temperature and multiply connected conjugate heat transfer problems, where the primary result of interest is the heat flux in or out of each obstacle under different flow conditions. First, we examine the heat flux out of simply connected regions with variable temperature data at both large and small Péclet numbers. We then examine several multiply connected problems.
4.1. Joukowsky–Fourier decomposition of variable temperature data
To determine the heat flux from a single obstacle with variable temperature data, we first map the problem onto the canonical domain flow around a disk. For a single generic obstacle
$B_1$
in a uniform potential flow, there is a unique conformal map from the fluid region to the exterior of a disk that maps infinity to infinity with unit derivative. The radius of the resulting disk
$R$
is the conformal radius of
$B_1$
. Since this map is unique, it must be identical to the map formed by composing the conformal map
$W(z)$
with the inverse of the standard Joukowsky map stretched by
$R$
. The original heat transfer problem therefore corresponds identically to that of a disk of radius
$R$
with boundary data given by the Joukowsky parametrisation of the original boundary data, i.e.
$\tilde {c}_1(\eta )$
. We can therefore study the properties of any simply connected variable-temperature problem via its decomposition into Fourier modes on the disk given by the Joukowsky reparametrisation of the surface temperature. We can further restrict our study to Fourier modes on the unit disk, noting that the conformal radius
$R$
can be absorbed into a rescaling of the Péclet number when
$R \neq 1$
.
We therefore consider the canonical case of potential flow with
$\varTheta =0$
around a unit disk with boundary
$z_1(\theta ) = e^{\mathrm{i} \theta }$
and boundary temperature given by
$c_1(\theta )= \sum _{m=0}^\infty a_m \cos (m\theta ) + \sum _{m=1}^\infty b_m \sin (m\theta )$
. Since the heat transfer problem is linear in the Dirichlet data, we can consider each Fourier component independently. For each sine mode, the temperature data are odd under reflections across the real axis, meaning the temperature values along the top and bottom of the slit in the mapped domain are equal and opposite. Therefore, the single-layer density
$f_1=0$
for all sine components, and these cases contribute no net heat flux out of the obstacle. The cosine components are therefore the only part of the generic temperature distribution that contribute to the heat flux. In figure 4, we plot the heat flux out of the unit disk as a function of the Péclet number for several cosine temperature distributions. Power law scalings are apparent in both the low- and high-Péclet-number limits, and we will explore these further in the following sections. The temperature distributions for cosine and sine distributions with
$m=4$
are also shown, and it is visually apparent that the antisymmetric cancellation in the sine case prevents net heat transfer to infinity by the potential flow.
(a) Heat flux out of a unit disk with boundary temperature
$\cos {(m \theta )}$
at various
${\textit{Pe}}$
numbers, computed using our numerical boundary integral framework (coloured circles). The black line on the right indicates the high-
${\textit{Pe}}$
flux scaling
$\sim {\textit{Pe}}^{1/2}$
derived by Choi et al. (Reference Choi, Margetis, Squires and Bazant2005) for isothermal obstacles, (1.3). Our low-
${\textit{Pe}}$
asymptotic predictions for the flux given by (4.38) – and written explicitly in table 1– are plotted as solid black lines. The numerical results follow these scalings well until numerical errors appear past the adaptive quadrature precision level of
$10^{-12}$
. In panels (b) and (c), we plot the temperature field induced by boundary temperatures
$\cos {(4\theta )}$
and
$\sin {(4\theta )}$
with
${\textit{Pe}} = 1$
.
4.2. High-Péclet-number flux from single obstacles
Choi et al. (Reference Choi, Margetis, Squires and Bazant2005) find that the flux from a single isothermal obstacle scales as
${\textit{Pe}}^{(1/2)}$
, in the limit of
${\textit{Pe}} \rightarrow \infty$
. In agreement with this result, we have shown numerically in figure 4 that the flux from a circle with a cosinusoidal boundary temperature also follows the same scaling. For boundary data of the form
$\cos { (m\theta )}$
, we observe a scaling of
${\textit{Nu}} \sim C_m {\textit{Pe}}^{(1/2)}$
as
${\textit{Pe}}\rightarrow \infty$
; we have not determined the constants,
$C_m$
. We therefore expect this scaling to apply to generic temperature distributions, whose flux can be determined from the distribution’s cosine series by linearity, other than the exceptional cases where this leading order contribution cancels. Specifically, boundary data of the form
$\sum _{k=1}^{M}\alpha _k \cos { (k\theta )}$
should lead to a Nusselt number scaling of
${\textit{Pe}}^{(1/2)}\sum _{k=1}^{M}\alpha _k C_k$
. Thus, it is possible to cancel the leading order behaviour through an appropriate choice of coefficients
$\alpha _k$
, in which case, we expect the flux to scale like
${\textit{Pe}}^{-(1/2)}$
as
${\textit{Pe}}\rightarrow \infty$
. We note that further analytical progress may be possible to compute the exact limiting flux for each cosine mode, but we leave such calculations to future work.
4.3. Low-Péclet-number flux from single obstacles
In the case of low-Péclet-number flows, there are significant differences between variable temperature problems and the isothermal problem analysed by Choi et al. (Reference Choi, Margetis, Squires and Bazant2005). In the Laplace problem exterior to a unit circle,
$c(\infty )$
is simply the average of the Dirichlet data on the circle boundary. It follows that any boundary data on the circle possessing a non-zero average cannot satisfy the boundary condition
$c(\infty )=0$
. Therefore, the constant boundary data problem examined by Choi et al. (Reference Choi, Margetis, Squires and Bazant2005) has the property that the limit of
${\textit{Pe}}\rightarrow 0$
is singular: there is no Laplace solution in the case of
${\textit{Pe}}=0$
exactly. For our canonical cases with boundary data
$\cos (m\theta )$
on the unit disk with
$m \gt 0$
, a Laplace solution does exist when
${\textit{Pe}}=0$
and the limit
${\textit{Pe}}\rightarrow 0$
is not singular. We thus deduce that the heat flux in these cases must decay at least as fast as
$\textit {O} ({\textit{Pe}} )$
as
${\textit{Pe}}\rightarrow 0$
. This rate is very fast compared with the slow logarithmic decay of the flux for isothermal obstacles (1.4). In subsequent sections, we will confirm this prediction by deriving a multipole expansion for the Nusselt number in the
${\textit{Pe}}\rightarrow 0$
limit.
4.3.1. Boundary integral analysis
Recall that the Joukowsky transformation,
$W(z) \equiv \phi +\mathrm{i}\psi = (z+({1}/{z}) )$
, takes the circle
$e^{i \theta }$
in the physical domain to the horizontal slit defined by
$\phi \in [-2,2 ]$
and
$\psi =0$
, where
$\phi =2\cos (\theta )$
. Noting the identity
$\cos { (m \theta )}=T_m(\cos { (\theta )})$
, where
$T_m(x)$
is the Chebyshev polynomial of the first kind of degree
$m$
, the boundary condition
$\cos { (m \theta )}$
on the circle is translated to
$T_m(\phi /2)$
on the slit. Our goal is now to deduce the flux associated with each such mode. In the slit domain, the single-layer equation is
\begin{equation} \frac {1}{2\pi }\int _{-2}^2 e^{\frac {{\textit{Pe}}}{2}(\phi -s')}K_0 \left (\frac {{\textit{Pe}}}{2}|\phi -s'|\right ) f(s') \,{\rm d}s' = T_m\left (\frac {\phi }{2}\right ). \end{equation}
For convenience in computing forthcoming integrals, we rescale the slit to the interval
$(-1,1)$
by scaling the coordinates such that
$x = \phi /2$
and defining
$\sigma (s ) = 2 f(2s)$
. The integral equation then becomes
\begin{equation} \int _{-1}^1 e^{{\textit{Pe}}(x-s)}K_0({\textit{Pe}}|x-s|) \sigma (s) \,{\rm d} s\equiv \int _{-1}^1 \tilde {G}(x-s) \sigma (s)\, {\rm d} s = 2\pi T_m(x), \end{equation}
while the net flux becomes
$F_m =\int _{-1}^1\sigma (s)\,{\rm d}s$
. This kernel can be expanded for small
${\textit{Pe}}$
in the form
\begin{align} \tilde {G}&=\tilde {G}_0+{\textit{Pe}} \tilde {G}_1+{\textit{Pe}}^2 \tilde {G}_2+\textit {o}\left ({\textit{Pe}}^2\right )\!, \end{align}
\begin{align} \tilde {G}_0(x-s)&=-\log {\left (\frac {{\textit{Pe}}}{2}|x-s|\right )-\gamma }, \end{align}
\begin{align} \tilde {G}_1(x-s)&=\left (x-s\right )\left (-\log {\left (\frac {{\textit{Pe}}}{2}|x-s|\right )-\gamma }\right )\!, \end{align}
\begin{align} \tilde {G}_2(x-s)&=\frac {\left (x-s\right )^2}{4}\left (C-3\log {\left (|x-s|\right )}\right )\!, \end{align}
\begin{align} C &=1-3\gamma -3\log {\left (\frac {{\textit{Pe}}}{2}\right )}. \end{align}
Note that each of
$\tilde {G}_k$
is of order at most
$\log { ({\textit{Pe}} )}$
, so that each subsequent term in our kernel expansion is indeed asymptotically smaller than its predecessor. We proceed by solving for the flux density perturbatively in
${\textit{Pe}}$
, where
$\sigma ^{ (k )}(x)=\sigma _0(x)+{\textit{Pe}}\sigma _1(x)+{\cdots} + {\textit{Pe}}^k \sigma _k$
represents the asymptotic solution to the integral equation with the approximate kernel
$\tilde {G}_0+{\textit{Pe}} \tilde {G}_1+{\cdots} + {\textit{Pe}}^k \tilde {G}_k$
. In the following section, we consider the solution
$\sigma ^{ (1 )}{ (x )}=\sigma _0+{\textit{Pe}} \sigma _1$
when
$m=1$
, in which case the right-hand side of (4.2) becomes
$2\pi T_1(x)$
.
4.3.2. Leading order solution for
$T_1(x)$
data
In the leading order problem governing
$\sigma _0(x)$
, the kernel simplifies to
$\tilde {G}_0$
and the integral equation takes the form
\begin{eqnarray} \int _{-1}^1 \left (-\log {\left (\frac {{\textit{Pe}}}{2}|x-s|\right )-\gamma }\right )\sigma _0(s)\,{\rm d}s=2\pi T_1(x), \end{eqnarray}
which can be solved exactly as follows. We begin by expanding
$\sigma _0$
in a Chebyshev series,
\begin{eqnarray} \sigma _0(x)=\frac {\sum _{k=0}^{\infty }\alpha _k^{\left (0\right )} T_k(x)}{\sqrt {1-x^2}}, \end{eqnarray}
and make note of the following important identities:
\begin{align} \int _{-1}^1\frac {\log {\left (|x-s|\right )} T_{\!j}(s)}{\sqrt {1-s^2}}&=-\frac {\pi }{j}T_{\!j}{\left (x\right )}, \end{align}
\begin{align} \int _{-1}^1\frac {\log {\left (|x-s|\right )} }{\sqrt {1-s^2}}&=-\pi \log {\left (2\right )}, \end{align}
valid for
$j$
being a positive integer. The solution to (4.8) can be found by substituting our Chebyshev expansion into the integral equation, which then becomes
\begin{eqnarray} \pi \left (\log {\left (\frac {{\textit{Pe}}}{4}\right )}+\gamma \right )\alpha _0-\sum _{k=1}^{\infty }\frac {\pi }{k}\alpha _k T_k(x)=-2\pi T_1(x). \end{eqnarray}
We immediately deduce the coefficients by inspection as follows:
\begin{align} \alpha _{0}^{\left (0\right )}&=0, \end{align}
\begin{align} \alpha _1^{\left (0\right )}&=2, \end{align}
\begin{align} \alpha _{k\neq 1}^{\left (0\right )}&=0, \end{align}
so that
$\sigma _0$
is given explicitly by
\begin{equation} \sigma _0(x)=\frac {2T_1(x)}{\sqrt {1-x^2}}. \end{equation}
A key observation is that
$\sigma _0$
has identically zero flux,
$\int _{-1}^{1}\sigma_{0}(x)\,{\rm d}x=0$
, since
$\alpha ^{ (0 )}_0=0$
in its Chebyshev expansion. This means that flux can only appear at higher order in
${\textit{Pe}}$
.
4.3.3. First-order solution for
$T_1(x)$
data
The solution at order
${\textit{Pe}}$
is found by considering now the kernel
$\tilde {G}_0+{\textit{Pe}} \tilde {G}_1$
. Seeking the solution to the equation of the form
$\sigma ^{(1)}=\sigma _0+{\textit{Pe}} \sigma _1$
, the governing equation for
$\sigma _1$
becomes
\begin{equation} \int _{-1}^1\tilde {G}_0(x-s)\sigma _1(s)\,{\rm d}s=-\int _{-1}^1\tilde {G}_1(x-s)\sigma _0(s)\,{\rm d}s, \end{equation}
where the right-hand side of the equation can be computed using the known form of
$\sigma _0$
from (4.16). We again proceed by expanding
$\sigma _1$
as a Chebyshev series,
\begin{eqnarray} \sigma _1(x)=\frac {\sum _{k=0}^{\infty }\alpha _k^{\left (1\right )} T_k(x)}{\sqrt {1-x^2}}, \end{eqnarray}
and our goal will be to deduce the coefficients
$\alpha _k^{ (1 )}$
. After performing integrals and manipulations that have been relegated to Appendix A, we find that the integral equation reduces to the following form:
\begin{align} &-\pi \left (\log {\left (\frac {{\textit{Pe}}}{4}\right )}+\gamma \right )\alpha _0+\sum _{k=1}^{\infty }\frac {\pi }{k}\alpha _k T_k(x)\nonumber\\ &=-\frac {\pi }{2}T_2(x)-\pi \left (\log {\left (\frac {{\textit{Pe}}}{4}\right )}+\gamma +1\right )T_0(x). \end{align}
By inspection, we deduce the coefficients of the expansion for
$\sigma _1$
as follows:
\begin{align} \alpha _0^{\left (1\right )}&=\frac {\left (\log {\left (\dfrac {{\textit{Pe}}}{4}\right )}+\gamma +1\right )}{\left (\log {\left (\dfrac {{\textit{Pe}}}{4}\right )}+\gamma \right )}, \end{align}
\begin{align} \alpha _1^{\left (1\right )}&=0, \end{align}
\begin{align} \alpha _2^{\left (1\right )}&=-1, \end{align}
\begin{align} \alpha _{k\gt 2}^{\left (1\right )}&=0, \end{align}
leading to the following expression for the flux density,
\begin{equation} \sigma ^{\left (1\right )}(s)= \frac {2T_1(s)+{\textit{Pe}}\left (\dfrac {\left (\log {\left (\frac {{\textit{Pe}}}{4}\right )}+\gamma +1\right )}{\left (\log {\left (\frac {{\textit{Pe}}}{4}\right )}+\gamma \right )}T_0(s)-T_2(s)\right )}{\sqrt {1-s^2}}. \end{equation}
The flux is then given by the integral of
$\sigma ^{ (1 )}(x)$
over the interval
$(-1,1)$
, such that the flux at this order is
\begin{equation} F_{m=1}\sim\pi {\textit{Pe}} \frac {\left (\log {\left (\dfrac {{\textit{Pe}}}{4}\right )}+\gamma +1\right )}{\left (\log {\left (\dfrac {{\textit{Pe}}}{4}\right )}+\gamma \right )}. \end{equation}
Note that when the logarithmic correction is ignored, the flux (4.25) reduces to
$\pi {\textit{Pe}}$
. In practice, the logarithmic terms are significant unless
${\textit{Pe}}$
is very small. We proceed by computing the flux in the case of the boundary condition
$T_2(x)$
.
4.3.4. Leading order solution for
$T_2(x)$
data
By an analysis identical to that presented in § 4.3.2, we find the flux density at leading order for boundary data
$T_2(x)$
on the slit is given by
\begin{equation} \sigma _0(x)=\frac {4T_2(x)}{\sqrt {1-x^2}}. \end{equation}
4.3.5. First-order solution for
$T_2(x)$
data
The equation governing the first-order solution is then,
\begin{equation} \int _{-1}^1\tilde {G}_0(x-s)\sigma _1(s)\,{\rm d}s=-\int _{-1}^1\tilde {G}_1(x-s)\sigma _0(s)\,{\rm d}s, \end{equation}
which we solve by expanding
$\sigma _1$
in a Chebyshev series and explicitly computing the right-hand side of (4.27), the details of which can be found in Appendix B. After comparing the Chebyshev coefficients in the expansion of
$\sigma _1$
and the terms from the explicit computation of the right side of (4.27), we find
\begin{equation} \sigma _1(x)=\frac {T_1(x)-T_3(x)}{\sqrt {1-x^2}}. \end{equation}
Again, we observe no flux at this order of the equation. To find the flux which appears at
$\textit {o} ({\textit{Pe}} )$
, we must examine the next order kernel equation.
4.3.6. Second-order solution for
$T_2(x)$
data
The second-order kernel equation takes the following form:
\begin{equation} \int _{-1}^1\tilde {G}_0(x-s)\sigma _2(s)\,{\rm d}s=\underbrace {-\int _{-1}^1\tilde {G}_1(x-s)\sigma _1(s)\,{\rm d}s}_{J_{2,1}}+\underbrace {-\int _{-1}^1\tilde {G}_2(x-s)\sigma _0(s)\,{\rm d}s}_{J_{2,0}}, \end{equation}
where the integrals
$J_{2,1}$
and
$J_{2,0}$
can be computed explicitly since we have already computed
$\sigma _0$
and
$\sigma _1$
. Instead of computing both sides of (4.29) as a function of
$x$
, we can integrate both sides of the equation against
$1/\sqrt {1-x^2}$
to obtain the flux,
\begin{equation} -\pi \left (\log {\left (\frac {{\textit{Pe}}}{4}\right )}+\gamma \right )F_{m=2}\sim\int _{-1}^1\frac {J_{2,1}(x)}{\sqrt {1-x^2}}\,{\rm d}x+\int _{-1}^1\frac {J_{2,0}(x)}{\sqrt {1-x^2}}\,{\rm d}x. \end{equation}
Performing the integration over
$x$
first greatly simplifies the integral computation once one exploits the orthogonality of the Chebyshev polynomials in the subsequent integrals over
$s$
. Details of the integral computations are suppressed for brevity. Such integrals yield the following leading order expression for the total flux in the case of
$m=2$
:
\begin{eqnarray} F_{m=2}\sim -\frac {\pi {\textit{Pe}}^2}{4}\left (\frac {\log {\left (\dfrac {{\textit{Pe}}}{4}\right )}+\gamma +\dfrac {3}{2}}{\log {\left (\dfrac {{\textit{Pe}}}{4}\right )+\gamma }}\right )\!, \end{eqnarray}
in the limit of
${\textit{Pe}}\rightarrow 0$
.
4.3.7. Generic boundary conditions
$T_m(x)$
Expressions for the leading order flux in the case of boundary data
$T_m(x)$
with
$m\gt 2$
may be computed in a perturbative manner similar to our analysis which yielded the fluxes when
$m=1$
and
$m=2$
. However, the perturbation equations become increasingly cumbersome as
$m$
is increased. By studying and exploiting the structure of the perturbation equations, we can derive a closed form expression for the flux valid for any
$m\gt 0$
.
The first important deduction is that when the boundary data are given by
$T_m(x)$
, the flux scales as
${\textit{Pe}}^m$
asymptotically as
${\textit{Pe}}\rightarrow 0$
. This fact can be seen as follows. The leading order equation has solution for
$m\gt 0$
of the form,
\begin{equation} \sigma _0(x)=\frac {2m T_m(x)}{\sqrt {1-x^2}}. \end{equation}
The next order equation can also be solved for arbitrary
$m$
and one finds for
$m\gt 1$
,
\begin{equation} \sigma _1(x)=\frac {T_{m-1}\left (x\right )-T_{m+1}\left (x\right )}{\sqrt {1-x^2}}. \end{equation}
In general, by examining the form of the
$k{\mathrm{th}}$
-order integral equation and noting the identity
$x T_k(x)= (T_{k-1} (x )+T_{k+1} (x ) )/2$
, it can be seen that Chebyshev polynomials are introduced at each order in a manner such that
\begin{equation} \sigma _k(x)=\frac {\sum _{p=m-k}^{m+k}\alpha _p^{\left (k\right )} T_p(x)}{\sqrt {1-x^2}}, \end{equation}
where the specific values of the coefficients
$\alpha _p^{ (k )}$
may be deduced by performing the asymptotic analysis to order
${\textit{Pe}}^k$
. The key takeaway from this observation is that the flux vanishes at all orders
${\textit{Pe}}^k$
for
$k\lt m$
, since the flux associated with
$\sigma _k$
is simply given by
$ \alpha _0^{ (k )} \pi {\textit{Pe}}^k$
. Structurally, it can be seen that the
$\alpha _{m}^{ (0 )}$
coefficient feeds the
$\alpha _{m-1}^{ (1 )}$
coefficient in the
$\sigma _1$
equation. Furthermore, the
$\alpha _{m-1}^{ (1 )}$
coefficient feeds
$\alpha _{m-2}^{ (2 )}$
, and generally the
$\alpha _{m-k+1}^{ (k-1 )}$
coefficient feeds the
$\alpha _{m-k}^{ (k )}$
for
$k\lt m$
. By observing this fact, it is possible to derive the following recursive summation expression for
$\alpha _{m-k}^{ (k )}$
in terms of terms
$\alpha _{m-k+j}^{ (k-j )}$
with
$j \geqslant 1$
:
\begin{align} \alpha _{m-k}^{\left (k\right )}&=-\sum _{j=1}^{k}\frac {\left (-1\right )^{j}C_{\!j} }{2^j}\frac {j! \left (m-k\right )!}{\left (m-k+j\right )!}\alpha _{m-k+j}^{\left (k-j\right )};\quad 0\lt k\lt m, \end{align}
\begin{align} \alpha _m^{\left (0\right )}&=2m, \end{align}
\begin{align} C_{\!j}&=\frac {1}{ j!} {}_2F_1\left (\frac {1-j}{2},-\frac {j}{2};1;1\right )\!, \end{align}
where
${}_2F_1$
is the generalised hypergeometric function, the details of the derivation being relegated to Appendix C. Furthermore, by examining the kernel expansion equation at order
${\textit{Pe}}^m$
, it can be shown that the flux associated with
$\sigma _m$
is only fed by one term in each of the lower-order flux densities,
$\sigma _{k\lt m}$
. Specifically, the flux of
$\sigma _m$
depends only on the terms associated with the coefficients
$\alpha _{m-k}^{ (k )}$
, where
$0\leqslant k \lt m$
. It is thus possible to show that the flux associated with the boundary condition
$T_m(x)$
is given by
\begin{align} F_{m}&=-\frac {\sum _{k=0}^{m-1} \dfrac {\left (-1\right )^{m-k}}{2^{m-k}}\left (C_{m-k}\left (\log {\left (\dfrac {{\textit{Pe}}}{4}\right )}+H_{m-k}\right )+\tilde {C}_{m-k}\right )\alpha ^{\left (k\right )}_{m-k}}{\log {\left (\dfrac {{\textit{Pe}}}{4}\right )}+\gamma }\pi {\textit{Pe}}^m,\end{align}
\begin{align} H_{\!j}&=\sum _{k=1}^j\frac {1}{k},\end{align}
\begin{align} \tilde {C}_{k}&=-\sum _{j=0}^{\left \lfloor \frac {k}{2} \right \rfloor }\frac {1}{\left (k-2j\right )!}\frac {\psi \left (j+1\right )}{\left (j!\right )^2 2^{2j}},\end{align}
where
$\psi$
is the Digamma function and
$H_{\!j}$
is the
$j{\mathrm{th}}$
harmonic number. Again, the details of the derivation can be found in Appendix C. We numerically evaluate the coefficients in (4.38) and tabulate the results for the cases
$m\in \{1,2,3,4\}$
in table 1. For these values of
$m$
, we have verified the expressions against our numerical results and find precise agreement. Specifically, we observe relative errors between our numerically computed fluxes and the scaling predictions of
$\lesssim 10\, {\textit{Pe}}^2$
over the range
$10^{-6} \leqslant {\textit{Pe}} \leqslant 10^{-2}$
and until the predicted fluxes are below our selected quadrature precision of
$10^{-12}$
.
Leading order expressions for the flux in the case of boundary conditions
$\cos {(m \theta )}$
on the unit circle, as a function of Péclet number, in the asymptotic limit
${\textit{Pe}}\rightarrow 0$
, as computed using (4.38).
It is worthwhile to emphasise the following interesting aspect of the analysis that led to (4.38). The flux associated with the
$m{\mathrm{th}}$
order perturbation equation was obtained without having to explicitly compute
$\sigma _k$
at all lower orders
$k\lt m$
. This was made possible through the examination of only the structure of the perturbation problem at each order and noticing a recursive relationship between a particular coefficient at each order of the problem.
4.4. Multiply connected configurations
(a) Numerically computed heat fluxes from three obstacles with varying surface temperatures at various Péclet numbers. The dashed lines indicate the written scalings. In panels (b) and (c), we plot the temperature field for
${\textit{Pe}}=10^{-3}$
and
${\textit{Pe}}=10^3$
. At low Péclet number, the obstacles exchange heat and generally behave like a collective monopole in the far field. At high Péclet number, the obstacles are thermally isolated as their wakes occupy separate streamlines.
First column shows the numerical heat fluxes from two stream-aligned disks with constant surface temperatures at various Péclet numbers. In panel (a), the disks have the same temperatures,
$(c_1, c_2) = (1, 1)$
. In panel (b), the disks have opposite temperatures,
$(c_1, c_2) = (1, -1)$
. The dashed lines indicate the written scalings. The centre and right columns show the temperature fields for
${\textit{Pe}}=10^{-2}$
and
${\textit{Pe}}=10^2$
. At low Péclet number, the top configuration behaves like a collective monopole, while the bottom configuration approaches a Laplace solution and behaves like a collective dipole. At high Péclet number, the obstacles all emit a flux following the singly connected scaling of
${\textit{Pe}}^{\frac {1}{2}}$
. The leading order total flux at high Péclet number can be cancelled by choosing
$c_2/c_1 \approx -0.467$
for this geometry.
In multiply connected configurations, such as that presented in figure 3, obstacles can interact so as to partially cancel or enhance the heat flux out of other obstacles. Compared with the simply connected case that has been studied by Choi et al. (Reference Choi, Margetis, Squires and Bazant2005), this produces a rich new set of inter-obstacle interactions to be studied which may have interesting implications in growth dynamics problems such as ADLA. Placing an obstacle downstream of another, as compared with it being placed transversely relative to the free stream direction, can particularly complicate the relative and total heat fluxes at high Péclet number. Nonetheless, scaling regimes remain apparent in the low- and high-
${\textit{Pe}}$
limits in multiply connected problems.
Figure 5 depicts the fluxes at various Péclet numbers, as well as example temperature fields, for three generic obstacles occupying different streamlines and with varying surface temperatures. This represents a generic case of interacting obstacles without any particular symmetries or alignment. We see that at low Péclet number, the fluxes from each obstacle asymptote towards constant values which cancel at leading order, and their sum asymptotes towards a multiple of the monopole (
$m=0$
) solution which scales like
$-1/\log ({\textit{Pe}})$
. For particular surface temperature distributions, for which a Laplace solution exists (with zero temperature at infinity), this leading contribution cancels for low
${\textit{Pe}}$
and the flux would follow a power law scaling in Péclet. At high Péclet number, the flux from each obstacle and the total flux all scale as
${\textit{Pe}}^{1/2}$
, as in the simply connected case examined by Choi et al. (Reference Choi, Margetis, Squires and Bazant2005). Physically, we can see in figure 5(c), which depicts the temperature field for
${\textit{Pe}} = 10^3$
, that heat from each obstacle is confined to thin non-interacting wakes behind each obstacle. In the high-Péclet-number limit, the obstacles become thermally isolated from other obstacles on different streamlines.
In figure 6, we consider two unit disks centred at
$z_1=-2$
and
$z_2=2$
, aligned with the potential flow (
$\varTheta =0$
). Since the heat transfer problem is linear in the imposed obstacle temperatures, we solve the two canonical isothermal problems with
$(c_1,c_2) = (1,1)$
, depicted in panel (a), and
$(c_1,c_2)=(1,-1)$
, depicted in panel (b). When the disks have the same temperature, they behave as a collective monopole (
$m=0$
) at low Péclet number. When the disks have opposite temperatures, there exists a Laplace solution, so the leading total flux is cancelled and they behave as a collective dipole (
$m=1$
) at low Péclet number. At high Péclet number, all of the fluxes scale as
${{\textit{Pe}}}^{(1/2)}$
. However, there exists a certain temperature ratio (for this geometry,
$c_2/c_1 \approx -0.467$
) for which the leading order fluxes will cancel. Typically, however, even when obstacles occupy the same streamline, the total flux will follow the
${{\textit{Pe}}}^{(1/2)}$
scaling.
5. Discussion and conclusion
In this paper, we developed an efficient numerical procedure for solving advection–diffusion problems involving a uniform potential flow past a collection of impermeable obstacles with generic boundary temperature distributions. Our work thus extends prior studies, which examined advection–diffusion in flows past a single obstacle with isothermal boundary conditions. Our method leverages the AAA-LS method for solving Laplace problems (Costa & Trefethen Reference Costa and Trefethen2023) to efficiently compute the complex potential of the corresponding potential flow, which doubles as a conformal mapping to our computational domain. In the mapped domain, comprising a collection of horizontal slits in an otherwise infinite plane, we formulate a mixed boundary integral solution to the heat transfer problem with arbitrary Dirichlet data on each boundary. We thus enable an efficient solution to advection–diffusion in potential flows past arbitrary collections of obstacles. We note that our numerical approach is fast and highly accurate when obstacles are smooth and not too close in proximity to one another.
We then applied this method to examine the heat flux from obstacles in various heat-transfer regimes. In the limit of high Péclet number, we demonstrate in figure 4(a) that the flux exiting each obstacle scales as
${\textit{Pe}}^{1/2}$
, in agreement with the scaling for a simply connected isothermal obstacle given in (1.4). A notable exception occurs in the simply connected problem for boundary conditions of the type
$\sin (m\eta )$
in the Joukowsky parametrisation, in which case,
$\mathcal{F}=0$
independent of Péclet number.
We showed numerically that the Nusselt number at low Péclet number, with variable boundary data, no longer scales in the same manner as the simply connected isothermal problem (1.4). We then derived asymptotic Nusselt number scalings in the limit of low Péclet number for flow past a circle with arbitrary boundary data. In particular, we showed that boundary data of the form
$\cos (m\theta )$
on the unit circle produces a net flux in proportion to
${\textit{Pe}}^m$
for
$m \geqslant 1$
. Expressions for the leading coefficient of the scaling were derived for all cases
$m\geqslant 1$
using a perturbative approach leading to (4.38), with explicit expressions for the cases
$m\in \{1,2,3,4\}$
given in table 1. These expressions include logarithmic corrections and are found to agree to high precision with our numerical results, specifically with relative errors
$\approx {\textit{Pe}}^2$
for
${\textit{Pe}} \leqslant 10^{-2}$
and above our set precision limit (
$|{\textit{Nu}}| \gt 10^{-12}$
). We have thus shown that, in the limit of small
${\textit{Pe}}$
, the flux out of an arbitrary obstacle is determined by its Chebyshev representation in the corresponding slit domain or, equivalently, its cosine series under the Joukowsky parametrisation.
Our results might prove useful in a variety of future studies of advection–diffusion within potential flows. Because we employ a boundary integral approach, our solution procedure is highly accurate and is substantially faster than volume-based methods. In the future, we hope to extend our results to more general types of potential flows, which might include circulation and flow singularities in addition to a free stream. Such extensions would make possible investigations of (1.2) within realistic flows such as the lifting flow around a wing and other aerodynamic applications. Other interesting directions might include the investigation of growth, melting and artificial freezing problems, which have mainly been analysed in the context of simply connected domains with isothermal (constant Dirichlet) boundary conditions. It would be interesting to see how fluid dynamically mediated interactions between obstacles might influence their growth.
Another potential application is in Hele-Shaw flows and microfluidic systems, where a common objective is to promote mixing. For example, one might study the evolution of solutes within Hele-Shaw flows; one might be interested in the concentration of important biological species in lab-on-a-chip technologies or possibly the evolution of a simple dye. Our method may also be useful in the design of cold plate systems meant to optimise heat transfer, in systems where the approximation of potential flow is reasonable, as our method obviates the need for solving the full conjugate heat transfer equations or resorting to surrogate models during optimisation. Again, extensions of the present work to generalised potential flows – e.g. flows with circulation or point sources/sinks/vortices – would be advantageous. While some work has investigated steady advection–diffusion in point source/sink flows, the effects of obstacles have not yet to be studied (Boulais & Gervais Reference Boulais and Gervais2020). Magnetohydrodynamic Hele-Shaw flows (McKee Reference McKee2024) represent singularity-free potential flows which, in general, comprise both a free stream as well as circulation around obstacles. Since one possible application of such flows is in promoting mixing, monitoring advection diffusion processes is certainly of interest. Thus, we intend to extend our present framework to treat such flows in the future.
Acknowledgements
The authors are grateful for many valuable discussions with N. Trefethen, D. Fortunato, A. Barnett, L. Greengard, S. Olver, R. Rosales, T. Driscoll, D. Crowdy, H. Zhu and M. Bazant.
Funding
K.M. was generously funded by a MathWorks fellowship during this work. This work was supported by a grant from the Simons Foundation (SSRFA-6826, K.J.B.).
Declaration of interests
The authors report no conflict of interest.
Data availability statement
The simulation code is available on GitHub at https://github.com/kburns/ConjugateHeatTransfer.jl.
Appendix A. Integrals from
$\boldsymbol{T}_{\boldsymbol{1}}\boldsymbol{(x)}$
perturbation analysis
The right-hand side of (4.17), which we define by
$J_{1,0}= (-\int _{-1}^1\tilde {G}_1(x-s)\sigma _0(s)\,{\rm d}s )$
, is written explicitly as follows:
\begin{equation} J_{1,0} =2\int _{-1}^1 \left (x-s\right )\left (\log {\left (|x-s|\right )}+\log {\left (\frac {{\textit{Pe}}}{2}\right )+\gamma }\right ) \frac {T_1(s)}{\sqrt {1-s^2}}\,{\rm d}s. \end{equation}
The first term of the expression becomes
\begin{align} &2\int _{-1}^1 \frac {\left (x-s\right )T_1(s)}{\sqrt {1-s^2}}\log {\left (|x-s|\right )}\,{\rm d}s=-2\pi x T_1(x)-2\int _{-1}^1\frac {T_1(s)T_1(s)\log {\left (|x-s|\right )}}{\sqrt {1-s^2}}\,{\rm d}s\nonumber\\ &=-\pi \left (T_0\left (x\right )+T_2\left (x\right )\right )+\pi \left ( \log {\left (2\right )}T_0(x)+\frac {1}{2}T_{2}\left (x\right )\right )\nonumber\\ &=\pi \left (\log {\left (2\right )}-1\right )-\frac {\pi }{2}T_2(x), \end{align}
and the second term becomes
\begin{eqnarray} 2\left (\log {\left (\frac {{\textit{Pe}}}{2}\right )}+\gamma \right )\int _{-1}^1 \frac {\left (x-s\right )T_1(s)}{\sqrt {1-s^2}}\,{\rm d}s=-\pi \left (\log {\left (\frac {{\textit{Pe}}}{2}\right )}+\gamma \right )T_0(x), \end{eqnarray}
so that in sum,
\begin{eqnarray} J_{1,0}=-\frac {\pi }{2} T_2(x)-\pi \left (\log {\left (\frac {{\textit{Pe}}}{4}\right )}+\gamma +1\right )T_0(x). \end{eqnarray}
Appendix B. Integrals from
$\boldsymbol{T}_{\boldsymbol{2}}\boldsymbol{(x)}$
perturbation analysis
B.1. Integrals appearing in first order equation
After noting that all
$s$
-independent terms of the kernel vanish due to the orthogonality of Chebyshev polynomials, the right-hand term of (4.27) becomes
\begin{align} 4\int _{-1}^1\left (x-s\right )\frac {T_2(s)\log {\left (|x-s|\right )}}{\sqrt {1-s^2}}\,{\rm d}s&=-4x\frac {\pi }{2}T_2(x)+2\pi \left (\frac {1}{3}T_3(x)+T_1(x)\right )\nonumber\\ &=-\pi \left (T_1(x)+T_3(x)\right )+2\pi \left (\frac {1}{3}T_3(x)+T_1(x)\right )\nonumber\\ & =\pi \left (-\frac {1}{3}T_3(x)+T_1(x)\right )\!. \end{align}
Appendix C. Derivation of flux equation for
$\boldsymbol{T}_{\boldsymbol{m}}(\boldsymbol{x})$
data
C.1. Determination of expression for
$\alpha _{m-k}^{(k)}$
for
$k\lt m$
We begin by examining the
$k{\mathrm{th}}$
-order integral equation that governs the flux density
$\sigma _k$
, which takes the form
\begin{equation} \int _{-1}^1\tilde {G}_0(x-s)\sigma _k(s)\,{\rm d}s=-\sum _{p=1}^k\int _{-1}^1\tilde {G}_p(x-s)\sigma _{k-p}(s)\,{\rm d}s, \end{equation}
where we note that
$\sigma _p$
for
$p\lt k$
are known, at this order of perturbation theory, via the solution of lower order perturbation equations. We proceed by considering the Chebyshev expansions of
$\sigma _{\!j}$
for all
$j$
and isolating the terms associated with
$T_{m-k}$
on both sides of the equation so that
$\alpha _{m-k}^{ (k )}$
may be deduced. To this end, for the remainder of this section of the paper, we shall suppress all terms that do not contribute to terms of the form
$T_{m-k}$
.
Let us isolate each term in the right-hand side of the equation by defining
\begin{equation} J_{k,k-p}\equiv -\int _{-1}^1 \tilde {G}_p(x-s)\sigma_{k-p}\,{\rm d}s, \end{equation}
so the right-hand side of (C1) is equal to the sum of all the
$J_{k,k-p}$
integrals over the range
$1\leqslant p \leqslant k$
. Writing out the kernel explicity gives the following expression for the portion of
$J_{k,k-p}$
that can contribute terms proportional to
$T_{m-k} (x )$
, which we denote
$I_{k,k-p}$
:
\begin{equation} I_{k,k-p} = \int _{-1}^1\left (x-s\right )^p\left (C_p \log {\left (\frac {{\textit{Pe}}}{2}|x-s|\right )}+\tilde {C}_p\right )\frac {\alpha _{m-k+p}^{\left (k-p\right )} T_{m-k+p}\left (s\right )}{\sqrt {1-s^2}}\,{\rm d}s, \end{equation}
where
$C_p$
and
$\tilde {C}_p$
represents the
$p{\mathrm{th}}$
coefficients of the following expansions:
\begin{align} e^{{\textit{Pe}} \, x} I_0{\left ({\textit{Pe}} \, x\right )}&=\sum _{p=0}^{\infty }C_p {\textit{Pe}}^p x^p, \end{align}
\begin{align} -e^{{\textit{Pe}} \, x} \sum _{j=0}^{\infty } \frac {\psi (j+1)}{(j!)^2} \left (\frac {{\textit{Pe}}}{2} x\right )^{2j} &= \sum _{p=0}^{\infty }\tilde {C}_p{\textit{Pe}}^p x^p, \end{align}
which are valid for small
${\textit{Pe}}$
, and where
$I_0$
is the modified Bessel function of the first kind of order zero. The forms of
$C_p$
and
$\tilde {C}_p$
are also summarised in the main text in (4.37) and (4.40). Once again, we note that (C3) only contains the contribution to the integral that is relevant in producing terms of the type
$T_{m-k} (x )$
.
Now, since
$m-k\geqslant 1$
, all terms in the integral not associated with the
$\log (|x-s| )$
term must vanish by the orthogonality of Chebyshev polynomials of unequal indices. It follows that
$I_{k,k-p}$
reduces to the following form:
\begin{equation} I_{k,k-p} = C_p \alpha _{m-k+p}^{\left (k-p\right )} \int _{-1}^1 \frac { \left (x-s\right )^p T_{m-k+p}\left (s\right )\log {\left (|x-s|\right )}}{\sqrt {1-s^2}}\,{\rm{d}}s. \end{equation}
Every term in the expansion of
$ (x-s )^p$
is capable of generating terms of the form
$T_{m-k}(x)$
through integration. By expanding
$ (x-s )^p$
and employing integral identities of the Chebyshev polynomials, one can show that
\begin{align} I_{k,k-p}&=C_p \alpha _{m-k+p}^{\left (k-p\right )} \sum _{l=0}^{p}\frac {\left (-1\right )^l p! }{l! \left (p-l\right )!}x^{p-l}\int _{-1}^1\frac { s^l T_{m-k+p}\left (s\right )\log {\left (|x-s|\right )}}{\sqrt {1-s^2}} \,{\rm{d}}s\end{align}
\begin{align} &=C_p \alpha _{m-k+p}^{\left (k-p\right )} \sum _{l=0}^{p}\frac {\left (-1\right )^l p! }{l! \left (p-l\right )!}x^{p-l}\frac {-\pi }{2^l \left (m-k+p-l\right )}T_{m-k+p-l}\left (x\right ) \end{align}
\begin{align} &=-\frac {\pi C_p T_{m-k}\left (x\right ) \alpha _{m-k+p}^{\left (k-p\right )}}{2^p} \sum _{l=0}^{p}\frac {\left (-1\right )^l p! }{l! \left (p-l\right )!}\frac {1}{ \left (m-k+p-l\right )}+ H.O.T., \end{align}
where the last summation can be evaluated explicitly and is given as
$ (-1 )^p p! (m-k-1 )!/ (m-k+p )!$
and where
$H.O.T.$
only includes Chebyshev polynomials with degrees greater than
$m-k$
. Noting that the relevant term in the left-hand side of (C1) reduces to
$(\pi / (m - k))\alpha _{m-k}^{ (k )} T_{m-k} (x )$
, we find the following expression for
$\alpha _{m-k}^{ (k )}$
:
\begin{align}\alpha _{m-k}^{\left (k\right )}&=-\sum _{p=1}^{k}\frac { C_p \alpha _{m-k+p}^{\left (k-p\right )}}{2^p}\frac {\left (-1\right )^p p! \left (m-k\right )!}{\left (m-k+p\right )!}\alpha _{m-k+p}^{\left (k-p\right )};\quad 0\lt k\lt m,\end{align}
\begin{align} \alpha _{m}^{\left (0\right )}&=2m,\end{align}
the latter equation setting the base case of the recursive relation for the coefficients
$\alpha _{m-k}^{ (k )}$
with
$k\gt 0$
. This expression can be evaluated symbolically using just a few lines of code in standard software such as Mathematica.
C.2. Determination of flux at order
${\textit{Pe}}^m$
The
$\textit {O} ({\textit{Pe}}^m )$
equations read as follows:
\begin{equation} \int _{-1}^1\tilde {G}_0(x-s)\sigma _m(s)\,{\rm d}s=\sum _{k=0}^{m-1}\underbrace {\int _{-1}^1{-\tilde {G}_{m-k}(x-s)\sigma _{k}(s)\,{\rm d}s}}_{J_{m,k}}, \end{equation}
where the portion of
$J_{m,k}$
that may contribute to the flux, which we denote
$I_{m,k}$
, takes the form
\begin{equation} I_{m,k} = \int _{-1}^1\left (C_{m-k} \log {\left (\frac {{\textit{Pe}}}{2}|x-s|\right )}+\tilde {C}_{m-k}\right )\frac {\alpha _{m-k}^{\left (k\right )} T_{m-k}\left (s\right )\left (x-s\right )^{m-k}}{\sqrt {1-s^2}}\,{\rm{d}}s, \end{equation}
which can be decomposed as follows:
\begin{align} I_{m,k}&=C_{m-k}\alpha _{m-k}^{\left (k\right )} I'_{m,k}+\left (\tilde {C}_{m-k}+C_{m-k}\log {\left (\frac {{\textit{Pe}}}{2}\right )}\right ) \alpha _{m-k}^{\left (k\right )} I''_{m,k}, \end{align}
\begin{align} I_{m,k}^{\prime} &= \int _{-1}^1\log {\left (|x-s|\right )}\frac { T_{m-k}\left (s\right )\left (x-s\right )^{m-k}}{\sqrt {1-s^2}}\,{\rm{d}}s, \end{align}
\begin{align} I_{m,k}^{\prime\prime} &= \int _{-1}^1\frac { T_{m-k}\left (s\right )\left (x-s\right )^{m-k}}{\sqrt {1-s^2}}\,{\rm{d}}s. \end{align}
By the orthogonality of Chebyshev polynomials,
$I_{m,k}''$
vanishes for every term in the expansion of
$ (x{-}s )^{m-k}$
except the
$s^{m-k}$
term such that
\begin{equation} I_{m,k}^{\prime\prime}=\frac {\pi }{2^{m-k}}\left (-1\right )^{m-k}, \end{equation}
whereas all the terms in
$I_{m,k}^{\prime}$
will need to be summed. Writing out the expansion of
$I_{m,k}^{\prime}$
explicitly, we find that
\begin{equation} I_{m,k}^{\prime}=\sum _{p=0}^{m-k}\frac {\left (-1\right )^{p}\left (m-k\right )! x^{m-k-p}}{p! \left (m-k-p\right )!}\int _{-1}^1\log {\left (|x-s|\right )}\frac { T_{m-k}\left (s\right )s^p}{\sqrt {1-s^2}}\,{\rm{d}}s, \end{equation}
the
$T_0$
component of which is
\begin{align} \frac{1}{\pi} \int_{-1}^{1} \frac{I'_{m,k}(x)}{\sqrt{1-x^2}} {\rm{d}}x=-\frac {\pi }{2^{m-k}} \sum _{p=0}^{m-k-1}\frac {\left (-1\right )^{p}\left (m-k\right )!}{p! \left (m-k-p\right )!}\frac {1}{m-k-p}-\frac {\left (-1\right )^{m-k}\pi \log {2}}{2^{m-k}}. \end{align}
The first sum can be computed exactly as
$ (-1 )^{m-k+1} H_{m-k}$
, where
$H_k$
is the
$k{\mathrm{th}}$
harmonic number. We thus find that
\begin{equation} \frac{1}{\pi} \int_{-1}^{1} \frac{I'_{m,k}(x)}{\sqrt{1-x^2}} {\rm{d}}x=\frac {\pi \left (-1\right )^{m-k}H_{m-k}-\left (-1\right )^{m-k}\pi \log {2}}{2^{m-k}}. \end{equation}
Summing this integral and
$I_{m,k}^{\prime\prime}$
, we find that the leading Chebyshev coefficient of the right-hand side of (C12) is given by
\begin{equation} \sum _{k=0}^{m-1} \frac {\pi \left (-1\right )^{m-k}}{2^{m-k}}\left (C_{m-k}\left (\log {\left (\frac {{\textit{Pe}}}{4}\right )}+H_{m-k}\right )+\tilde {C}_{m-k}\right )\alpha _{m-k}^{\left (k\right )}. \end{equation}
After comparing this constant with the left-hand side of (C12), we find that
\begin{equation} \alpha _{0}^{\left (m\right )}=-\frac {\sum _{k=0}^{m-1}\frac { \left (-1\right )^{m-k}}{2^{m-k}}\left (C_{m-k}\left (\log {\left (\frac {{\textit{Pe}}}{4}\right )}+H_{m-k}\right )+\tilde {C}_{m-k}\right )\alpha _{m-k}^{\left (k\right )}}{\log {\left (\frac {{\textit{Pe}}}{4}\right )}+\gamma }. \end{equation}
The flux is then simply given by
$\pi \alpha _{0}^{ (m )} {\textit{Pe}}^m$
.
