Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Penguin Huddling: A Continuum Model

Penguin Huddling: A Continuum Model Penguins huddling in a cold wind are represented by a two-dimensional, continuum model. The huddle boundary evolves due to heat loss to the huddle exterior and through the re- organisation of penguins as they seek to regulate their heat production within the huddle. These two heat transfer mechanisms, along with area, or penguin number, conservation, gives a free boundary problem whose dynamics depend on both the dynamics interior and exterior to the huddle. Assuming the huddle shape evolves slowly compared to the advective timescale of the exterior wind, the interior temperature is governed by a Poisson equation and the exterior temperature by the steady advection-diffusion equation. The exterior, ad- vective wind velocity is the gradient of a harmonic, scalar field. The conformal invariance of the exterior governing equations is used to convert the system to a Polubarinova-Galin type equation, with forcing depending on both the interior and exterior temperature gradi- ents at the huddle boundary. The interior Poisson equation is not conformally invariant, so the interior temperature gradient is found numerically using a combined adaptive Antoulas- Anderson and least squares algorithm. The results show that, irrespective of the starting shape, penguin huddles evolve into an egg-like steady shape. This shape is dependent on the wind strength, parameterised by the Péclet number Pe, and a parameter β which effectively measures the strength of the interior self-generation of heat by the penguins. The numerical method developed is applicable to a further five free boundary problems. Keywords Free boundary problems · Conformal mapping · AAA-least squares algorithm · Fluid mechanics · Collective behaviour · Huddling 1Introduction Among all the penguin species living in the Antarctic region [1], the emperor penguin (Aptenodytes forsteri) experiences the most extreme weather conditions, enduring wind −1 ◦ speeds over 30 ms at temperatures below −40 C[1–4]. Unlike other species such as the gentoo and the Adélie [5], emperor penguins have no fixed nest [6, 7], allowing the birds to huddle together for warmth in severe wind conditions. This strategy is very effective; ambient temperatures within the huddle often exceed 20 C, with a maximum temperature of 37.5 C having been recorded [8]. Such huddling is crucial for the penguins’ survival, S.J. Harris sam.harris.16@ucl.ac.uk Department of Mathematics, University College London, Gower Street, London, WC1E 6BT, UK 7 Page 2of17 S.J.Harris, N.R. McDonald especially during foraging trips [9] and throughout their breeding season, which coincides with the Antarctic winter [3]. Over the few hours that huddling events typically last [6, 8], the penguins are contin- uously (albeit slowly) on the move. The birds constantly reorganise their position in the huddle, where [1] observes that “birds that at first are in the center of the huddle become members of the rear flank [the side most exposed to the wind] and move, in their turn, up the sideline”. This process benefits all members of the huddle and allows the heat to be equally shared: those at the edge do not remain cold and those at the centre do not overheat. These individual penguin movements affect the propagation of the huddle as a whole - the huddle is seen to propagate downwind [10] and the shape traced out by the edge of the huddle, or huddle boundary, may change over time [4, 6]. It is the propagation and evolution of the huddle boundary shape that is the focus of this work. The modelling of penguin huddle dynamics has received increased attention over the last decade, but still remains sparse in the literature, despite numerous observations and findings from the field, see [1–3, 10]. The works [4, 7, 11] draw the analogy between penguin hud- dle/colony formation and condensed matter physics, with Lennard-Jones-like forces used to represent attractive and repulsive interactions between penguins that allow complex, lattice structures (such as huddles) to form. This idea is expanded upon in [12], where Gaussian processes are used to test if robots could adopt this penguin inspired huddling strategy. Sim- ilarly, [13] liken the rearrangement of the penguin huddle to a phase transition from solid (dense huddle) to liquid (loose huddle) to gas (no huddle) depending on environmental fac- tors such as temperature and wind speed. A rigorous fluid dynamics approach is taken in [14], where a finite difference method is used to solve the full Navier-Stokes equations of the surrounding fluid (the exterior wind), giving a two-way coupling between the environ- ment and the behaviour of the huddle. Using a simpler fluids based model, [6] assume the exterior wind to be a two-dimensional, irrotational flow of an incompressible fluid and from this model, the temperature profile exterior to the huddle and individual penguin movements can be found. In this work, the approach of [6] is expanded upon with two key differences. First, a simple model for the effect of interior penguin reorganisation as a result of self-generation of heat throughout the entire huddle, as observed in [1], is included. In [6], while penguin reorganisation around the huddle boundary was acknowledged, it was assumed that interior penguins were so tightly packed that they could not change their position in the huddle. That assumption is here dropped, yet it is still assumed that the penguins are packed tightly enough such that wind does not permeate through the huddle and so the wind flow remains entirely in the exterior region. Throughout this work, the phrases “exterior only penguin problem” and “full penguin problem” are used to refer to the penguin huddle problem with- out and with interior penguin reorganisation, respectively. Second, the penguin huddle itself is treated as a continuum. Existing huddle models treat penguins as discrete elements (or particles) such that individual penguin movements can easily be tracked. For sufficiently large huddles on the order of a thousand penguins [1], these discrete models quickly become computationally expensive. Continuum models, as used here, are often faster, cheaper, and preferred when studying the change in shape and position of the overall huddle. Using a continuum to model the interior and exterior of the full penguin problem draws analogues to other two-dimensional free boundary problems. Foremost among these is the classical Hele-Shaw free-boundary problem and its extensions, including Saffman-Taylor fingering (a consequence of injecting a fluid into another, more viscous fluid) [15–17], the expansion and contraction of bubbles in a Hele-Shaw cell [18–21] and the dissolution of soluble objects in two-dimensional potential flow [22]. Using conformal mapping to re- formulate the free boundary problem as a Polubarinova-Galin (PG) equation has proved Penguin Huddling: A Continuum Model Page 3 of 17 7 effective in finding exact and numerical solutions for this class of problems. Solidification and melting problems have also been extensively studied using analogous techniques, see [23–28]and [29]; the latter of these examines the freezing of porous media flow by freeze pipes, with the velocity of the free boundary dependent on both interior and exterior temper- ature gradients, similar to the full penguin problem. Free boundary problems arise in other natural processes, including, but not limited to, dendritic crystal growth [30, 31]and flame front propagation [32–34]. Continua may also be used to model other collective behaviours [4, 35, 36] such as fish schooling [37], bird flocking [38], bat swarming [39, 40]and ant colony formation [41, 42]. There are two heat transfer mechanisms in the full penguin problem, which formulate a suitable boundary condition for the propagation of the huddle boundary. The first is the cooling effect of the wind on the huddle. As in [6], the incompressible wind flow is assumed irrotational and flowing sufficiently faster than the propagation of the huddle, thus, the flow can be treated as steady. As such, the wind flow and temperature transport in the region exterior to the penguin huddle can be modelled by the Laplace and the steady advection- diffusion equations, respectively. This coupled pair of partial differential equations (PDEs) is conformally invariant [43] - a property which motivates the use of a conformal mapping approach in the numerical method [21–24, 34]. The exterior region and the huddle bound- ary in the physical z−plane can be conformally mapped to the exterior and boundary of the unit disk in some “mathematical” (or “canonical”) ζ −plane, and vice versa. Therefore, the physical coordinates can be expressed in terms of the conformal map z = f(ζ, t) and the condition on the huddle boundary rewritten as a PG type equation [44]. The task then becomes that of finding the conformal map. The second heat transfer mechanism is that of the uniform, self-generation of heat by the penguins in the interior. Its steady diffusion amongst the densely packed penguins is modelled by a Poisson equation for the interior temperature with a constant forcing term. By the Riemann mapping theorem (see e.g. [43]), the interior of the penguin huddle will not also be mapped to the interior of the ζ −disk by the conformal map z = f(ζ, t) used to map the exteriors. Further, the Poisson equation is not conformally invariant [45]. Hence, to solve for temperature in the interior using a conformal mapping approach, a different conformal map f would need to be found numerically and the transformed Poisson equation in the ζ − plane solved. Alternatively, the interior temperature can be found directly in the z−plane without the use of conformal mapping. Here, the methods developed by Trefethen and colleagues [46–50] are used to find a rational approximation of a harmonic function in some simply connected domain as the real part of some analytic function, subject to given boundary conditions. This function consists of a power series in z and a series of rational functions involving poles in the exterior clustered near singularity points [47]. Both the adaptive Antoulas-Anderson algorithm (AAA, pronounced triple-A) [49] and a least- squares (LS) algorithm, such as the lightning Laplace solver [48] can be used to numerically find this analytic function [50] and hence the required harmonic function. In this work, by first finding a particular solution, the problem of solving the interior Poisson equation then becomes one of solving the Laplace equation with a modified boundary condition. In turn, this is readily solved numerically by the AAA-LS method. The structure of this paper is as follows. The two-dimensional, continuum penguin hud- dle model is formulated in Sect. 2, with assumptions stated and dimensionless quantities introduced. A suitable boundary condition for the normal boundary velocity v is derived, incorporating effects of exterior wind cooling, internal heat generation by the penguins and area (penguin) conservation. The numerical method is then developed in Sect. 3,where the equation for v is rewritten as a PG type equation using a conformal mapping approach. n 7 Page 4of17 S.J.Harris, N.R. McDonald Attention here is given to the exterior heat flux σ in the ζ −plane and the interior heat flux which is found in the physical plane. Results are presented in Sect. 4, illustrating the evo- lution of penguin huddles of different starting shapes under a variety of wind and interior heating effect strengths. Conclusions are discussed in Sect. 5 and extensions to the work, and applications to further (free boundary) problems, noted. 2ModelSetup Consider a huddle of penguins with horizontal length scale L where the penguins are as- sumed homogeneous in size and shape. The huddle is considered as a continuum with pen- guins packed as tightly as possible while still allowing for huddle reorganisation [11]and hence uniform spacing between each penguin can be assumed constant across the entire huddle and for all time. Emperor penguins have a typical height of H = 1m [1]and their tight packing means that background winds are diverted around the huddle, allowing the dy- namics to be considered as two-dimensional, see Fig. 1. The boundary of the penguin huddle γ(t) is defined as the two-dimensional, non-overlapping curve enclosing the huddle; the fi- nite, simply connected interior of the curve occupied by the huddling penguins is labelled R and the penguin-free exterior as . The boundary curve γ is orientated anti-clockwise such that the unit normal n ˆ points from the interior to the exterior. The huddle R is buffeted by a unidirectional and constant strength wind of magnitude U which, without loss of generality, flows in the positive x ˆ direction. The wind is the irrota- tional flow of air, which is considered to be an inviscid and incompressible fluid which flows sufficiently rapidly compared to the evolution of the penguin huddle as a whole, so that the transport of temperature can be treated as steady [6]. In reality, viscous effects taking the form of thin boundary layers around the exterior of the penguin huddle are present. In this work, these are ignored in order to facilitate the accurate computation of the exterior po- tential flow using conformal mapping based methods. Denote the temperature far from the huddle as T , the temperature on the huddle boundary as T and the maximum temperature ∞ f in the huddle interior as T . The equations governing the wind and exterior temperature transport are, respectively, the Laplace equation and the steady advection-diffusion equation ∇ φ =0in , (1) u ·∇T = D ∇ T in , (2) where φ is the wind velocity potential, u = ∇φ is the wind velocity, T is the temperature in  and D is the diffusivity of temperature in air. The associated boundary and far field Fig. 1 Schematic diagram of the full penguin problem showing the interior region R and exterior region . a) Plan view where T denotes the maximum temperature T inside the huddle, b) side view c Penguin Huddling: A Continuum Model Page 5 of 17 7 conditions are ∂φ T = T , =0on γ, (3) ∂n T → T ,φ → Ux as r →∞. (4) Now, consider the interior R of the huddle where there is no wind. Owing to their ability to self-generate heat, which subsequently diffuses within the huddle, penguins nearer the huddle centre are warmer than those at the boundary [1, 11], i.e. T >T . This source of c f heat provides a further contribution to the heat flux at the boundary γ . Consistent with the continuum assumption, it is assumed that heat is generated uniformly across the huddle. Further, as in the exterior, it is assumed the interior heat diffuses quickly compared to the evolution of the boundary. Thus the temperature distribution in R is modelled by the steady heat equation ∇ T =− in R, (5) where D is the diffusivity of the interior temperature T through the penguin huddle and R R Q is the (constant) source term representing the self-generation of heat by the penguins. In order to obtain parameter estimates for the interior problem, it is helpful to first con- sider a circular penguin huddle of radius L, noting that T (0) = T . The Poisson equation R c (5) can then be solved exactly, giving Q Q 4(T − T ) c f T = T(r) = T − r , where = . (6) R c 4D D L R R The relation (6)isusedinSect. 2.1 to non-dimensionalise the problem for cases when the region R is non-circular. Finally, a condition on γ(t) is sought by specifying its outward normal velocity v .Fol- lowing [6], the idea that the boundary evolves owing to the heat loss of penguins on the boundary is used. Consequently, the problem has some similarity to the two-phase Stefan problem [23, 24, 29, 51] with the boundary temperature gradient driving its evolution, yet, unlike such problems, no phase change occurs. Waters et al. [6] move individual boundary penguins according to their heat loss to the exterior, whereas here both exterior and interior heat fluxes are considered. First, the heat loss to the exterior wind is the flux κ n ˆ ·∇T eval- uated on γ(t),where κ is the thermal conductivity of air. Second, the heat gain due to the warm huddle interior is −κ n ˆ ·∇T evaluated on γ(t),where κ is the thermal conductivity R R R of huddling penguins. Equating the normal velocity with the net heat flux gives Av = κ n ˆ ·∇T − κ n ˆ ·∇T , on γ(t), (7) n   R R where the constant term A is chosen to ensure a match in dimensions between both sides of (7), with the particular choice of A affecting the timescale for the boundary evolu- tion. For convenience in the non-dimensionalisation process upcoming in Sect. 2.1,let A = λρ c (T − T ),where ρ and c are the density and specific heat capacity of the f ∞ air, respectively, and λ is a non-dimensional constant to be chosen later. The right hand side of (7) is the jump in heat flux at the boundary. This observation reinforces the analogy with the Stefan problem where typically the normal velocity of the boundary depends on the jump in the heat flux [51]. For example, analogous boundary con- ditions arise in other natural free boundary problems, including the dissolution of finite 7 Page 6of17 S.J.Harris, N.R. McDonald objects, where only the gradient of the agent undergoing advection and diffusion exterior to the object is important [22, 24], or in porous media flow about freeze pipes, in which the temperature gradient in both the exterior fluid and interior ice regions are significant [29]. However it is emphasised that, in this work, (7) is purely a modelling assumption. Rewriting (7) gives the following condition on γ ∂x λρ c (T − T ) n ˆ · = κ n ˆ ·∇T − κ n ˆ ·∇T + C(t). (8) f ∞   R R ∂t The additional term C (t) in (8) is needed to enforce conservation of penguins, or, equiva- lently, area conservation [6], namely that ∂x n ˆ · ds = v ds = 0. (9) ∂t γ γ 2.1 Non-dimensionalisation and Parameter Values In order to identify the parameters of the problem, the model (1)-(5), (8) is non- dimensionalised, with non-dimensional quantities temporarily labelled as starred variables. Let L characterise the length scale of region R i.e. x = Lx , and use the far-field wind magnitude U as the characteristic velocity scale, u = U u . The following scalings are also used 1 T − T 1 U ∗ ∗ ∗ ∇ = ∇ ,T = ,t = t ,τ = . (10) L T − T τ L f ∞ Using (6), the non-dimensional form of the Poisson equation (5) can be written (now dropping the stars) as 4(T − T ) c f ∇ T =−α, α = . (11) T − T f ∞ The boundary condition (8) becomes ∂x λPe n ˆ · = n ˆ ·∇T − βn ˆ ·∇T + C(t), (12) ∂t where Pe = UL/D is the Péclet number, with D = κ /ρ c the thermal diffusivity of air, β = κ /κ is the ratio of thermal conductivities and C(t) is the dimensionless area conserving term, an explicit expression for which can be found by substituting (12)into(9) (see (13g) below). Note that the appearance of Pe on the left hand side of (12) affects the time scaling reported in the results. The real world dimensional time scale “t ” is therefore t = λPe (L/U )t . In what follows, the choice λ = 1 is made and justified a posteriori by comparing the time for the huddle to evolve with observations - see Sect. 5. There are three parameters: Pe, β and α. Following the justification in [6], the wind around the penguin huddle is likely turbulent and hence the Péclet number can be approxi- mated by a turbulent Reynolds number, which is O(100) for Antarctic winds at this scale. Consistent with [6], the present work considers Péclet values between 1 and 100. The values of the remaining two parameters, α and β , are estimated as follows. The air temperature far upstream is taken to be T =−40 C[2, 3] while ambient temperatures in the huddle are above freezing [11]. It seems reasonable then to let T = 0 C. The centre of a huddle can f Penguin Huddling: A Continuum Model Page 7 of 17 7 ◦ ◦ exceed 20 C, with a maximum temperature of 37.5 C having been recorded [8]. Using, therefore, T = 30 Cgives α = 3. From table C.5 of [52], the thermal conductivity of air at −40 C is approximately κ ≈ −1 −1 −2 −1 0.022Wm K . The thermal conductance of penguins is 1.3Wm K from [1], however this is not the same as the thermal conductivity - it has different units for example. To find the thermal conductivity, the conductance must be multiplied by the thickness of the −1 −1 penguin insulation layer which is d = 0.024m [53]. Therefore, κ = 0.0312Wm K and so β = 1.4182 ≈ 1.5 = O(1). This estimate for β depends on the wind flow being either laminar or turbulent, with β decreasing as the flow becomes more turbulent and the effective thermal conductivity of air increases [54]. In this work, a range of β values between β = 0 −1 and β = 1.5 are considered with β = 0.14 ≈ 0.1 = O(10 ) taken to be a representative value. To summarise, the full dimensionless system for the penguin problem is Pe v = n ˆ ·∇T − βn ˆ ·∇T + C(t), on γ, (13a) n  R ∇ φ =0in , (13b) Pe u ·∇T =∇ T in , (13c) ∇ T =−α in R, (13d) ∂φ T = T = 1, =0on γ, (13e) ∂n T → 0,φ → x as r →∞, (13f) − (n ˆ ·∇T − βn ˆ ·∇T )ds C(t) = . (13g) ds Thetaskistosolve thesystem(13a)–(13g) and, in particular, find the evolution of the γ(t) for some given initial huddle shape γ(0). This is done in Sect. 4 with the following choices: Pe = 1 → 100, β = 0 → 1.5and α = 3. 3 Numerical Method Following the work of [22, 29], the free boundary problem (13a)–(13g) can be written as a PG type equation [43] in terms of a conformal map z = f(ζ, t) from the exterior of the unit disk in the canonical ζ −plane to the exterior region  in the physical z−plane [6, 26]. However, the interior R is not mapped using f(ζ, t), thus care must be taken with the term involving the interior temperature T . First, note the following which hold on the unit disk r =|ζ|= 1[20, 22] Re f ζf t ζ 1 ∂T σ(θ) v = , n ˆ ·∇T = Re[n∇T ]= = , (14) |f | |f | ∂r |f | ζ ζ ζ iθ where ∇= ∂ − i∂ , ζ = e , f = ∂f /∂ ζ and n = n + in = ζf /|f | is the complex x y ζ x y ζ ζ representation of the normal vector in the z−plane [34]. The function σ(θ) is found by solving the following integral equation [22, 26, 29] π = K Pe(cos θ − cos θ ) σ(θ )dθ , (15) 0 7 Page 8of17 S.J.Harris, N.R. McDonald where K(x) = e K (|x|) and K is the modified Bessel function. The function σ is the heat 0 0 flux (due to the exterior effect alone) around the unit ζ −disk [22]. A numerical solution to (15) can be found at M equally spaced points on the upper half of the unit disk θ ∈ (0,π),see ∗ ∗ appendix A of [22]. To evaluate σ on the remainder of the disk, note that σ(θ ) = σ(π + θ ), where 0 ≤ θ <π,and σ(θ) → 1/π as θ → 0[22, 26]. Second, consider the contribution from the interior n ˆ ·∇T . This is found by direct com- putation in the z−plane. To begin, introduce the function u = T + αr /4; from (13d)and (13e), u satisfies ∇ u =0in R, (16) αr u = 1 + on γ. (17) This is a Laplace problem in the interior of a smooth, simply connected domain with Dirich- let boundary condition u = h(z) on γ . Efficient methods developed by Trefethen and col- leagues are here used to solve the Laplace problem directly, see [47, 55] for further details. These are based on a rational approximation of u in the complex z−plane: there exists an analytic function F(z) such that u = Re[F(z)] which can be approximated by the series solution N N 1 2 F(z) = A + A (z − c ) + , (18) 0 k 0 z − p k=1 k=1 where A , A and B are constants, N is the power series truncation, c is a point in the 0 k k 1 0 interior R (in this work, the conformal centre of the polygon is used [22]) and p are some N poles in the exterior . The first sum is known as the “Runge” (or smooth) part and the second as the “Newman” (or singular) part [49, 55]. Suitable poles p are found by the “adaptive Antoulas-Anderson” (AAA, pronounced ‘triple-A’) algorithm [49, 56] and the unknown coefficients A , A and B by a least-squares 0 k k (LS) algorithm [46, 48, 57]: combining these into a AAA-LS method [50] can thus find F(z) numerically. It follows [46]that ∇u =∇ Re[F(z)]= F (z), (19) where the primed notation indicates differentiation, in this case with respect to z. Therefore, ∇T = ∇[u − αr /4]= F (z) − z, (20) recalling that ∇= 2∂ and r = zz.Thisgives 1 ω(θ) n ˆ ·∇T = Re[n∇T ]= Re[ζf (F (z) − αz/2)]≡ , (21) R R ζ |f | |f | ζ ζ where the boundary points z in (21) are considered in the ζ −plane so that the right hand side of (21) is purely a function of θ . Drawing an analogue with σ , the function −βω is the heat flux around the unit ζ −disk due to the interior effect alone. Third, the formula (13g) for the area conserving term C(t) can be transformed to the ζ −plane, noting that ds =|f |dθ . Substituting (14)and (21)into(13g)gives − σ − βω dθ −π C(t) =  . (22) |f |dθ −π Penguin Huddling: A Continuum Model Page 9 of 17 7 Therefore, using (14), (21)and (22), the PG form of (13a)is Pe Re[f ζf ]= σ(θ) − βω(θ) + C(t)|f |, (23) t ζ ζ where σ(θ) satisfies (15)and ω(θ) = Re[ζf (F (z) − αz/2)], with F(z) given by (18). Note that σ(θ) − βω(θ) is the total heat flux across |ζ|= 1; recall that the heat flux across the huddle boundary γ is n ˆ ·∇T − βn ˆ ·∇T . It remains to find the conformal map f , the general form of which is [24] −k z = f(ζ, t) = a (t )ζ + c (t )ζ , (24) −1 k k=0 where a (t ) is a real function in time and c (t ) = a (t ) + ib (t ) are complex functions −1 k k k in time. Note that a is the conformal radius and c is the conformal centre [22]. The −1 0 same c is used as the interior point in constructing the power series solution (18)for the Poisson problem in R. Numerically, the Laurent series (24) is truncated at N terms, giving n = 2(N + 1) + 1 = 2N + 3 unknown real functions in time: a (t ), a (t ), b (t ), k = −1 k k 0, 1,...,N . Likewise, n equally spaced points around the unit ζ −disk are chosen; the first point is at θ = 0, the next N + 1 points are in the interval θ ∈ (0,π) and the final N + 1 points in θ ∈ (π, 2π). Hence the choice M = N + 1in (0,π) is made when numerically solving the integral equation (15)for σ(θ). The PG equation (23) becomes a system of n coupled ordinary differential equations (ODEs) determining the time evolution of the n Laurent coefficients (see [20, 21, 34]). In this work, the MATLAB routine ode15i was used to solve this ODE system. 4 Results In the results presented here, the evolution of the huddle boundary is plotted at equal time intervals in the range [0,t ],where t is the maximum time. The series truncation max max N = 128 is chosen, leading to a system of n = 259 ODEs to be solved. The initial start- ing shape of the penguin huddle is taken to be either a circle, a slanted ellipse, a triangle, an irregular pentagon or an hourglass. The initial Laurent coefficients a (0), c (0) for the −1 k conformal map (24) of these shapes are taken from [24] (see their Figures 2a, 2b, 2d, 5a and 7a, respectively) and scaled such that each shape encloses the same area. Area conserva- tion of the huddle was monitored during each numerical experiment and found to hold with −5 a relative error of O(10 ). As a test, experiments were run which neglected the interior (i.e. β = 0) and the area conserving (C(t)) terms; the geometry of the obtained solutions matched quantitatively with results in [22] who considered an analogous problem in this class. First, consider the “exterior only” penguin problem when β = 0. This is a continuum ana- log of the discrete model developed in [6], where the penguin huddle evolves in response to the wind but the interior reorganisation of individual penguins is not represented. Figure 2 shows results for various Péclet values and initial starting shapes of the huddle. In general, these show good agreement to the experiments of [6]: the huddle propagates in the windward direction and evolves into a more streamlined shape over time. For large time (t  t /2), max This scaling can be found analytically using the complex form of Green’s theorem [45] to find the exact area from the conformal map z = f(ζ, 0). 7 Page 10 of 17 S.J. Harris, N.R. McDonald Fig. 2 Exterior only (β = 0) penguin problem. Shapes are plotted at equal time intervals up to a maximum time t .The x and y axes (not shown) are scaled at a 1-to-1 aspect ratio. Each description refers to the max initial penguin huddle shape, the Péclet number and the value of t : a) circle with Pe = 1, t = 10, b) max max circle with Pe = 10, t = 20, c) slanted ellipse with Pe = 10, t = 20, d) irregular pentagon with Pe max max = 100, t = 50 max the huddle approaches an invariant egg-shape which continues to propagate downwind with some constant velocity. To test this observation, a root mean-squared error (RMSE) measur- ing the deviation between consecutive plotted huddle shapes is computed for the examples −4 in Fig. 2 and found to approach a constant value with a relative error of O(10 ).Thisper- mits the reasonable conclusion that huddles do evolve into a steady shape and thus in all experiments presented, the maximum time t is chosen such that the penguin huddle has max reached a steady shape prior to or at time t . max Now, consider the full penguin problem where exterior, interior (β = 0) and area con- serving terms are all included - recall from Sect. 2.1 that β = 0.1 is used. Figures 3a-e present the evolution of the same initial starting shape, a circle, with Pe = 1, 5, 10, 50 and 100. The final plot in each experiment is the steady shape and these are compared in Fig. 3f, where the shapes have been re-centered so that they share the same centre of mass. As the Péclet number is increased, the steady shape becomes less circular and more egg like. Similarly, the effect of changing the parameter β can be explored. Figures 4a-e shows the evolution of the circle with the values β = 0, 0.1, 0.5, 1 and 1.5, respectively, where Pe = 10 and t = 20 are taken for all plots. While a low value of β gives results that look max similar to the exterior only penguin problem, a high value of β causes the penguin huddle to retain a more circular shape. This is shown in Fig. 4f, where the steady shapes of Figs. 4a-e are compared. The dependence of the final steady shape of the huddle on the values of the Péclet num- ber and β is explored further. A given steady shape can be quantified by its aspect ratio (AR): the ratio between its major and minor axes with AR = 1 being a circle and increasing AR > 1 becoming more egg-shaped. Figures 5a,b show the dependence of AR on Pe and β , respectively. Each marker on the diagram represents an experiment run: each curve contains 21 markers, all run to a maximum time of t = 80, which is sufficiently long to reach the max steady state huddle shape. The results from Fig. 5 agree with the conclusions drawn from Figs. 3 and 4.First,asPe increases, AR increases and thus a higher Péclet number results in a more egg-like steady shape, with a lower value of β giving higher values of AR for all time. This is also consistent with the parameter study in [6] - see their Fig. 3 - which shows that the huddle thickness decreases with an increasing Péclet number. As Pe → 0, AR→ 1, consistent with the limit of a circular steady shape when there is no wind, as expected. Second, as β increases, AR decreases, where lower values of Pe result in lower values of the AR for all time. Penguin Huddling: A Continuum Model Page 11 of 17 7 Fig. 3 The effect of different Péclet numbers on the evolution of the huddle is shown by plotting its shape at equal time intervals. The initial starting shape of the penguin huddle is a circle, with β = 0.1and a) Pe = 1, t = 10, b) Pe = 5, t = 15, c) Pe = 10, t = 20, d) Pe = 50, t = 40, e) Pe = 100, t = 50. f) max max max max max Comparison of the steady shapes of (a)-(e), the initial (circular) huddle shape is also plotted for comparison Fig. 4 Interior β effect: initial starting shape of a circle with Pe = 10, t = 20 and a) β = 0, b) β = 0.1, max c) β = 0.5, d) β = 1, e) β = 1.5. f) Comparison of the steady shapes of (a)-(e), the initial (circular) huddle shape is also plotted for comparison 7 Page 12 of 17 S.J. Harris, N.R. McDonald Fig. 5 The aspect ratio of the final steady shape of the penguin huddle as a function of a) Péclet number, for the three given values β = 0, 0.1and 1, and b) β , for the three given values Pe = 10, 50 and 100. Each marker represents a numerical experiment Fig. 6 Initial shape effect. Shapes are plotted at equal time intervals with Pe = 10, β = 0.1, t = 35 and the max initial starting shape of a) a circle, b) a slanted ellipse, c) a triangle, d) an irregular pentagon, e) an hourglass. f) Comparison of the steady shapes of (a)-(e): all of the steady shapes are overlapped Figure 6 shows how huddles of different starting shapes evolve under the same Péclet number Pe = 10, the same interior effect β = 0.1 and up to the same maximum time t = max 35. Figures 6a-e show the evolution of a circle, a slanted ellipse, a triangle, an irregular pentagon and an hourglass, respectively, and Fig. 6f compares the steady shapes of these five experiments. The results show that, irrespective of the starting shape of the huddle and even those starting with significant concave boundaries, the same steady shape is reached. The along-boundary variation of the total heat flux across the huddle boundary of the steady shape of Fig. 6f is considered. By symmetry, only the upper half of the boundary need be considered; this is plotted in Fig. 7a. The arrows point in the normal direction of the Penguin Huddling: A Continuum Model Page 13 of 17 7 Fig. 7 a) Total heat flux across the boundary of the steady shape of Fig. 6f. The arrows point in the normal direction of the given boundary point and their magnitude is proportional to the local heat flux across the boundary. b) Canonical heat flux in the upper half ζ −disk vs θ . c) Physical heat flux in the upper half of the huddle boundary vs x given point on the boundary and their lengths represent the total physical heat flux at that point. Recall from Sect. 3 that the total heat flux in the physical domain across the huddle boundary γ is given by n ˆ ·∇T − βn ˆ ·∇T whereas the total heat flux in the canonical domain across the unit ζ −disk is given by σ(θ) − βω(θ). The canonical heat flux in the upper half ζ −disk vs θ is plotted in Fig. 7b and the physical heat flux vs x is plotted in Fig. 7c. There is a much larger heat flux on the windward side of the penguin huddle (the side of the penguin huddle directly facing the oncoming wind) compared to the leeward side. This is expected and is in agreement with the exterior temperature profile presented in [6] - see their Fig. 1 - which shows that there is a steeper temperature gradient, and hence a higher total heat flux, on the windward side. 5 Discussion In this continuum model, the evolution of the penguin huddle boundary γ(t), is affected by three factors: exterior cooling by the ambient wind, interior reorganisation of the penguins driven by their self-generation of heat, and area conservation of the huddle. The resulting two-dimensional model (13a)-(13g) can be rewritten as a PG-type equation (23)interms of the conformal map z = f(ζ, t) between the canonical (ζ ) and physical (z) planes. The map is approximated as a truncated Laurent series (24), whose unknown time-varying coefficients are found by numerically solving a system of coupled ODEs. Results from Sect. 4 show that the huddle shape eventually reaches a steady shape. As the wind strength increases, the huddle becomes increasingly egg-shaped, minimising the portion of the huddle bound- ary directly facing the oncoming wind. For a stronger interior heating effect, the huddle 7 Page 14 of 17 S.J. Harris, N.R. McDonald approaches a more circle-like steady shape, to optimise interior heat regulation. For given wind and interior effects, any initial huddle shape eventually approaches the same steady shape. There is also a greater total heat flux on the windward side of the huddle than on the leeward side. While dimensionless quantities have been used throughout the paper, it is useful to esti- mate the dimensional time scale associated with the huddle evolution. As noted in Sect. 2.1, ∗ ∗ the appropriate dimensional time scale is t = λPe(L/U )t where t is dimensionless time, and the choice λ = 1 has been made. A large huddle of, for example, 6000 penguins [1], is assumed and taking the radius of each penguin to be r ≈ 0.1m [3]gives L ≈ 10 m. Now −1 assume a sufficiently strong wind of U ≈ 20 ms for Pe = 100. Considering Fig. 3cwhere t = t = 50, the total real time for this experiment is t = 2500 seconds or ≈ 42 min- max r utes. Thus it is concluded that the penguin huddle motion described in this work occurs over roughly a one hour period, in agreement with previous modelling and field observations [6, 8, 9], and justifying the choice λ = 1. In this paper, a continuum model has been used to find the shape and propagation of the huddle boundary at the expense of tracking individual penguin movements. While there are pre-existing discrete models for penguin huddling (and other animal grouping dynam- ics) [4, 6, 7, 36], treating a penguin huddle as a continuum appears to be new. One of the advantages of the continuum model developed here is its computational efficiency, with a typical runtime of around 5-10 minutes on a standard laptop for a huddle to reach a steady shape. This can be sped up further, to a matter of seconds, by reducing the Laurent series truncation N , though care should be taken to ensure that the accuracy of the results is not compromised. When considering the “exterior only” (β = 0) penguin problem, our results matched well with those from [6], who did not include the interior effect of penguin reorganisation in their model. As discussed in Sect. 2.1, fully turbulent wind flows have β ≈ 0 and hence are well modelled by considering exterior effects alone. The representative value β = 0.1, which is likely to be an overestimate, chosen in this work demonstrates that the interior heating of the colony only provides a weak effect in determining the steady shape of the huddle, see Fig. 4f. Larger values of β resulting in more circular huddle shapes have been included here for completeness. Further, the results in this work match well with real world observations: penguin huddles are seen to march downwind over time and observed huddle shapes resemble the egg-like steady shapes found here - see [1, 4, 10, 12, 13]. Extensions to the penguin problem are of interest. A time-varying wind could be consid- ered, where the variations are on the order of an hour; when the wind direction changes, the huddle will seek to propagate back to the optimal (egg) steady shape in the new windward direction. Stochastic terms could also be incorporated into the model, as in [6]. This includes random variations in the direction and magnitude of the wind and stochasticity in the interior and exterior heat fluxes. The numerical method would need to be adapted to handle these stochastic terms. The evolution of any arbitrary starting huddle shape could also be simulated. In this work, the initial huddle boundary was taken such that its conformal map was readily represented in the form of a truncated Laurent series (24). For more arbitrary starting shapes, the hud- dle boundary can be discretised as a polygon of m vertices, then its conformal map to the ζ −plane found numerically at each time step as the penguin huddle evolves. Pre-exisiting numerical schemes for finding the conformal map between a polygon and the unit ζ −disk include the Schwarz-Christoffel toolbox [58] and a variation of the AAA-LS algorithm [50] from Sect. 3. It would also be possible to consider multi-connected huddles such as those with penguin-free holes or the interaction, and possible merger, of initially disjoint huddles. Penguin Huddling: A Continuum Model Page 15 of 17 7 While this would require a significant modification of the numerical method used here, it is worth remarking that AAA-LS type methods are able to handle such geometries [47, 49, 50]. Finally, consider the boundary condition (13a), and its associated PG equation (23). The right-hand side is constructed by exterior (n ˆ ·∇T ), interior (βn ˆ ·∇T ) and area conserva- tion (C(t)) effects. By inclusion and exclusion of these three terms, (13a) encapsulates six possible free boundary problems: 1. The dissolution problem: exterior effect only. A permeable object is placed in a uniform stream of some dissolving agent - see [22, 59]. 2. The Poisson growth problem: interior effect only. Applicable to squeeze flow and the evaporation of thin liquid films - see [45, 60, 61]. 3. The two-phase melting/freezing problem: exterior and interior effects. Porous media flow about freeze pipes where interior and exterior temperature gradients are in effect - see [24, 29]. 4. The “exterior only” penguin problem: exterior and area conservation effects. Huddle problems where interior reorganisation is not represented, such as if individuals in the centre are unable to move or there exists a “pecking order” in the population - see [6]. 5. The “bat huddle” problem: interior and area conservation effects. Huddle problems with no wind effect, for example bats huddling in caves sheltered from exterior winds. To the authors’ knowledge, there is no mathematical work on this problem, but for biological context see [39, 40]. 6. The “full” penguin problem: exterior, interior and area conservation effects. This is the full huddle continuum problem covered in this work. Therefore, the numerical method developed in this work can be used to model and simulate additional free boundary problems, such as those from the fields of fluid mechanics and mathematical biology. Acknowledgements The authors thank Carole Hall, Stony Brook University, for their inspiring talk on Adélie penguins and the useful discussion thereafter. The authors also thank the anonymous reviewers for their useful comments and suggestions. Author Contribution Both authors contributed equally to this work. Funding Samuel J. Harris was supported by a UK Engineering and Physical Sciences Research Council PhD studentship, grant numbers EP/N509577/1 and EP/T517793/1. Data Availability All data generated or analysed during this study are included in this published article and its supplementary information files. Code Availability All relevant code files are publicly available at the following GitHub link: https://github. com/Sam-J-Harris/Penguin-huddling-a-continuum-model.git. Declarations Competing Interests Not applicable. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 7 Page 16 of 17 S.J. Harris, N.R. McDonald References 1. Le Maho, Y.: The emperor penguin: a strategy to live and breed in the cold: morphology, physiology, ecology, and behavior distinguish the polar emperor penguin from other penguin species, particularly from its close relative, the king penguin. Am. Sci. 65(6), 680–693 (1977) 2. McCafferty, D.J., Gilbert, C., Thierry, A.-M., Currie, J., Le Maho, Y., Ancel, A.: Emperor penguin body surfaces cool below air temperature. Biol. Lett. 9(3), 20121192 (2013) 3. Williams, C.L., Hagelin, J.C., Kooyman, G.L.: Hidden keys to survival: the type, density, pattern and functional role of emperor penguin body feathers. Proc. R. Soc. Lond. B 282(1817), 20152033 (2015) 4. Gerum, R.C., Fabry, B., Metzner, C., Beaulieu, M., Ancel, A., Zitterbart, D.P.: The origin of traveling waves in an emperor penguin huddle. New J. Phys. 15(12), 125022 (2013) 5. Kooyman, G.L., Gentry, R.L., Bergman, W.P., Hammel, H.T.: Heat loss in penguins during immersion and compression. Comp. Biochem. Physiol. 54(1), 75–80 (1976) 6. Waters, A., Blanchette, F., Kim, A.D.: Modeling huddling penguins. PLoS ONE 7(11), 50277 (2012) 7. Gerum, R., Richter, S., Fabry, B., Le Bohec, C., Bonadonna, F., Nesterova, A., Zitterbart, D.P.: Structural organisation and dynamics in king penguin colonies. J. Phys. D, Appl. Phys. 51(16), 164004 (2018) 8. Gilbert, C., Robertson, G., Le Maho, Y., Naito, Y., Ancel, A.: Huddling behavior in emperor penguins: dynamics of huddling. Physiol. Behav. 88(4–5), 479–488 (2006) 9. Kirkwood, R., Robertson, G.: The occurrence and purpose of huddling by emperor penguins during foraging trips. Emu 99(1), 40–45 (1999) 10. Ancel, A., Gilbert, C., Poulin, N., Beaulieu, M., Thierry, B.: New insights into the huddling dynamics of emperor penguins. Anim. Behav. 110, 91–98 (2015) 11. Zitterbart, D.P., Wienecke, B., Butler, J.P., Fabry, B.: Coordinated movements prevent jamming in an emperor penguin huddle. PLoS ONE 6(6), 20260 (2011) 12. Mina, T., Min, B.-C.: Penguin huddling inspired distributed boundary movement for group survival in multi-robot systems using Gaussian processes. In: 2018 IEEE International Conference on Robotics and Biomimetics, pp. 2177–2183. IEEE (2018) 13. Richter, S., Gerum, R., Winterl, A., Houstin, A., Seifert, M., Peschel, J., Fabry, B., Le Bohec, C., Zit- terbart, D.P.: Phase transitions in huddling emperor penguins. J. Phys. D, Appl. Phys. 51(21), 214002 (2018) 14. Gu, W., Christian, J.K., Woodson, C.B.: A novel coupled fluid-behavior model for simulating dynamic huddle formation. PLoS ONE 13(8), 0203231 (2018) 15. Saffman, P.G., Taylor, G.I.: The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc. R. Soc. A 245(1242), 312–329 (1958) 16. Howison, S.: Fingering in Hele-Shaw cells. J. Fluid Mech. 167, 439–453 (1986) 17. Paterson, L.: Radial fingering in a Hele Shaw cell. J. Fluid Mech. 113, 513–529 (1981) 18. Taylor, G., Saffman, P.G.: A note on the motion of bubbles in a Hele-Shaw cell and porous medium. Q. J. Mech. Appl. Math. 12(3), 265–279 (1959) 19. Entov, V.M., Etingof, P.I.: Bubble contraction in Hele-Shaw cells. Q. J. Mech. Appl. Math. 44(4), 507–535 (1991) 20. Dallaston, M.C., McCue, S.W.: An accurate numerical scheme for the contraction of a bubble in a Hele– Shaw cell. ANZIAM J. 54, 309–326 (2013) 21. Dallaston, M.C., McCue, S.W.: A curve shortening flow rule for closed embedded plane curves with a prescribed rate of change in enclosed area. Proc. R. Soc. A 472(2185), 20150629 (2016) 22. Ladd, A.J.C., Yu, L., Szymczak, P.: Dissolution of a cylindrical disk in Hele-Shaw flow: a conformal- mapping approach. J. Fluid Mech. 903, 46 (2020) 23. Cummings, L.M., Hohlov, Y.E., Howison, S.D., Kornev, K.: Two-dimensional solidification and melting in potential flows. J. Fluid Mech. 378, 1–18 (1999) 24. Rycroft, C.H., Bazant, M.Z.: Asymmetric collapse by dissolution or melting in a uniform flow. Proc. R. Soc. A 472(2185), 20150531 (2016) 25. Mullins, W.W., Sekerka, R.F.: Stability of a planar interface during solidification of a dilute binary alloy. J. Appl. Phys. 35(2), 444–451 (1964) 26. Choi, J., Margetis, D., Squires, T.M., Bazant, M.Z.: Steady advection–diffusion around finite absorbers in two-dimensional potential flows. J. Fluid Mech. 536, 155–184 (2005) 27. Tsai, V.C., Wettlaufer, J.S.: Star patterns on lake ice. Phys. Rev. E 75(6), 066105 (2007) 28. Grodzki, P., Szymczak, P.: Reactive-infiltration instability in radial geometry: from dissolution fingers to star patterns. Phys. Rev. E 100(3), 033108 (2019) 29. Goldstein, M.E., Reid, R.L.: Effect of fluid flow on freezing and thawing of saturated porous media. Proc. R. Soc. A 364(1716), 45–73 (1978) 30. Langer, J.S.: Instabilities and pattern formation in crystal growth. Rev. Mod. Phys. 52(1), 1 (1980) Penguin Huddling: A Continuum Model Page 17 of 17 7 31. Brower, R.C., Kessler, D.A., Koplik, J., Levine, H.: Geometrical models of interface evolution. Phys. Rev. A 29(3), 1335 (1984) 32. Sethian, J.A.: Curvature and the evolution of fronts. Commun. Math. Phys. 101(4), 487–499 (1985) 33. Hilton, J.E., Miller, C., Sharples, J.J., Sullivan, A.L.: Curvature effects in the dynamic propagation of wildfires. Int. J. Wildland Fire 25(12), 1238–1251 (2016) 34. Harris, S.J., McDonald, N.R.: Fingering instability in wildfire fronts. J. Fluid Mech. 943, 34 (2022) 35. Burger, M., Haškovec, J., Wolfram, M.-T.: Individual based and mean-field modeling of direct aggrega- tion. Physica D 260, 145–158 (2013) 36. Bernardi, S., Scianna, M.: An agent-based approach for modelling collective dynamics in animal groups distinguishing individual speed and orientation. Philos. Trans. R. Soc. Lond. B 375(1807), 20190383 (2020) 37. Katz, Y., Tunstrøm, K., Ioannou, C.C., Huepe, C., Couzin, I.D.: Inferring the structure and dynamics of interactions in schooling fish. Proc. Natl. Acad. Sci. 108(46), 18720–18725 (2011) 38. Bhattacharya, K., Vicsek, T.: Collective decision making in cohesive flocks. New J. Phys. 12(9), 093019 (2010) 39. Herreid, C.F.: Temperature regulation of Mexican free-tailed bats in cave habitats. J. Mammal. 44(4), 560–573 (1963) 40. Ryan, C.C., Burns, L.E., Broders, H.G.: Changes in underground roosting patterns to optimize energy conservation in hibernating bats. Can. J. Zool. 97(11), 1064–1070 (2019) 41. Nave, G.K. Jr, Mitchell, N.T., Chan Dick, J.A., Schuessler, T., Lagarrigue, J.A., Peleg, O.: Attraction, dynamics, and phase transitions in fire ant tower-building. Front. Robot. AI 7, 25 (2020) 42. Ko, T.-Y.H., Yu, H.D.L.: Fire ant rafts elongate under fluid flows. Bioinspir. Biomim. 17(4), 045007 (2022) 43. Bazant, M.Z., Crowdy, D.: Conformal mapping methods for interfacial dynamics. In: Handbook of Ma- terials Modeling, pp. 1417–1451. Springer, Dordrecht (2005) 44. Gustafsson, B., Vasil’ev, A.: Conformal and Potential Analysis in Hele-Shaw Cells. Springer, Berlin (2006) 45. McDonald, R., Mineev-Weinstein, M.: Poisson growth. Anal. Math. Phys. 5(2), 193–205 (2015) 46. Trefethen, L.N.: Series solution of Laplace problems. ANZIAM J. 60(1), 1–26 (2018) 47. Trefethen, L.N.: Numerical conformal mapping with rational functions. Comput. Methods Funct. Theory 20(3), 369–387 (2020) 48. Baddoo, P.J.: Lightning solvers for potential flows. Fluids 5(4), 227 (2020) 49. Costa, S.: Solving Laplace problems with the AAA algorithm (2020). arXiv:2001.09439. ArXiv preprint 50. Costa, S., Trefethen, L.N.: AAA-least squares rational approximation and solution of Laplace problems (2021). arXiv:2107.01574. ArXiv preprint 51. Gupta, S.C.: The Classical Stefan Problem: Basic Concepts, Modelling and Analysis with Quasi- Analytical Solutions and Methods, vol. 45. Elsevier, Amsterdam (2017) 52. Basu, P.: Biomass Gasification, Pyrolysis and Torrefaction: Practical Design and Theory. Academic Press, Boston (2018) 53. Dawson, C., Vincent, J.F.V., Jeronimidis, G., Rice, G., Forshaw, P.: Heat transfer through penguin feath- ers. J. Theor. Biol. 199(3), 291–295 (1999) 54. Kittel, C., McEuen, P., McEuen, P.: Introduction to Solid State Physics, vol. 8. Wiley, New York (1996) 55. Gopal, A., Trefethen, L.N.: Solving Laplace problems with corner singularities via rational functions. SIAM J. Numer. Anal. 57(5), 2074–2094 (2019) 56. Nakatsukasa, Y., Sète, O., Trefethen, L.N.: The AAA algorithm for rational approximation. SIAM J. Sci. Comput. 40(3), 1494–1522 (2018) 57. Brubeck, P.D., Nakatsukasa, Y., Trefethen, L.N.: Vandermonde with Arnoldi. SIAM Rev. 63(2), 405–415 (2021) 58. Driscoll, T.A.: Algorithm 756: a MATLAB toolbox for Schwarz-Christoffel mapping. ACM Trans. Math. Softw. 22(2), 168–186 (1996) 59. Dutka, F., Starchenko, V., Osselin, F., Magni, S., Szymczak, P., Ladd, A.J.: Time-dependent shapes of a dissolving mineral grain: comparisons of simulations with microfluidic experiments. Chem. Geol. 540, 119459 (2020) 60. Agam, O.: Viscous fingering in volatile thin films. Phys. Rev. E 79(2), 021603 (2009) 61. Crowdy, D., Kang, H.: Squeeze flow of multiply-connected fluid domains in a Hele-Shaw cell. J. Non- linear Sci. 11(4), 279–304 (2001) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Acta Applicandae Mathematicae Springer Journals

Penguin Huddling: A Continuum Model

Download PDF
17 pages

Loading next page...
 
/lp/springer-journals/penguin-huddling-a-continuum-model-xCHuzCm2Ya

References (70)

Publisher
Springer Journals
Copyright
Copyright © The Author(s) 2023
ISSN
0167-8019
eISSN
1572-9036
DOI
10.1007/s10440-023-00578-2
Publisher site
See Article on Publisher Site

Abstract

Penguins huddling in a cold wind are represented by a two-dimensional, continuum model. The huddle boundary evolves due to heat loss to the huddle exterior and through the re- organisation of penguins as they seek to regulate their heat production within the huddle. These two heat transfer mechanisms, along with area, or penguin number, conservation, gives a free boundary problem whose dynamics depend on both the dynamics interior and exterior to the huddle. Assuming the huddle shape evolves slowly compared to the advective timescale of the exterior wind, the interior temperature is governed by a Poisson equation and the exterior temperature by the steady advection-diffusion equation. The exterior, ad- vective wind velocity is the gradient of a harmonic, scalar field. The conformal invariance of the exterior governing equations is used to convert the system to a Polubarinova-Galin type equation, with forcing depending on both the interior and exterior temperature gradi- ents at the huddle boundary. The interior Poisson equation is not conformally invariant, so the interior temperature gradient is found numerically using a combined adaptive Antoulas- Anderson and least squares algorithm. The results show that, irrespective of the starting shape, penguin huddles evolve into an egg-like steady shape. This shape is dependent on the wind strength, parameterised by the Péclet number Pe, and a parameter β which effectively measures the strength of the interior self-generation of heat by the penguins. The numerical method developed is applicable to a further five free boundary problems. Keywords Free boundary problems · Conformal mapping · AAA-least squares algorithm · Fluid mechanics · Collective behaviour · Huddling 1Introduction Among all the penguin species living in the Antarctic region [1], the emperor penguin (Aptenodytes forsteri) experiences the most extreme weather conditions, enduring wind −1 ◦ speeds over 30 ms at temperatures below −40 C[1–4]. Unlike other species such as the gentoo and the Adélie [5], emperor penguins have no fixed nest [6, 7], allowing the birds to huddle together for warmth in severe wind conditions. This strategy is very effective; ambient temperatures within the huddle often exceed 20 C, with a maximum temperature of 37.5 C having been recorded [8]. Such huddling is crucial for the penguins’ survival, S.J. Harris sam.harris.16@ucl.ac.uk Department of Mathematics, University College London, Gower Street, London, WC1E 6BT, UK 7 Page 2of17 S.J.Harris, N.R. McDonald especially during foraging trips [9] and throughout their breeding season, which coincides with the Antarctic winter [3]. Over the few hours that huddling events typically last [6, 8], the penguins are contin- uously (albeit slowly) on the move. The birds constantly reorganise their position in the huddle, where [1] observes that “birds that at first are in the center of the huddle become members of the rear flank [the side most exposed to the wind] and move, in their turn, up the sideline”. This process benefits all members of the huddle and allows the heat to be equally shared: those at the edge do not remain cold and those at the centre do not overheat. These individual penguin movements affect the propagation of the huddle as a whole - the huddle is seen to propagate downwind [10] and the shape traced out by the edge of the huddle, or huddle boundary, may change over time [4, 6]. It is the propagation and evolution of the huddle boundary shape that is the focus of this work. The modelling of penguin huddle dynamics has received increased attention over the last decade, but still remains sparse in the literature, despite numerous observations and findings from the field, see [1–3, 10]. The works [4, 7, 11] draw the analogy between penguin hud- dle/colony formation and condensed matter physics, with Lennard-Jones-like forces used to represent attractive and repulsive interactions between penguins that allow complex, lattice structures (such as huddles) to form. This idea is expanded upon in [12], where Gaussian processes are used to test if robots could adopt this penguin inspired huddling strategy. Sim- ilarly, [13] liken the rearrangement of the penguin huddle to a phase transition from solid (dense huddle) to liquid (loose huddle) to gas (no huddle) depending on environmental fac- tors such as temperature and wind speed. A rigorous fluid dynamics approach is taken in [14], where a finite difference method is used to solve the full Navier-Stokes equations of the surrounding fluid (the exterior wind), giving a two-way coupling between the environ- ment and the behaviour of the huddle. Using a simpler fluids based model, [6] assume the exterior wind to be a two-dimensional, irrotational flow of an incompressible fluid and from this model, the temperature profile exterior to the huddle and individual penguin movements can be found. In this work, the approach of [6] is expanded upon with two key differences. First, a simple model for the effect of interior penguin reorganisation as a result of self-generation of heat throughout the entire huddle, as observed in [1], is included. In [6], while penguin reorganisation around the huddle boundary was acknowledged, it was assumed that interior penguins were so tightly packed that they could not change their position in the huddle. That assumption is here dropped, yet it is still assumed that the penguins are packed tightly enough such that wind does not permeate through the huddle and so the wind flow remains entirely in the exterior region. Throughout this work, the phrases “exterior only penguin problem” and “full penguin problem” are used to refer to the penguin huddle problem with- out and with interior penguin reorganisation, respectively. Second, the penguin huddle itself is treated as a continuum. Existing huddle models treat penguins as discrete elements (or particles) such that individual penguin movements can easily be tracked. For sufficiently large huddles on the order of a thousand penguins [1], these discrete models quickly become computationally expensive. Continuum models, as used here, are often faster, cheaper, and preferred when studying the change in shape and position of the overall huddle. Using a continuum to model the interior and exterior of the full penguin problem draws analogues to other two-dimensional free boundary problems. Foremost among these is the classical Hele-Shaw free-boundary problem and its extensions, including Saffman-Taylor fingering (a consequence of injecting a fluid into another, more viscous fluid) [15–17], the expansion and contraction of bubbles in a Hele-Shaw cell [18–21] and the dissolution of soluble objects in two-dimensional potential flow [22]. Using conformal mapping to re- formulate the free boundary problem as a Polubarinova-Galin (PG) equation has proved Penguin Huddling: A Continuum Model Page 3 of 17 7 effective in finding exact and numerical solutions for this class of problems. Solidification and melting problems have also been extensively studied using analogous techniques, see [23–28]and [29]; the latter of these examines the freezing of porous media flow by freeze pipes, with the velocity of the free boundary dependent on both interior and exterior temper- ature gradients, similar to the full penguin problem. Free boundary problems arise in other natural processes, including, but not limited to, dendritic crystal growth [30, 31]and flame front propagation [32–34]. Continua may also be used to model other collective behaviours [4, 35, 36] such as fish schooling [37], bird flocking [38], bat swarming [39, 40]and ant colony formation [41, 42]. There are two heat transfer mechanisms in the full penguin problem, which formulate a suitable boundary condition for the propagation of the huddle boundary. The first is the cooling effect of the wind on the huddle. As in [6], the incompressible wind flow is assumed irrotational and flowing sufficiently faster than the propagation of the huddle, thus, the flow can be treated as steady. As such, the wind flow and temperature transport in the region exterior to the penguin huddle can be modelled by the Laplace and the steady advection- diffusion equations, respectively. This coupled pair of partial differential equations (PDEs) is conformally invariant [43] - a property which motivates the use of a conformal mapping approach in the numerical method [21–24, 34]. The exterior region and the huddle bound- ary in the physical z−plane can be conformally mapped to the exterior and boundary of the unit disk in some “mathematical” (or “canonical”) ζ −plane, and vice versa. Therefore, the physical coordinates can be expressed in terms of the conformal map z = f(ζ, t) and the condition on the huddle boundary rewritten as a PG type equation [44]. The task then becomes that of finding the conformal map. The second heat transfer mechanism is that of the uniform, self-generation of heat by the penguins in the interior. Its steady diffusion amongst the densely packed penguins is modelled by a Poisson equation for the interior temperature with a constant forcing term. By the Riemann mapping theorem (see e.g. [43]), the interior of the penguin huddle will not also be mapped to the interior of the ζ −disk by the conformal map z = f(ζ, t) used to map the exteriors. Further, the Poisson equation is not conformally invariant [45]. Hence, to solve for temperature in the interior using a conformal mapping approach, a different conformal map f would need to be found numerically and the transformed Poisson equation in the ζ − plane solved. Alternatively, the interior temperature can be found directly in the z−plane without the use of conformal mapping. Here, the methods developed by Trefethen and colleagues [46–50] are used to find a rational approximation of a harmonic function in some simply connected domain as the real part of some analytic function, subject to given boundary conditions. This function consists of a power series in z and a series of rational functions involving poles in the exterior clustered near singularity points [47]. Both the adaptive Antoulas-Anderson algorithm (AAA, pronounced triple-A) [49] and a least- squares (LS) algorithm, such as the lightning Laplace solver [48] can be used to numerically find this analytic function [50] and hence the required harmonic function. In this work, by first finding a particular solution, the problem of solving the interior Poisson equation then becomes one of solving the Laplace equation with a modified boundary condition. In turn, this is readily solved numerically by the AAA-LS method. The structure of this paper is as follows. The two-dimensional, continuum penguin hud- dle model is formulated in Sect. 2, with assumptions stated and dimensionless quantities introduced. A suitable boundary condition for the normal boundary velocity v is derived, incorporating effects of exterior wind cooling, internal heat generation by the penguins and area (penguin) conservation. The numerical method is then developed in Sect. 3,where the equation for v is rewritten as a PG type equation using a conformal mapping approach. n 7 Page 4of17 S.J.Harris, N.R. McDonald Attention here is given to the exterior heat flux σ in the ζ −plane and the interior heat flux which is found in the physical plane. Results are presented in Sect. 4, illustrating the evo- lution of penguin huddles of different starting shapes under a variety of wind and interior heating effect strengths. Conclusions are discussed in Sect. 5 and extensions to the work, and applications to further (free boundary) problems, noted. 2ModelSetup Consider a huddle of penguins with horizontal length scale L where the penguins are as- sumed homogeneous in size and shape. The huddle is considered as a continuum with pen- guins packed as tightly as possible while still allowing for huddle reorganisation [11]and hence uniform spacing between each penguin can be assumed constant across the entire huddle and for all time. Emperor penguins have a typical height of H = 1m [1]and their tight packing means that background winds are diverted around the huddle, allowing the dy- namics to be considered as two-dimensional, see Fig. 1. The boundary of the penguin huddle γ(t) is defined as the two-dimensional, non-overlapping curve enclosing the huddle; the fi- nite, simply connected interior of the curve occupied by the huddling penguins is labelled R and the penguin-free exterior as . The boundary curve γ is orientated anti-clockwise such that the unit normal n ˆ points from the interior to the exterior. The huddle R is buffeted by a unidirectional and constant strength wind of magnitude U which, without loss of generality, flows in the positive x ˆ direction. The wind is the irrota- tional flow of air, which is considered to be an inviscid and incompressible fluid which flows sufficiently rapidly compared to the evolution of the penguin huddle as a whole, so that the transport of temperature can be treated as steady [6]. In reality, viscous effects taking the form of thin boundary layers around the exterior of the penguin huddle are present. In this work, these are ignored in order to facilitate the accurate computation of the exterior po- tential flow using conformal mapping based methods. Denote the temperature far from the huddle as T , the temperature on the huddle boundary as T and the maximum temperature ∞ f in the huddle interior as T . The equations governing the wind and exterior temperature transport are, respectively, the Laplace equation and the steady advection-diffusion equation ∇ φ =0in , (1) u ·∇T = D ∇ T in , (2) where φ is the wind velocity potential, u = ∇φ is the wind velocity, T is the temperature in  and D is the diffusivity of temperature in air. The associated boundary and far field Fig. 1 Schematic diagram of the full penguin problem showing the interior region R and exterior region . a) Plan view where T denotes the maximum temperature T inside the huddle, b) side view c Penguin Huddling: A Continuum Model Page 5 of 17 7 conditions are ∂φ T = T , =0on γ, (3) ∂n T → T ,φ → Ux as r →∞. (4) Now, consider the interior R of the huddle where there is no wind. Owing to their ability to self-generate heat, which subsequently diffuses within the huddle, penguins nearer the huddle centre are warmer than those at the boundary [1, 11], i.e. T >T . This source of c f heat provides a further contribution to the heat flux at the boundary γ . Consistent with the continuum assumption, it is assumed that heat is generated uniformly across the huddle. Further, as in the exterior, it is assumed the interior heat diffuses quickly compared to the evolution of the boundary. Thus the temperature distribution in R is modelled by the steady heat equation ∇ T =− in R, (5) where D is the diffusivity of the interior temperature T through the penguin huddle and R R Q is the (constant) source term representing the self-generation of heat by the penguins. In order to obtain parameter estimates for the interior problem, it is helpful to first con- sider a circular penguin huddle of radius L, noting that T (0) = T . The Poisson equation R c (5) can then be solved exactly, giving Q Q 4(T − T ) c f T = T(r) = T − r , where = . (6) R c 4D D L R R The relation (6)isusedinSect. 2.1 to non-dimensionalise the problem for cases when the region R is non-circular. Finally, a condition on γ(t) is sought by specifying its outward normal velocity v .Fol- lowing [6], the idea that the boundary evolves owing to the heat loss of penguins on the boundary is used. Consequently, the problem has some similarity to the two-phase Stefan problem [23, 24, 29, 51] with the boundary temperature gradient driving its evolution, yet, unlike such problems, no phase change occurs. Waters et al. [6] move individual boundary penguins according to their heat loss to the exterior, whereas here both exterior and interior heat fluxes are considered. First, the heat loss to the exterior wind is the flux κ n ˆ ·∇T eval- uated on γ(t),where κ is the thermal conductivity of air. Second, the heat gain due to the warm huddle interior is −κ n ˆ ·∇T evaluated on γ(t),where κ is the thermal conductivity R R R of huddling penguins. Equating the normal velocity with the net heat flux gives Av = κ n ˆ ·∇T − κ n ˆ ·∇T , on γ(t), (7) n   R R where the constant term A is chosen to ensure a match in dimensions between both sides of (7), with the particular choice of A affecting the timescale for the boundary evolu- tion. For convenience in the non-dimensionalisation process upcoming in Sect. 2.1,let A = λρ c (T − T ),where ρ and c are the density and specific heat capacity of the f ∞ air, respectively, and λ is a non-dimensional constant to be chosen later. The right hand side of (7) is the jump in heat flux at the boundary. This observation reinforces the analogy with the Stefan problem where typically the normal velocity of the boundary depends on the jump in the heat flux [51]. For example, analogous boundary con- ditions arise in other natural free boundary problems, including the dissolution of finite 7 Page 6of17 S.J.Harris, N.R. McDonald objects, where only the gradient of the agent undergoing advection and diffusion exterior to the object is important [22, 24], or in porous media flow about freeze pipes, in which the temperature gradient in both the exterior fluid and interior ice regions are significant [29]. However it is emphasised that, in this work, (7) is purely a modelling assumption. Rewriting (7) gives the following condition on γ ∂x λρ c (T − T ) n ˆ · = κ n ˆ ·∇T − κ n ˆ ·∇T + C(t). (8) f ∞   R R ∂t The additional term C (t) in (8) is needed to enforce conservation of penguins, or, equiva- lently, area conservation [6], namely that ∂x n ˆ · ds = v ds = 0. (9) ∂t γ γ 2.1 Non-dimensionalisation and Parameter Values In order to identify the parameters of the problem, the model (1)-(5), (8) is non- dimensionalised, with non-dimensional quantities temporarily labelled as starred variables. Let L characterise the length scale of region R i.e. x = Lx , and use the far-field wind magnitude U as the characteristic velocity scale, u = U u . The following scalings are also used 1 T − T 1 U ∗ ∗ ∗ ∇ = ∇ ,T = ,t = t ,τ = . (10) L T − T τ L f ∞ Using (6), the non-dimensional form of the Poisson equation (5) can be written (now dropping the stars) as 4(T − T ) c f ∇ T =−α, α = . (11) T − T f ∞ The boundary condition (8) becomes ∂x λPe n ˆ · = n ˆ ·∇T − βn ˆ ·∇T + C(t), (12) ∂t where Pe = UL/D is the Péclet number, with D = κ /ρ c the thermal diffusivity of air, β = κ /κ is the ratio of thermal conductivities and C(t) is the dimensionless area conserving term, an explicit expression for which can be found by substituting (12)into(9) (see (13g) below). Note that the appearance of Pe on the left hand side of (12) affects the time scaling reported in the results. The real world dimensional time scale “t ” is therefore t = λPe (L/U )t . In what follows, the choice λ = 1 is made and justified a posteriori by comparing the time for the huddle to evolve with observations - see Sect. 5. There are three parameters: Pe, β and α. Following the justification in [6], the wind around the penguin huddle is likely turbulent and hence the Péclet number can be approxi- mated by a turbulent Reynolds number, which is O(100) for Antarctic winds at this scale. Consistent with [6], the present work considers Péclet values between 1 and 100. The values of the remaining two parameters, α and β , are estimated as follows. The air temperature far upstream is taken to be T =−40 C[2, 3] while ambient temperatures in the huddle are above freezing [11]. It seems reasonable then to let T = 0 C. The centre of a huddle can f Penguin Huddling: A Continuum Model Page 7 of 17 7 ◦ ◦ exceed 20 C, with a maximum temperature of 37.5 C having been recorded [8]. Using, therefore, T = 30 Cgives α = 3. From table C.5 of [52], the thermal conductivity of air at −40 C is approximately κ ≈ −1 −1 −2 −1 0.022Wm K . The thermal conductance of penguins is 1.3Wm K from [1], however this is not the same as the thermal conductivity - it has different units for example. To find the thermal conductivity, the conductance must be multiplied by the thickness of the −1 −1 penguin insulation layer which is d = 0.024m [53]. Therefore, κ = 0.0312Wm K and so β = 1.4182 ≈ 1.5 = O(1). This estimate for β depends on the wind flow being either laminar or turbulent, with β decreasing as the flow becomes more turbulent and the effective thermal conductivity of air increases [54]. In this work, a range of β values between β = 0 −1 and β = 1.5 are considered with β = 0.14 ≈ 0.1 = O(10 ) taken to be a representative value. To summarise, the full dimensionless system for the penguin problem is Pe v = n ˆ ·∇T − βn ˆ ·∇T + C(t), on γ, (13a) n  R ∇ φ =0in , (13b) Pe u ·∇T =∇ T in , (13c) ∇ T =−α in R, (13d) ∂φ T = T = 1, =0on γ, (13e) ∂n T → 0,φ → x as r →∞, (13f) − (n ˆ ·∇T − βn ˆ ·∇T )ds C(t) = . (13g) ds Thetaskistosolve thesystem(13a)–(13g) and, in particular, find the evolution of the γ(t) for some given initial huddle shape γ(0). This is done in Sect. 4 with the following choices: Pe = 1 → 100, β = 0 → 1.5and α = 3. 3 Numerical Method Following the work of [22, 29], the free boundary problem (13a)–(13g) can be written as a PG type equation [43] in terms of a conformal map z = f(ζ, t) from the exterior of the unit disk in the canonical ζ −plane to the exterior region  in the physical z−plane [6, 26]. However, the interior R is not mapped using f(ζ, t), thus care must be taken with the term involving the interior temperature T . First, note the following which hold on the unit disk r =|ζ|= 1[20, 22] Re f ζf t ζ 1 ∂T σ(θ) v = , n ˆ ·∇T = Re[n∇T ]= = , (14) |f | |f | ∂r |f | ζ ζ ζ iθ where ∇= ∂ − i∂ , ζ = e , f = ∂f /∂ ζ and n = n + in = ζf /|f | is the complex x y ζ x y ζ ζ representation of the normal vector in the z−plane [34]. The function σ(θ) is found by solving the following integral equation [22, 26, 29] π = K Pe(cos θ − cos θ ) σ(θ )dθ , (15) 0 7 Page 8of17 S.J.Harris, N.R. McDonald where K(x) = e K (|x|) and K is the modified Bessel function. The function σ is the heat 0 0 flux (due to the exterior effect alone) around the unit ζ −disk [22]. A numerical solution to (15) can be found at M equally spaced points on the upper half of the unit disk θ ∈ (0,π),see ∗ ∗ appendix A of [22]. To evaluate σ on the remainder of the disk, note that σ(θ ) = σ(π + θ ), where 0 ≤ θ <π,and σ(θ) → 1/π as θ → 0[22, 26]. Second, consider the contribution from the interior n ˆ ·∇T . This is found by direct com- putation in the z−plane. To begin, introduce the function u = T + αr /4; from (13d)and (13e), u satisfies ∇ u =0in R, (16) αr u = 1 + on γ. (17) This is a Laplace problem in the interior of a smooth, simply connected domain with Dirich- let boundary condition u = h(z) on γ . Efficient methods developed by Trefethen and col- leagues are here used to solve the Laplace problem directly, see [47, 55] for further details. These are based on a rational approximation of u in the complex z−plane: there exists an analytic function F(z) such that u = Re[F(z)] which can be approximated by the series solution N N 1 2 F(z) = A + A (z − c ) + , (18) 0 k 0 z − p k=1 k=1 where A , A and B are constants, N is the power series truncation, c is a point in the 0 k k 1 0 interior R (in this work, the conformal centre of the polygon is used [22]) and p are some N poles in the exterior . The first sum is known as the “Runge” (or smooth) part and the second as the “Newman” (or singular) part [49, 55]. Suitable poles p are found by the “adaptive Antoulas-Anderson” (AAA, pronounced ‘triple-A’) algorithm [49, 56] and the unknown coefficients A , A and B by a least-squares 0 k k (LS) algorithm [46, 48, 57]: combining these into a AAA-LS method [50] can thus find F(z) numerically. It follows [46]that ∇u =∇ Re[F(z)]= F (z), (19) where the primed notation indicates differentiation, in this case with respect to z. Therefore, ∇T = ∇[u − αr /4]= F (z) − z, (20) recalling that ∇= 2∂ and r = zz.Thisgives 1 ω(θ) n ˆ ·∇T = Re[n∇T ]= Re[ζf (F (z) − αz/2)]≡ , (21) R R ζ |f | |f | ζ ζ where the boundary points z in (21) are considered in the ζ −plane so that the right hand side of (21) is purely a function of θ . Drawing an analogue with σ , the function −βω is the heat flux around the unit ζ −disk due to the interior effect alone. Third, the formula (13g) for the area conserving term C(t) can be transformed to the ζ −plane, noting that ds =|f |dθ . Substituting (14)and (21)into(13g)gives − σ − βω dθ −π C(t) =  . (22) |f |dθ −π Penguin Huddling: A Continuum Model Page 9 of 17 7 Therefore, using (14), (21)and (22), the PG form of (13a)is Pe Re[f ζf ]= σ(θ) − βω(θ) + C(t)|f |, (23) t ζ ζ where σ(θ) satisfies (15)and ω(θ) = Re[ζf (F (z) − αz/2)], with F(z) given by (18). Note that σ(θ) − βω(θ) is the total heat flux across |ζ|= 1; recall that the heat flux across the huddle boundary γ is n ˆ ·∇T − βn ˆ ·∇T . It remains to find the conformal map f , the general form of which is [24] −k z = f(ζ, t) = a (t )ζ + c (t )ζ , (24) −1 k k=0 where a (t ) is a real function in time and c (t ) = a (t ) + ib (t ) are complex functions −1 k k k in time. Note that a is the conformal radius and c is the conformal centre [22]. The −1 0 same c is used as the interior point in constructing the power series solution (18)for the Poisson problem in R. Numerically, the Laurent series (24) is truncated at N terms, giving n = 2(N + 1) + 1 = 2N + 3 unknown real functions in time: a (t ), a (t ), b (t ), k = −1 k k 0, 1,...,N . Likewise, n equally spaced points around the unit ζ −disk are chosen; the first point is at θ = 0, the next N + 1 points are in the interval θ ∈ (0,π) and the final N + 1 points in θ ∈ (π, 2π). Hence the choice M = N + 1in (0,π) is made when numerically solving the integral equation (15)for σ(θ). The PG equation (23) becomes a system of n coupled ordinary differential equations (ODEs) determining the time evolution of the n Laurent coefficients (see [20, 21, 34]). In this work, the MATLAB routine ode15i was used to solve this ODE system. 4 Results In the results presented here, the evolution of the huddle boundary is plotted at equal time intervals in the range [0,t ],where t is the maximum time. The series truncation max max N = 128 is chosen, leading to a system of n = 259 ODEs to be solved. The initial start- ing shape of the penguin huddle is taken to be either a circle, a slanted ellipse, a triangle, an irregular pentagon or an hourglass. The initial Laurent coefficients a (0), c (0) for the −1 k conformal map (24) of these shapes are taken from [24] (see their Figures 2a, 2b, 2d, 5a and 7a, respectively) and scaled such that each shape encloses the same area. Area conserva- tion of the huddle was monitored during each numerical experiment and found to hold with −5 a relative error of O(10 ). As a test, experiments were run which neglected the interior (i.e. β = 0) and the area conserving (C(t)) terms; the geometry of the obtained solutions matched quantitatively with results in [22] who considered an analogous problem in this class. First, consider the “exterior only” penguin problem when β = 0. This is a continuum ana- log of the discrete model developed in [6], where the penguin huddle evolves in response to the wind but the interior reorganisation of individual penguins is not represented. Figure 2 shows results for various Péclet values and initial starting shapes of the huddle. In general, these show good agreement to the experiments of [6]: the huddle propagates in the windward direction and evolves into a more streamlined shape over time. For large time (t  t /2), max This scaling can be found analytically using the complex form of Green’s theorem [45] to find the exact area from the conformal map z = f(ζ, 0). 7 Page 10 of 17 S.J. Harris, N.R. McDonald Fig. 2 Exterior only (β = 0) penguin problem. Shapes are plotted at equal time intervals up to a maximum time t .The x and y axes (not shown) are scaled at a 1-to-1 aspect ratio. Each description refers to the max initial penguin huddle shape, the Péclet number and the value of t : a) circle with Pe = 1, t = 10, b) max max circle with Pe = 10, t = 20, c) slanted ellipse with Pe = 10, t = 20, d) irregular pentagon with Pe max max = 100, t = 50 max the huddle approaches an invariant egg-shape which continues to propagate downwind with some constant velocity. To test this observation, a root mean-squared error (RMSE) measur- ing the deviation between consecutive plotted huddle shapes is computed for the examples −4 in Fig. 2 and found to approach a constant value with a relative error of O(10 ).Thisper- mits the reasonable conclusion that huddles do evolve into a steady shape and thus in all experiments presented, the maximum time t is chosen such that the penguin huddle has max reached a steady shape prior to or at time t . max Now, consider the full penguin problem where exterior, interior (β = 0) and area con- serving terms are all included - recall from Sect. 2.1 that β = 0.1 is used. Figures 3a-e present the evolution of the same initial starting shape, a circle, with Pe = 1, 5, 10, 50 and 100. The final plot in each experiment is the steady shape and these are compared in Fig. 3f, where the shapes have been re-centered so that they share the same centre of mass. As the Péclet number is increased, the steady shape becomes less circular and more egg like. Similarly, the effect of changing the parameter β can be explored. Figures 4a-e shows the evolution of the circle with the values β = 0, 0.1, 0.5, 1 and 1.5, respectively, where Pe = 10 and t = 20 are taken for all plots. While a low value of β gives results that look max similar to the exterior only penguin problem, a high value of β causes the penguin huddle to retain a more circular shape. This is shown in Fig. 4f, where the steady shapes of Figs. 4a-e are compared. The dependence of the final steady shape of the huddle on the values of the Péclet num- ber and β is explored further. A given steady shape can be quantified by its aspect ratio (AR): the ratio between its major and minor axes with AR = 1 being a circle and increasing AR > 1 becoming more egg-shaped. Figures 5a,b show the dependence of AR on Pe and β , respectively. Each marker on the diagram represents an experiment run: each curve contains 21 markers, all run to a maximum time of t = 80, which is sufficiently long to reach the max steady state huddle shape. The results from Fig. 5 agree with the conclusions drawn from Figs. 3 and 4.First,asPe increases, AR increases and thus a higher Péclet number results in a more egg-like steady shape, with a lower value of β giving higher values of AR for all time. This is also consistent with the parameter study in [6] - see their Fig. 3 - which shows that the huddle thickness decreases with an increasing Péclet number. As Pe → 0, AR→ 1, consistent with the limit of a circular steady shape when there is no wind, as expected. Second, as β increases, AR decreases, where lower values of Pe result in lower values of the AR for all time. Penguin Huddling: A Continuum Model Page 11 of 17 7 Fig. 3 The effect of different Péclet numbers on the evolution of the huddle is shown by plotting its shape at equal time intervals. The initial starting shape of the penguin huddle is a circle, with β = 0.1and a) Pe = 1, t = 10, b) Pe = 5, t = 15, c) Pe = 10, t = 20, d) Pe = 50, t = 40, e) Pe = 100, t = 50. f) max max max max max Comparison of the steady shapes of (a)-(e), the initial (circular) huddle shape is also plotted for comparison Fig. 4 Interior β effect: initial starting shape of a circle with Pe = 10, t = 20 and a) β = 0, b) β = 0.1, max c) β = 0.5, d) β = 1, e) β = 1.5. f) Comparison of the steady shapes of (a)-(e), the initial (circular) huddle shape is also plotted for comparison 7 Page 12 of 17 S.J. Harris, N.R. McDonald Fig. 5 The aspect ratio of the final steady shape of the penguin huddle as a function of a) Péclet number, for the three given values β = 0, 0.1and 1, and b) β , for the three given values Pe = 10, 50 and 100. Each marker represents a numerical experiment Fig. 6 Initial shape effect. Shapes are plotted at equal time intervals with Pe = 10, β = 0.1, t = 35 and the max initial starting shape of a) a circle, b) a slanted ellipse, c) a triangle, d) an irregular pentagon, e) an hourglass. f) Comparison of the steady shapes of (a)-(e): all of the steady shapes are overlapped Figure 6 shows how huddles of different starting shapes evolve under the same Péclet number Pe = 10, the same interior effect β = 0.1 and up to the same maximum time t = max 35. Figures 6a-e show the evolution of a circle, a slanted ellipse, a triangle, an irregular pentagon and an hourglass, respectively, and Fig. 6f compares the steady shapes of these five experiments. The results show that, irrespective of the starting shape of the huddle and even those starting with significant concave boundaries, the same steady shape is reached. The along-boundary variation of the total heat flux across the huddle boundary of the steady shape of Fig. 6f is considered. By symmetry, only the upper half of the boundary need be considered; this is plotted in Fig. 7a. The arrows point in the normal direction of the Penguin Huddling: A Continuum Model Page 13 of 17 7 Fig. 7 a) Total heat flux across the boundary of the steady shape of Fig. 6f. The arrows point in the normal direction of the given boundary point and their magnitude is proportional to the local heat flux across the boundary. b) Canonical heat flux in the upper half ζ −disk vs θ . c) Physical heat flux in the upper half of the huddle boundary vs x given point on the boundary and their lengths represent the total physical heat flux at that point. Recall from Sect. 3 that the total heat flux in the physical domain across the huddle boundary γ is given by n ˆ ·∇T − βn ˆ ·∇T whereas the total heat flux in the canonical domain across the unit ζ −disk is given by σ(θ) − βω(θ). The canonical heat flux in the upper half ζ −disk vs θ is plotted in Fig. 7b and the physical heat flux vs x is plotted in Fig. 7c. There is a much larger heat flux on the windward side of the penguin huddle (the side of the penguin huddle directly facing the oncoming wind) compared to the leeward side. This is expected and is in agreement with the exterior temperature profile presented in [6] - see their Fig. 1 - which shows that there is a steeper temperature gradient, and hence a higher total heat flux, on the windward side. 5 Discussion In this continuum model, the evolution of the penguin huddle boundary γ(t), is affected by three factors: exterior cooling by the ambient wind, interior reorganisation of the penguins driven by their self-generation of heat, and area conservation of the huddle. The resulting two-dimensional model (13a)-(13g) can be rewritten as a PG-type equation (23)interms of the conformal map z = f(ζ, t) between the canonical (ζ ) and physical (z) planes. The map is approximated as a truncated Laurent series (24), whose unknown time-varying coefficients are found by numerically solving a system of coupled ODEs. Results from Sect. 4 show that the huddle shape eventually reaches a steady shape. As the wind strength increases, the huddle becomes increasingly egg-shaped, minimising the portion of the huddle bound- ary directly facing the oncoming wind. For a stronger interior heating effect, the huddle 7 Page 14 of 17 S.J. Harris, N.R. McDonald approaches a more circle-like steady shape, to optimise interior heat regulation. For given wind and interior effects, any initial huddle shape eventually approaches the same steady shape. There is also a greater total heat flux on the windward side of the huddle than on the leeward side. While dimensionless quantities have been used throughout the paper, it is useful to esti- mate the dimensional time scale associated with the huddle evolution. As noted in Sect. 2.1, ∗ ∗ the appropriate dimensional time scale is t = λPe(L/U )t where t is dimensionless time, and the choice λ = 1 has been made. A large huddle of, for example, 6000 penguins [1], is assumed and taking the radius of each penguin to be r ≈ 0.1m [3]gives L ≈ 10 m. Now −1 assume a sufficiently strong wind of U ≈ 20 ms for Pe = 100. Considering Fig. 3cwhere t = t = 50, the total real time for this experiment is t = 2500 seconds or ≈ 42 min- max r utes. Thus it is concluded that the penguin huddle motion described in this work occurs over roughly a one hour period, in agreement with previous modelling and field observations [6, 8, 9], and justifying the choice λ = 1. In this paper, a continuum model has been used to find the shape and propagation of the huddle boundary at the expense of tracking individual penguin movements. While there are pre-existing discrete models for penguin huddling (and other animal grouping dynam- ics) [4, 6, 7, 36], treating a penguin huddle as a continuum appears to be new. One of the advantages of the continuum model developed here is its computational efficiency, with a typical runtime of around 5-10 minutes on a standard laptop for a huddle to reach a steady shape. This can be sped up further, to a matter of seconds, by reducing the Laurent series truncation N , though care should be taken to ensure that the accuracy of the results is not compromised. When considering the “exterior only” (β = 0) penguin problem, our results matched well with those from [6], who did not include the interior effect of penguin reorganisation in their model. As discussed in Sect. 2.1, fully turbulent wind flows have β ≈ 0 and hence are well modelled by considering exterior effects alone. The representative value β = 0.1, which is likely to be an overestimate, chosen in this work demonstrates that the interior heating of the colony only provides a weak effect in determining the steady shape of the huddle, see Fig. 4f. Larger values of β resulting in more circular huddle shapes have been included here for completeness. Further, the results in this work match well with real world observations: penguin huddles are seen to march downwind over time and observed huddle shapes resemble the egg-like steady shapes found here - see [1, 4, 10, 12, 13]. Extensions to the penguin problem are of interest. A time-varying wind could be consid- ered, where the variations are on the order of an hour; when the wind direction changes, the huddle will seek to propagate back to the optimal (egg) steady shape in the new windward direction. Stochastic terms could also be incorporated into the model, as in [6]. This includes random variations in the direction and magnitude of the wind and stochasticity in the interior and exterior heat fluxes. The numerical method would need to be adapted to handle these stochastic terms. The evolution of any arbitrary starting huddle shape could also be simulated. In this work, the initial huddle boundary was taken such that its conformal map was readily represented in the form of a truncated Laurent series (24). For more arbitrary starting shapes, the hud- dle boundary can be discretised as a polygon of m vertices, then its conformal map to the ζ −plane found numerically at each time step as the penguin huddle evolves. Pre-exisiting numerical schemes for finding the conformal map between a polygon and the unit ζ −disk include the Schwarz-Christoffel toolbox [58] and a variation of the AAA-LS algorithm [50] from Sect. 3. It would also be possible to consider multi-connected huddles such as those with penguin-free holes or the interaction, and possible merger, of initially disjoint huddles. Penguin Huddling: A Continuum Model Page 15 of 17 7 While this would require a significant modification of the numerical method used here, it is worth remarking that AAA-LS type methods are able to handle such geometries [47, 49, 50]. Finally, consider the boundary condition (13a), and its associated PG equation (23). The right-hand side is constructed by exterior (n ˆ ·∇T ), interior (βn ˆ ·∇T ) and area conserva- tion (C(t)) effects. By inclusion and exclusion of these three terms, (13a) encapsulates six possible free boundary problems: 1. The dissolution problem: exterior effect only. A permeable object is placed in a uniform stream of some dissolving agent - see [22, 59]. 2. The Poisson growth problem: interior effect only. Applicable to squeeze flow and the evaporation of thin liquid films - see [45, 60, 61]. 3. The two-phase melting/freezing problem: exterior and interior effects. Porous media flow about freeze pipes where interior and exterior temperature gradients are in effect - see [24, 29]. 4. The “exterior only” penguin problem: exterior and area conservation effects. Huddle problems where interior reorganisation is not represented, such as if individuals in the centre are unable to move or there exists a “pecking order” in the population - see [6]. 5. The “bat huddle” problem: interior and area conservation effects. Huddle problems with no wind effect, for example bats huddling in caves sheltered from exterior winds. To the authors’ knowledge, there is no mathematical work on this problem, but for biological context see [39, 40]. 6. The “full” penguin problem: exterior, interior and area conservation effects. This is the full huddle continuum problem covered in this work. Therefore, the numerical method developed in this work can be used to model and simulate additional free boundary problems, such as those from the fields of fluid mechanics and mathematical biology. Acknowledgements The authors thank Carole Hall, Stony Brook University, for their inspiring talk on Adélie penguins and the useful discussion thereafter. The authors also thank the anonymous reviewers for their useful comments and suggestions. Author Contribution Both authors contributed equally to this work. Funding Samuel J. Harris was supported by a UK Engineering and Physical Sciences Research Council PhD studentship, grant numbers EP/N509577/1 and EP/T517793/1. Data Availability All data generated or analysed during this study are included in this published article and its supplementary information files. Code Availability All relevant code files are publicly available at the following GitHub link: https://github. com/Sam-J-Harris/Penguin-huddling-a-continuum-model.git. Declarations Competing Interests Not applicable. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. 7 Page 16 of 17 S.J. Harris, N.R. McDonald References 1. Le Maho, Y.: The emperor penguin: a strategy to live and breed in the cold: morphology, physiology, ecology, and behavior distinguish the polar emperor penguin from other penguin species, particularly from its close relative, the king penguin. Am. Sci. 65(6), 680–693 (1977) 2. McCafferty, D.J., Gilbert, C., Thierry, A.-M., Currie, J., Le Maho, Y., Ancel, A.: Emperor penguin body surfaces cool below air temperature. Biol. Lett. 9(3), 20121192 (2013) 3. Williams, C.L., Hagelin, J.C., Kooyman, G.L.: Hidden keys to survival: the type, density, pattern and functional role of emperor penguin body feathers. Proc. R. Soc. Lond. B 282(1817), 20152033 (2015) 4. Gerum, R.C., Fabry, B., Metzner, C., Beaulieu, M., Ancel, A., Zitterbart, D.P.: The origin of traveling waves in an emperor penguin huddle. New J. Phys. 15(12), 125022 (2013) 5. Kooyman, G.L., Gentry, R.L., Bergman, W.P., Hammel, H.T.: Heat loss in penguins during immersion and compression. Comp. Biochem. Physiol. 54(1), 75–80 (1976) 6. Waters, A., Blanchette, F., Kim, A.D.: Modeling huddling penguins. PLoS ONE 7(11), 50277 (2012) 7. Gerum, R., Richter, S., Fabry, B., Le Bohec, C., Bonadonna, F., Nesterova, A., Zitterbart, D.P.: Structural organisation and dynamics in king penguin colonies. J. Phys. D, Appl. Phys. 51(16), 164004 (2018) 8. Gilbert, C., Robertson, G., Le Maho, Y., Naito, Y., Ancel, A.: Huddling behavior in emperor penguins: dynamics of huddling. Physiol. Behav. 88(4–5), 479–488 (2006) 9. Kirkwood, R., Robertson, G.: The occurrence and purpose of huddling by emperor penguins during foraging trips. Emu 99(1), 40–45 (1999) 10. Ancel, A., Gilbert, C., Poulin, N., Beaulieu, M., Thierry, B.: New insights into the huddling dynamics of emperor penguins. Anim. Behav. 110, 91–98 (2015) 11. Zitterbart, D.P., Wienecke, B., Butler, J.P., Fabry, B.: Coordinated movements prevent jamming in an emperor penguin huddle. PLoS ONE 6(6), 20260 (2011) 12. Mina, T., Min, B.-C.: Penguin huddling inspired distributed boundary movement for group survival in multi-robot systems using Gaussian processes. In: 2018 IEEE International Conference on Robotics and Biomimetics, pp. 2177–2183. IEEE (2018) 13. Richter, S., Gerum, R., Winterl, A., Houstin, A., Seifert, M., Peschel, J., Fabry, B., Le Bohec, C., Zit- terbart, D.P.: Phase transitions in huddling emperor penguins. J. Phys. D, Appl. Phys. 51(21), 214002 (2018) 14. Gu, W., Christian, J.K., Woodson, C.B.: A novel coupled fluid-behavior model for simulating dynamic huddle formation. PLoS ONE 13(8), 0203231 (2018) 15. Saffman, P.G., Taylor, G.I.: The penetration of a fluid into a porous medium or Hele-Shaw cell containing a more viscous liquid. Proc. R. Soc. A 245(1242), 312–329 (1958) 16. Howison, S.: Fingering in Hele-Shaw cells. J. Fluid Mech. 167, 439–453 (1986) 17. Paterson, L.: Radial fingering in a Hele Shaw cell. J. Fluid Mech. 113, 513–529 (1981) 18. Taylor, G., Saffman, P.G.: A note on the motion of bubbles in a Hele-Shaw cell and porous medium. Q. J. Mech. Appl. Math. 12(3), 265–279 (1959) 19. Entov, V.M., Etingof, P.I.: Bubble contraction in Hele-Shaw cells. Q. J. Mech. Appl. Math. 44(4), 507–535 (1991) 20. Dallaston, M.C., McCue, S.W.: An accurate numerical scheme for the contraction of a bubble in a Hele– Shaw cell. ANZIAM J. 54, 309–326 (2013) 21. Dallaston, M.C., McCue, S.W.: A curve shortening flow rule for closed embedded plane curves with a prescribed rate of change in enclosed area. Proc. R. Soc. A 472(2185), 20150629 (2016) 22. Ladd, A.J.C., Yu, L., Szymczak, P.: Dissolution of a cylindrical disk in Hele-Shaw flow: a conformal- mapping approach. J. Fluid Mech. 903, 46 (2020) 23. Cummings, L.M., Hohlov, Y.E., Howison, S.D., Kornev, K.: Two-dimensional solidification and melting in potential flows. J. Fluid Mech. 378, 1–18 (1999) 24. Rycroft, C.H., Bazant, M.Z.: Asymmetric collapse by dissolution or melting in a uniform flow. Proc. R. Soc. A 472(2185), 20150531 (2016) 25. Mullins, W.W., Sekerka, R.F.: Stability of a planar interface during solidification of a dilute binary alloy. J. Appl. Phys. 35(2), 444–451 (1964) 26. Choi, J., Margetis, D., Squires, T.M., Bazant, M.Z.: Steady advection–diffusion around finite absorbers in two-dimensional potential flows. J. Fluid Mech. 536, 155–184 (2005) 27. Tsai, V.C., Wettlaufer, J.S.: Star patterns on lake ice. Phys. Rev. E 75(6), 066105 (2007) 28. Grodzki, P., Szymczak, P.: Reactive-infiltration instability in radial geometry: from dissolution fingers to star patterns. Phys. Rev. E 100(3), 033108 (2019) 29. Goldstein, M.E., Reid, R.L.: Effect of fluid flow on freezing and thawing of saturated porous media. Proc. R. Soc. A 364(1716), 45–73 (1978) 30. Langer, J.S.: Instabilities and pattern formation in crystal growth. Rev. Mod. Phys. 52(1), 1 (1980) Penguin Huddling: A Continuum Model Page 17 of 17 7 31. Brower, R.C., Kessler, D.A., Koplik, J., Levine, H.: Geometrical models of interface evolution. Phys. Rev. A 29(3), 1335 (1984) 32. Sethian, J.A.: Curvature and the evolution of fronts. Commun. Math. Phys. 101(4), 487–499 (1985) 33. Hilton, J.E., Miller, C., Sharples, J.J., Sullivan, A.L.: Curvature effects in the dynamic propagation of wildfires. Int. J. Wildland Fire 25(12), 1238–1251 (2016) 34. Harris, S.J., McDonald, N.R.: Fingering instability in wildfire fronts. J. Fluid Mech. 943, 34 (2022) 35. Burger, M., Haškovec, J., Wolfram, M.-T.: Individual based and mean-field modeling of direct aggrega- tion. Physica D 260, 145–158 (2013) 36. Bernardi, S., Scianna, M.: An agent-based approach for modelling collective dynamics in animal groups distinguishing individual speed and orientation. Philos. Trans. R. Soc. Lond. B 375(1807), 20190383 (2020) 37. Katz, Y., Tunstrøm, K., Ioannou, C.C., Huepe, C., Couzin, I.D.: Inferring the structure and dynamics of interactions in schooling fish. Proc. Natl. Acad. Sci. 108(46), 18720–18725 (2011) 38. Bhattacharya, K., Vicsek, T.: Collective decision making in cohesive flocks. New J. Phys. 12(9), 093019 (2010) 39. Herreid, C.F.: Temperature regulation of Mexican free-tailed bats in cave habitats. J. Mammal. 44(4), 560–573 (1963) 40. Ryan, C.C., Burns, L.E., Broders, H.G.: Changes in underground roosting patterns to optimize energy conservation in hibernating bats. Can. J. Zool. 97(11), 1064–1070 (2019) 41. Nave, G.K. Jr, Mitchell, N.T., Chan Dick, J.A., Schuessler, T., Lagarrigue, J.A., Peleg, O.: Attraction, dynamics, and phase transitions in fire ant tower-building. Front. Robot. AI 7, 25 (2020) 42. Ko, T.-Y.H., Yu, H.D.L.: Fire ant rafts elongate under fluid flows. Bioinspir. Biomim. 17(4), 045007 (2022) 43. Bazant, M.Z., Crowdy, D.: Conformal mapping methods for interfacial dynamics. In: Handbook of Ma- terials Modeling, pp. 1417–1451. Springer, Dordrecht (2005) 44. Gustafsson, B., Vasil’ev, A.: Conformal and Potential Analysis in Hele-Shaw Cells. Springer, Berlin (2006) 45. McDonald, R., Mineev-Weinstein, M.: Poisson growth. Anal. Math. Phys. 5(2), 193–205 (2015) 46. Trefethen, L.N.: Series solution of Laplace problems. ANZIAM J. 60(1), 1–26 (2018) 47. Trefethen, L.N.: Numerical conformal mapping with rational functions. Comput. Methods Funct. Theory 20(3), 369–387 (2020) 48. Baddoo, P.J.: Lightning solvers for potential flows. Fluids 5(4), 227 (2020) 49. Costa, S.: Solving Laplace problems with the AAA algorithm (2020). arXiv:2001.09439. ArXiv preprint 50. Costa, S., Trefethen, L.N.: AAA-least squares rational approximation and solution of Laplace problems (2021). arXiv:2107.01574. ArXiv preprint 51. Gupta, S.C.: The Classical Stefan Problem: Basic Concepts, Modelling and Analysis with Quasi- Analytical Solutions and Methods, vol. 45. Elsevier, Amsterdam (2017) 52. Basu, P.: Biomass Gasification, Pyrolysis and Torrefaction: Practical Design and Theory. Academic Press, Boston (2018) 53. Dawson, C., Vincent, J.F.V., Jeronimidis, G., Rice, G., Forshaw, P.: Heat transfer through penguin feath- ers. J. Theor. Biol. 199(3), 291–295 (1999) 54. Kittel, C., McEuen, P., McEuen, P.: Introduction to Solid State Physics, vol. 8. Wiley, New York (1996) 55. Gopal, A., Trefethen, L.N.: Solving Laplace problems with corner singularities via rational functions. SIAM J. Numer. Anal. 57(5), 2074–2094 (2019) 56. Nakatsukasa, Y., Sète, O., Trefethen, L.N.: The AAA algorithm for rational approximation. SIAM J. Sci. Comput. 40(3), 1494–1522 (2018) 57. Brubeck, P.D., Nakatsukasa, Y., Trefethen, L.N.: Vandermonde with Arnoldi. SIAM Rev. 63(2), 405–415 (2021) 58. Driscoll, T.A.: Algorithm 756: a MATLAB toolbox for Schwarz-Christoffel mapping. ACM Trans. Math. Softw. 22(2), 168–186 (1996) 59. Dutka, F., Starchenko, V., Osselin, F., Magni, S., Szymczak, P., Ladd, A.J.: Time-dependent shapes of a dissolving mineral grain: comparisons of simulations with microfluidic experiments. Chem. Geol. 540, 119459 (2020) 60. Agam, O.: Viscous fingering in volatile thin films. Phys. Rev. E 79(2), 021603 (2009) 61. Crowdy, D., Kang, H.: Squeeze flow of multiply-connected fluid domains in a Hele-Shaw cell. J. Non- linear Sci. 11(4), 279–304 (2001) Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Journal

Acta Applicandae Mathematicae – Springer Journals

Published: Jun 1, 2023

Keywords: Free boundary problems; Conformal mapping; AAA-least squares algorithm; Fluid mechanics; Collective behaviour; Huddling

There are no references for this article.