Non-Local Formulation of Heat Transfer with Phase Change in Domains with Spherical and Axial Symmetries

Abstract

Describing heat transfer in domains with strong non-linearities and discontinuities, e.g. propagating fronts between different phases, or growing cracks, is a challenge for classical approaches, where conservation laws are formulated as partial differential equations subsequently solved by discretisation methods such as the finite element method (FEM). An alternative approach for such problems is based on the non-local formulation; a prominent example is peridynamics (PD). Its numerical implementation however demands substantial computational resources for problems of practical interest. In many engineering situations, the problems of interest may be considered with either axial or spherical symmetry. Specialising the non-local description to such situations would decrease the number of PD particles by several orders of magnitude with proportional decrease of the computational time, allowing for analyses of larger domains or with higher resolution as required. This work addresses the need for specialisation by developing bond-based peridynamic formulations for physical problems with axial and spherical symmetries. The development is focused on the problem of heat transfer with phase change. The accuracy of the new non-local description is verified by comparing the computational results for several test problems with analytical solutions where available, or with numerical solutions by the finite element method.

1 Introduction

The solutions of many engineering problems require an accurate description of heat transfer with phase change, e.g. manufacturing of structures and machine parts [1, 2], construction and geotechnics projects [3, 4], and latent heat thermal energy storage [5, 6].

The mathematical description of this problem needs to capture strong non-linearities and discontinuities, which are either challenging for, or outside the remit of, the classical local formulations, i.e. those using partial differential equations to describe the conservation laws. Peridynamics (PD) provides a description that can deal with strong non-linearities and capture the evolution of discontinuities [7]. In PD, the conservation laws are represented by a set of spatially integral and temporally differential equations [8].

The PD formulation of diffusion was first developed in [9], which led to the classical bond-based formulations of transient heat transfer in [10, 11]; see also [12, 13]. These have been used to develop coupled models that take into account various physical factors and that cannot be solved by other numerical methods, e.g. [14,15,16,17]. Specifically, bond-based PD formulations for heat transfer with phase change have been proposed recently: for modelling of soils [7] and for welding simulations [18].

To date, PD simulations have been predominantly performed with 1D and 2D models using rectangular grids of peridynamic particles because 3D models of engineering-scale problems require significant computational power. Even when modern numerical algorithms for solution of PD equations on GPUs realise almost linear dependence between the number of PD particles and the computational time, as was shown in [19], the computational time for 3D models that include millions of particles can be extremely large. However, many engineering and physical problems involve symmetries in geometry and boundary conditions, which allow for simplifying real 3D problems to corresponding 2D problems in the case of axial symmetry, and 1D problems in the case of spherical symmetry. Accounting for such symmetries can decrease the number of PD particles by several orders of magnitude with a proportional decrease in the required computational time. At present, such specialisations have been developed only for state-based peridynamics for mechanical problems. For example, in [20], an axisymmetric ordinary state-based peridynamic model for linear elastic solids was originally formulated. Later works considered shear deformations in brittle solids [21], crack deflection in ceramic composites [22], and thermal cracking [23]. Bond-based PD formulations for diffusion problems with symmetries have not been developed.

The contribution of this work is a bond-based PD description of transient heat transfer with phase change in domains with axial and spherical symmetries. The formulation can be used for analysis of heat transfer with and without the phase change. The specialisation to problems with axial symmetry can be used to study: the process of energy extraction from soils and rocks by borehole heat exchangers or energy piles [24]; artificial ground freezing by a single freeze pipe [3, 25, 26]; and heat transfer in insulated wires and buried tunnels and pipes [27, 28]. Since heat transfer in a substance and mass transport in porous media are described by similar mathematical expressions, the proposed formulation can be applied to study: the process of water flow in soils around a single water well or drainage in open pits and tunnels [29]; and the process of wood drying [30]. A prospective application of the model can be found in the description of oil or natural gas wells with and without hydraulic fracturing of rock formation [31]. A similar form of the peridynamic equation can be used to describe chemical diffusion. It should be noted that, for porous media such as soils, the present model can be extended to describe coupled physical problems that include the process of heat conduction or chemical diffusion associated with water flow, e.g. [7, 32, 33], i.e. problems involving heat convection and/or chemical advection. Specialisation to problems with spherical symmetry can be used to describe the heat (or with modification, moisture or solute) transfer in spherical bodies, such as bubbles and droplets, which can be important in problems of hydrodynamics and chemical engineering.

The development is based on the classical local formulation of heat transfer in the enthalpy form which considers thermo-physical properties of materials as discontinuous functions [34]. This allows for a bond-based PD formulation that can be used to describe materials with sharp interfaces between different phases, i.e. without a transition zone. The use of a transition zone can be numerically challenging, as discussed in a previous paper [7], where the model was based on the apparent heat capacity approach. For models with transition zones, the solution [7] can in principle be adapted to describe problems with axial and spherical symmetries.

The paper is structured as follows. Section 2 presents the bond-based PD formulation for heat conduction with phase change in domains with axial and spherical symmetries, and its numerical implementation. Section 3 presents the analysis results for several test problems and comparisons with analytical and numerical solutions to verify the developed model. Conclusions are drawn in Sect. 4.

2 Mathematical Formulations and Numerical Implementation

2.1 General Bond-Based Peridynamics for Heat Conduction with Phase Change

A body occupying a region $\mathbf {R}$ is considered to be a collection of distinct particles labelled by position vectors $\mathbf {r}$ with respect to a background coordinate system. Each particle has associated volume, mass, and enthalpy H, corresponding to temperature T. A particle at position r interacts with (and is connected to) all particles at positions $\mathbf {r}'$ within a horizon $\mathbf {\mathcal {H}_r}$. The term ’bond’ refers to the interactions between two particles at positions $\mathbf {r}$ and $\mathbf {r}' \in \mathbf {\mathcal {H}_r}$. Considering the processes to be modelled, the bonds will be called thermal or ’t-bonds’. The peridynamic heat flux density (flux per unit volume) along a t-bond depends on the distance between points ${\textbf {r}}$ and ${\textbf {r'}}$, represented by a distance vector $\varvec{\xi } =\left( \mathbf {r'} - \mathbf {r}\right)$.

The presence of water in two phases, liquid and solid/ice, separates $\mathbf {R}$ into two regions—solid and liquid. The boundary between the two regions is defined by the position of the isotherm with temperature $T_f$, which is the phase change temperature. Depending on the region, where connected particles are located, three different t-bonds are considered: solid t-bonds, interfacial t-bonds connecting particles from different regions, and liquid t-bonds.

The peridynamic formulation allows for using discontinuous material properties, specifically here to describe the jump in thermo-physical properties across the phase transition boundary. The thermal conductivity of a particle is defined by:

$$\begin{aligned} \lambda (\mathbf {r},t) = {\left\{ \begin{array}{ll} \lambda _f, &{} T<T_f \\ \frac{1}{2} (\lambda _f + \lambda _m), &{} T=T_f \\ \lambda _m, &{} T_f<T \end{array}\right. }\,. \end{aligned}$$

(1)

where $\lambda _f$ and $\lambda _m$ are the thermal conductivity values of the frozen and unfrozen (melted) materials, respectively.

The volumetric heat capacity of a particle is defined by:

$$\begin{aligned} \rho C (\mathbf {r},t) = {\left\{ \begin{array}{ll} \rho _f C_f, &{} T<T_f \\ \rho _m C_m, &{} T_f<T \end{array}\right. }\,. \end{aligned}$$

(2)

where $\rho _f$ and $\rho _m$ are the densities, and $C_f$ and $C_m$ are the specific heat capacities of the frozen and unfrozen (melted) materials, respectively.

If the enthalpy of a particle is known, Eq. (2) can be used to define the particle temperature by:

$$\begin{aligned} T (\mathbf {r},t) = {\left\{ \begin{array}{ll} \frac{H (\mathbf {r} ,t)}{\rho _f C_f} , &{} H (\mathbf {r} ,t)< \rho _f C_f T_f \\ T_f , &{} \rho _f C_f T_f \le H (\mathbf {r} ,t) \le \rho _f C_f T_f + \rho _f \Pi L \\ \frac{H (\mathbf {r} ,t) - (\rho _f C_f T_f + \rho _f \Pi L)}{\rho _m C_m}, &{} \rho _f C_f T_f + \rho _f \Pi L < H (\mathbf {r} ,t) \end{array}\right. }\,. \end{aligned}$$

(3)

where L is the specific latent heat of solidification of the material and $\Pi$ is the fraction of the material volume affected by the phase change defined for multi-component materials.

The conductive heat flux along a t-bond, $\mathbf {J}\left( T,\mathbf {r}',\mathbf {r},t\right)$, is defined by [32, 33]:

$$\begin{aligned} \mathbf {J} = \overline{\lambda }\left( \mathbf {r}', \mathbf {r}, t\right) \frac{T\left( \mathbf {r}',t\right) - T\left( \mathbf {r}, t\right) }{\left\| \varvec{\xi } \right\| } \cdot \frac{\varvec{\xi }}{\left\| \varvec{\xi } \right\| } \end{aligned}$$

(4)

where $T\left( \mathbf {r}', t\right)$ and $T\left( \mathbf {r}, t\right)$ are the temperatures at points $\mathbf {r}'$ and $\mathbf {r}$ at time t, respectively; $\overline{\lambda }\left( \mathbf {r}', \mathbf {r},t\right)$ is the thermal conductivity of the t-bond; $\varvec{\xi } / \left\| \varvec{\xi } \right\|$ is the unit vector along the t-bond ($\left\| \varvec{\xi } \right\|$ is the length of the distance vector).

The thermal conductivity of a t-bond depends on the locations of the two connected particles and is defined by:

$$\begin{aligned} \overline{\lambda }\left( \mathbf {r}', \mathbf {r}, t\right) =\frac{2 \lambda (\mathbf {r},t) \lambda (\mathbf {r}',t)}{\lambda (\mathbf {r},t) + \lambda (\mathbf {r}',t)}, \end{aligned}$$

(5)

where $\lambda (\mathbf {r},t)$ and $\lambda (\mathbf {r}',t)$ are the thermal conductivities at the peridynamic particles located at $\mathbf {r}$ and $\mathbf {r}'$, respectively.

2.2 Specialisation to Problems with Axial Symmetry

In the case of axisymmetric heat transfer, $\mathbf {R}$ is a cylindrical region described by coordinates $(r, \theta , z)$, where z is the axis of symmetry; see Fig. 1. The axial symmetry of the process allows for defining an axisymmetric peridynamic particle at position $\mathbf {r}$ as a circle with centre (0, 0, z) and radius r, i.e. such particle is labelled by two coordinates (r, z). The distance vector between two axisymmetric peridynamic particles is independent of $\theta$ with the length equal to the distance between their traces on a given $r-z$ plane. Following the tradition for selecting horizon shapes, the horizon of an axisymmetric particle at position $\mathbf {r}$ is a torus with a circular cross-section of radius $\delta$ (in $r-z$ planes) and a radius or revolution r (in $r-\theta$ planes). An axisymmetric t-bond between two axisymmetric peridynamic particles is a surface of revolution formed by rotating the distance vector between the particles around the z-axis, i.e. a frustum of a cone. Fixing an $r-z$ plane, the axisymmetric peridynamic particles (circles) are put in correspondence with their traces (points), the axisymmetric t-bonds (frusta of cones) correspond to their traces (line segments), and the horizons of axisymmetric particles (tori) are put in correspondence with their traces (circles) on that plane.

Fig. 1

Illustration of phase regions, particle horizons, and t-bond types for axisymmetric heat transfer with phase change

The energy conservation at t-bonds for the axisymmetric case is formulated following the classical approach [10]:

$$\begin{aligned} V(\mathbf {r}', \mathbf {r}) \frac{\partial \overline{H}(\mathbf {r}', \mathbf {r}, t)}{\partial t}= S(\mathbf {r}', \mathbf {r}) J\left( \mathbf {r}', \mathbf {r}, t\right) \end{aligned}$$

(6)

where $\overline{H}(\mathbf {r}', \mathbf {r}, t)$ is the average enthalpy of the t-bond, associated with the average temperature, heat capacity and density of the t-bond, but the volume $V(\mathbf {r}', \mathbf {r})$, and the cross-section area $S(\mathbf {r}', \mathbf {r})$ of the t-bond are specific to the axisymmetric case. Their calculation proceeds as follows.

The t-bond, being the frustum of a cone, is considered with thickness h. Its volume is then given by:

$$\begin{aligned} {\begin{matrix} V(\mathbf {r}', \mathbf {r}) = \pi (r' + r) {\left\| \varvec{\xi } \right\| } h \end{matrix}} \end{aligned}$$

(7)

where r and $r'$ are the radial coordinates of the particles $\mathbf {r}$ and $\mathbf {r}'$.

The amount of heat transferred to and accumulated in a t-bond can be assumed proportional to the area of the outer surface, so that the cross-section area of a t-bond $S(\mathbf {r}', \mathbf {r})$ is given by:

$$\begin{aligned} {\begin{matrix} S(\mathbf {r}', \mathbf {r}) = 2 \pi r' h \end{matrix}} \end{aligned}$$

(8)

Substitution of Eqs. (7) and (8) into (6) leads to:

$$\begin{aligned} \frac{\partial \overline{H}(\mathbf {r}', \mathbf {r}, t)}{\partial t}= \Omega (\mathbf {r}', \mathbf {r}) \frac{J\left( \mathbf {r}', \mathbf {r}, t\right) }{{\left\| \varvec{\xi } \right\| }} \end{aligned}$$

(9)

where $\Omega (\mathbf {r}', \mathbf {r})$ is a shape factor given by:

$$\begin{aligned} \Omega (\mathbf {r}', \mathbf {r}) = \frac{2 r'}{r' + r} \end{aligned}$$

(10)

Inserting Eq. (4) into Eq. (6) provides the energy conservation in the t-bond in terms of enthalpy in axisymmetric form. The energy conservation for a particle at position $\mathbf {r}$ involves the fluxes in all adjacent t-bonds (bonds to particles within the horizon $\mathbf {\mathcal {H}_r}$), and the energy conservation equation for a single bond is found by integrating over the horizon and given by:

$$\begin{aligned} {\begin{matrix} {\int _{\mathbf {\mathcal {H}_r}} \frac{\partial \overline{H}(\mathbf {r}', \mathbf {r}, t)}{\partial t}} \mathrm{d}V_{\mathbf {r}'} = \int _{\mathbf {\mathcal {H}_r}} \Omega (\mathbf {r}', \mathbf {r}) \overline{\lambda }\left( \mathbf {r}',\mathbf {r},t\right) \frac{T\left( \mathbf {r}',t\right) -T\left( \mathbf {r},t\right) }{\left\| \varvec{\xi }\right\| ^{2} } \mathrm{d}V_{\mathbf {r}'} \end{matrix}} \end{aligned}$$

(11)

where $V_{\mathbf {r}'}$ is the horizon volume of the particle at position $\mathbf {r}'$.

The relationship between the enthalpy at point $\mathbf {r}$ and time t and the enthalpy in all adjacent bonds is given by [10]:

$$\begin{aligned} {\int _\mathbf {\mathcal {H}_r} \frac{\overline{H}(\mathbf {r}', \mathbf {r}, t)}{\partial t}} \mathrm{d}V_{\mathbf {r}'} = \frac{\partial H(\mathbf {r},t)}{\partial t} V_{\mathbf {r} } \end{aligned}$$

(12)

where $V_{\mathbf {r}}$ is the horizon volume of particle $\mathbf {r}$.

Combining Eqs. (4)–(12), and taking into account that an additional heat sink source $Q \left( \mathbf {r},t\right)$ may exist within the horizon, the following equation for the evolution of enthalpy at particle $\mathbf {r}$ is derived:

$$\begin{aligned} \frac{\partial H(\mathbf {r},t)}{\partial t} = \int _{\mathbf {\mathcal {H}_r}} \Omega (\mathbf {r}', \mathbf {r}) \overline{\Lambda }\left( \mathbf {r}',\mathbf {r},t\right) \frac{T\left( \mathbf {r}',t\right) -T\left( \mathbf {r},t\right) }{\left\| \varvec{\xi }\right\| ^{2} }\mathrm{d}V_{\mathbf {r}'} + Q \left( \mathbf {r},t\right) \end{aligned}$$

(13)

where $\overline{\Lambda }\left( \mathbf {r}',\mathbf {r},t\right) = {\overline{\lambda }\left( \mathbf {r}',\mathbf {r},t\right) }/{V_{\mathbf {r} } }$ is the peridynamic micro-conductivity.

The last equation is the peridynamic energy conservation equation for axisymmetric heat transfer. It can be used to solve 1D and 2D axisymmetric problems.

The relations between the peridynamic micro-conductivity parameter and the macroscopic conductivity in the case of constant influence function have been derived for a 1D problem in [10] and for the 2D problem in [11]. These are:

$$\begin{aligned} \begin{aligned} \overline{\Lambda }&= \frac{\overline{\lambda }}{\delta }, \qquad (1\text {D case}) \\ \overline{\Lambda }&= \frac{4 \overline{ \lambda }}{\pi \delta ^2}, \qquad (2\text {D case}) \end{aligned} \end{aligned}$$

(14)

The approach presented in [11] and illustrated in Fig. 2a is taken to define the conductive heat flux $J^r(\mathbf {r},t)$ in the radial direction for a particle at position $\mathbf {r} = (R_{cs}, z)$, i.e. the heat flux through some cylindrical surface with the radius $R_{cs}$. This approach takes into account only the bonds that connect the particle $\mathbf {r}$ with particles $\mathbf {r}'$ with higher temperatures. These are particles within a restricted horizon denoted by $\mathbf {\mathcal {H}_r^{+}}$, and the heat flux between them and the particle $\mathbf {r}$ is given by:

$$\begin{aligned} {\begin{matrix}{ J^r (r = R_{cs}, z, t) = \int _{\mathbf {\mathcal {H}_r^{+}}} \Omega (\mathbf {r}', \mathbf {r}) \overline{\Lambda }\left( \mathbf {r}',\mathbf {r},t\right) \frac{T\left( \mathbf {r}',t\right) -T\left( \mathbf {r},t\right) }{\left\| \varvec{\xi }\right\| }cos(\phi ) \mathrm{d}V_{\mathbf {r}'}} \end{matrix}} \end{aligned}$$

(15)

where $\phi$ is the angle between the bond direction and the r-axis and $cos(\phi ) = \frac{r - r'}{\left\| {\textbf {r}} - {\textbf {r}}'\right\| }$, $(r - r')$ is the radial distance between particles ${\textbf {r}}$ and ${\textbf {r}}'$ and $\left\| {\textbf {r}} - {\textbf {r}}'\right\|$ is the Euclidean distance. If the particle $\mathbf {r} = (r, z)$ does not belong to the cylindrical surface, but is located close to it, i.e. $R_{cs} \ne r$, then, it can be assumed that $2 \pi R_{cs} J^r ( R_{cs}, z, t) = 2 \pi r J^r (r, z, t)$. This allows for establishing $J^r ( R_{cs} \ne r, z, t) = \frac{r}{R_{cs}} J^r (r, z, t)$.

The conductive heat flux $J^z (\mathbf {r}, t)$ in the vertical direction for a particle at position $\mathbf {r} = (r, Z_{hs})$ , i.e. the heat flux through a horizontal surface at the level $Z_{hs}$, is defined similarly, Fig. 2b, by:

$$\begin{aligned} {\begin{matrix}{ J^z (r, z = Z_{hs}, t) = \int _{\mathbf {\mathcal {H}_r^{+}}} \Omega (\mathbf {r}', \mathbf {r}) \overline{\Lambda }\left( \mathbf {r}',\mathbf {r},t\right) \frac{T\left( \mathbf {r}',t\right) -T\left( \mathbf {r},t\right) }{\left\| \varvec{\xi }\right\| }cos(\psi ) \mathrm{d}V_{\mathbf {r}'}} \end{matrix}} \end{aligned}$$

(16)

where $\psi$ is the angle between the bond direction and the z-axis.

Fig. 2

Radial (a) and vertical (b) heat fluxes at the particle ${\textbf {r}}$

2.3 Specialisation to Problems with Spherical Symmetry

In the case of spherical symmetry, a peridynamic particle represents a sphere of a given radius r, centred at the origin. The horizon of a particle is the shell of a ball enclosed between two spheres of radii $r+\delta$ and $r-\delta$. Thus, a 3D spherical domain is represented by a 1D set of particles with associated volumes and masses; see Fig. 3. Every particle is characterised by its own enthalpy H and temperature T. A particle at position r interacts with all particles at positions $r'$ within its horizon $\mathbf {\mathcal {H}_r}$.

Fig. 3

Illustration of phase regions in heat transfer problem with spherical symmetry: a a spherical body; b peridynamic bonds and horizon

The energy conservation at t-bonds for the spherically symmetric case is given by:

$$\begin{aligned} V(r', r) \frac{\partial \overline{H}(r', r, t)}{\partial t}= S(r', r) \mathbf {J}\left( r', r, t\right) \end{aligned}$$

(17)

where $\mathbf {J}\left( r',r,t\right)$, is defined by Eq. (4). The volume $V(r', r)$ and the cross-section area $S(r', r)$ of the t-bond are specific to the spherically symmetric case.

The volume of the t-bond $V(r', r)$ is given by:

$$\begin{aligned} {\begin{matrix} V (r', r) = \frac{4}{3} \pi (r'^3 - r^3) = \frac{4}{3} \pi (r'^2 + r' r + r^2) {\left\| \varvec{\xi } \right\| } \end{matrix}} \end{aligned}$$

(18)

Similarly to the case of axial symmetry, the cross-section of the t-bond is identified with the outer surface of the bond, i.e., $S(r', r)$ is given by:

$$\begin{aligned} {\begin{matrix} S(r', r) = 4 \pi r'^2 \end{matrix}} \end{aligned}$$

(19)

Substitution of Eqs. (18) and (19) into (17) gives:

$$\begin{aligned} \frac{\partial \overline{H}(r', r, t)}{\partial t} = \Omega _{sp} (r', r) \frac{\mathbf {J}\left( r', r, t\right) }{{\left\| \varvec{\xi } \right\| }} \end{aligned}$$

(20)

where $\Omega _{sp} (r', r)$ is the shape factor given by:

$$\begin{aligned} \Omega _{sp} (r', r) = \frac{3 r'^2}{r^2 + r' r + r'^2} \end{aligned}$$

(21)

Following the procedure explained in Sect. 2.2 for the axisymmetric case (see Eqs. (11) and (12)), the evolution of temperature at particle r for the case of spherical symmetry is:

$$\begin{aligned} \frac{\partial H(r,t)}{\partial t} = \int _{\mathbf {\mathcal {H}_r}} \Omega _{sp} (r', r) \overline{\Lambda }\left( r',r,t\right) \frac{T\left( r',t\right) -T\left( r,t\right) }{\left\| \varvec{\xi }\right\| ^{2} }\mathrm{d}V_{r'} + Q \left( r,t\right) \end{aligned}$$

(22)

2.4 Numerical Implementation

In bond-based peridynamics, the volumes associated with (occupied by) different particles can be identical or different. This means that the particles can be located at the nodes of a regular grid in the former case, or at the nodes of an irregular grid in the latter case. In the present paper, only the case of a uniform rectangular grid is considered, such that the volume of every particle $\mathbf {r}$ is determined by the grid spacing ${\Delta _\mathbf {r}}$. The radius of the horizon $\delta$ is defined as three grid spacings, e.g. $\delta = 3 {\Delta _\mathbf {r}}$. The temporal discretization is implemented by considering equally spaced time instances, $t^n$ ($n \in \mathbb {N}$), with time interval between the instances ${\Delta _t} = t^{n+1} - t^n$.

The spatial discretization of the conservation of energy, Eq. (13), for a particle at position $\mathbf {r}_{\alpha }$ is

$$\begin{aligned} \frac{\partial H(\mathbf {r}_{\alpha } ,t)}{\partial t} = \sum _{\mathbf {r}_{\beta }} \Omega (\mathbf {r}_{\beta }, \mathbf {r}_{\alpha }) {\overline{\Lambda }\left( \mathbf {r}_{\beta },\mathbf {r}_{\alpha },t\right) } \frac{T\left( \mathbf {r}_{\beta } ,t\right) -T\left( \mathbf {r}_{\alpha } ,t\right) }{\left\| \mathbf {r}_{\beta } -\mathbf {r}_{\alpha } \right\| ^{2} } V_{\alpha \beta } \end{aligned}$$

(23)

where $V_{\alpha \beta }$ is the part of the volume of particle $\beta$ that is located within the horizon of particle $\alpha$.

The time derivative of temperature in Eq. (23) is approximated by the forward Euler method. The fully discretised conservation of energy becomes

$$\begin{aligned} \begin{aligned}&\frac{H\left( \mathbf {r}_{\alpha },t^{n+1} \right) -H\left( \mathbf {r}_{\alpha },t^{n} \right) }{{\Delta _t}} = \\&\sum _{\mathbf {r}_{\beta }} \Omega (\mathbf {r}_{\beta }, \mathbf {r}_{\alpha }) {\overline{\Lambda }\left( \mathbf {r}_{\beta },\mathbf {r}_{\alpha },t^n\right) } \frac{T\left( \mathbf {r}_{\beta } ,t^n\right) -T\left( \mathbf {r}_{\alpha } ,t^n\right) }{\left\| \mathbf {r}_{\beta } -\mathbf {r}_{\alpha } \right\| ^{2} } V_{\alpha \beta } \end{aligned} \end{aligned}$$

(24)

After determining the enthalpy at time step $t^{n+1}$, the temperature of the particle is calculated according to Eq. (3). Based on the new temperature value, the thermophysical properties of all particles are modified following Eqs. (1) and (2), and the calculation for the next time increment is performed. Equation (24) represents the numerical implementation of the developed bond-based PD model.

The discrete version of the radial heat flux (15) through a particle $\mathbf {r_{\alpha }}$ at a time instance $t^n$ is:

$$\begin{aligned} \begin{aligned} J^r (\mathbf {r_{\alpha }}, z, t^n) &= \sum _{\mathbf {r}_{\beta }^{+}}{\Omega (\mathbf {r}_{\beta }, \mathbf {r}_{\alpha }) } {\overline{\Lambda }\left( \mathbf {r}_{\beta },\mathbf {r}_{\alpha },t^n\right) \frac{T\left( \mathbf {r}_{\beta } ,t^n\right) -T\left( \mathbf {r}_{\alpha } ,t^n\right) }{\left\| \mathbf {r}_{\beta } -\mathbf {r}_{\alpha } \right\| } cos(\phi ) V_{\alpha \beta }} \end{aligned} \end{aligned}$$

(25)

where $\mathbf {r}_{\beta }^{+}$ is the set of particles within the horizon of the particle $\mathbf {r_{\alpha }}$ that has a higher temperature.

The discrete version of the axial heat flux (16) through a particle $\mathbf {r_{\alpha }}$ at a time instance $t^n$ is:

$$\begin{aligned} \begin{aligned}{ J^z (\mathbf {r_{\alpha }}, z, t^n) = \sum _{\mathbf {r}_{\beta }^{+}}}&{\Omega (\mathbf {r}_{\beta }, \mathbf {r}_{\alpha }) } {\overline{\Lambda }\left( \mathbf {r}_{\beta },\mathbf {r}_{\alpha },t^n\right) \frac{T\left( \mathbf {r}_{\beta } ,t^n\right) -T\left( \mathbf {r}_{\alpha } ,t^n\right) }{\left\| \mathbf {r}_{\beta } -\mathbf {r}_{\alpha } \right\| } cos(\psi ) V_{\alpha \beta }} \end{aligned} \end{aligned}$$

(26)

A simple numerical scheme to calculate the parameter $V_{\alpha \beta }$ for the case of regular grids was proposed in [11]; see also e.g. [15, 35, 36].

Details on implementing boundary conditions are presented in [37], where three options are analysed: mirror-type, naïve-type, and inner-type. For the first two types, following [10, 11], several layers of boundary particles are added around the simulated domain. The temperature of these particles is defined according to some law to ensure either constant temperature (Dirichlet boundary condition) or constant heat flux (Neumann boundary condition) at the domain boundary. Implementation of a convective (Robin) boundary condition in peridynamic diffusion models can be performed by a specific definition of the heat source term for the boundary particles. For more details about the implementation of different types of boundary conditions, see [17, 38,39,40].

In the next section, the developed PD formulations and their computational implementation will be verified using several test problems which have either known analytical solutions or can be solved with a well-established numerical method, the finite element method.

3 Verification

3.1 Axisymmetric Heat Conduction: a 1D Problem

To illustrate the accuracy of the developed PD axisymmetric formulation, the problem of cooling of a thick pipe, illustrated in Fig. 4, is considered. The outer diameter of the pipe is 3 m and the hole diameter is 0.1 m. The pipe material properties are as follows: thermal conductivity 50 W m${}^{-1}$ K${}^{-1}$, density 7850 kg m${}^{-3}$, and specific heat capacity 420 J kg${}^{-1}$ K${}^{-1}$. The pipe has an initial constant temperature and is cooled from the inner and outer boundary surfaces. If the length of the pipe is much larger than its diameter, the problem can be simplified to 1D axisymmetric. The simulated domain is then 1D of length 1.45 m as shown by the red line in the figure. The PD formulation uses particle size ${\Delta _\mathbf {r}}= 0.01$ m, horizon radius $\delta =0.03$ m, and time step ${\Delta _t}=5$ s.

Fig. 4

1D axisymmetric heat transfer problem. The red line is the 1D simulated domain

The initial pipe temperature is set to 500 K. Three types of boundary conditions are considered: constant temperature, constant heat flux, and constant convective heat flux. For the first type, the left and right boundaries of the domain are set at temperature 300 K. For the second type, constant heat fluxes of 10,000 W m${}^{-2}$ are prescribed in the outward directions. For the third type, the convective heat transfer coefficient is set to 500 W m${}^{-2}$ K${}^{-1}$ and the temperature of the liquid is set to 300 K. For these problems, the mirror-type boundary conditions on the ’left’ and ’right’ boundaries of the simulated domain are realized by the introduction of three layers of boundary particles. The temperatures at these particles are explicitly defined at every time step.

Simulations are performed for 10 h physical time.

The same problems were solved by the finite element method. The results for the time instances 2, 4, 6, 8, and 10 hours are presented in Fig. 5.

Fig. 5

Temperature distribution in the 1D axisymmetric domain for different types of boundary conditions: a constant temperature, b constant heat flux, c constant convective heat transfer coefficient. The solid lines, FEM results. The dotted lines, PD results

The comparison shows full agreement between the two numerical methods, providing confidence in the implementation of the PD approach for solving axisymmetric heat transfer problems under different types of boundary conditions.

3.2 Spherically Symmetric Heat Conduction: a 1D Problem

To assess the proposed specialisation to cases of spherical symmetry, the temperature distribution within a sphere is calculated and compared to the results of the FEM simulation. The sphere has a diameter of 0.2 m. The thermophysical properties of the material are the following: thermal conductivity 400 W m${}^{-1}$ K${}^{-1}$, density 8920 kg m${}^{-3}$, and specific heat capacity 385 J kg${}^{-1}$ K${}^{-1}$. The initial temperature of the sphere is 1000 K and the applied temperature at its outer surface is 300 K.

The simulated domain is 1D with a length 0.1 m. The boundary conditions are defined using three additional layers of boundary particles at the outer boundary, i.e., at one end of the 1D domain. No additional layers are introduced at the centre of the sphere, i.e., at the second end of the 1D domain. The particle size is ${\Delta _\mathbf {r}} = 0.001$ m, the horizon radius is $\delta = 0.003$ m, and the time step is ${\Delta _t}=0.001$ s. The total simulation time is 10 s. In this problem, the mirror-type boundary condition was set by the introduction of three layers of boundary particles at the ’right’ side of the 1D simulated domain that represents the outer surface of the sphere. The temperature at these boundary particles is defined at every time step. The ’left’ boundary of the simulated domain must have zero heat flux. As there is no space to allocate boundary particles, the zero heat flux can be enforced only by considering the inner-type boundary condition. However, it can be assumed that the zero heat flux boundary condition is realised by default, as there are no bonds that connect the domain with the environment. To simplify the model, this assumption is made.

The computational results based on the PD and FEM simulations are presented in Fig. 6.

Fig. 6

Temperature distribution in the 1D spherical domain. The solid line, FEM results. The dotted lines, PD results

The comparison shows an excellent agreement between the two numerical solutions. It should be noted that even without the presence of boundary particles at the centre of the sphere, there is no difference between the solutions. Therefore, the developed PD bond-based model can be used to describe the heat transfer in the presence of point (spherical) symmetry.

3.3 Axisymmetric Heat Transfer with Phase Change: a 1D Problem

There is no analytical solution for the 1D axisymmetric heat transfer problem with phase change and Dirichlet boundary conditions. However, there is a well-known solution for the case with constant heat flux, presented, e.g. in [34]. The problem is formulated for a semi-infinite domain, $0 \le x \le \infty$, with initial condition $T(x,0)=T_{\infty }$ and boundary conditions $q(0,t)=q_{0}$ (constant heat source) and $T(\infty ,t)=T_{\infty }$.

The following example considers the freezing of water in a domain similar to the one shown in Fig. 4. The outer radius of the domain is 1.1 m, and the inner radius is 0.1 m. The thermophysical properties of water are the following: thermal conductivity $\lambda _m = 0.6$ W m${}^{-1}$ K${}^{-1}$, $\lambda _f = 2.14$ W m${}^{-1}$ K${}^{-1}$; heat capacity $C_m = 4182$ J kg${}^{-1}$ K${}^{-1}$, $C_f = 2060$ J kg${}^{-1}$ K${}^{-1}$; density $\rho = 1000$ kg m${}^{-3}$; and latent heat of solidification $L = 333$ kJ/kg. The temperature of the phase change is $T_f = 273.15$ K. The initial and boundary conditions are $T\left( x,0\right) = T\left( \infty ,t\right) = 8$ K and $Q_h^r (0.1, 0, t) = 1000$W/m$^2$, i.e. heat is extracted from the inner surface with constant flux.

The simulations are performed with particle size ${\Delta _\mathbf {r}} = 0.01$ m, horizon radius $\delta =0.03$ m, and time step ${\Delta _t}=10$ s. In this problem, the mirror-type boundary conditions on the ’left’ and ’right’ boundaries of the simulated domain are realized by introducing three layers of boundary particles, where temperatures are defined explicitly at every time step. The process is followed for 3 days. The temperature distributions in the computational domain at time instances 1, 2, and 3 days based on numerical and analytical solutions are presented in Fig. 7.

Fig. 7

Temperature distribution in the 1D axisymmetric domain in case of phase change material. The solid line, FEM results. The dotted lines, PD results

The agreement between the PD and analytical results demonstrates that the proposed formulation captures accurately the axisymmetric heat transfer with phase change under constant heat flux boundary conditions. A small difference is observed for the particles with temperature equal to the phase change temperature. However, the difference can be reduced by decreasing the mesh size. A similar behaviour was observed with the solution of the phase change problem in [7], but it was less pronounced. It can be explained by the fact that in [7], the thermo-physical properties of the material were defined as continuous functions and the temperature form of equations was used, whereas in this study they are defined discontinuously and the enthalpy form is used.

3.4 Axisymmetric Heat Conduction: a 2D Problem

A thick pipe, illustrated in Fig. 8, with an initial temperature of 200 K is heated by a source of constant temperature attached to a part of the pipe’s outer surface as shown with a red line. The material of the pipe is the same as the one used in Sect. 3.1. This constitutes a 2D axisymmetric problem formulated on the square of size 0.8 $\times$ 0.8 m depicted in the figure. The parts of the surface that are not in contact with the heat source are insulated (zero heat flux boundary condition). The simulations are performed with particle size ${\Delta _\mathbf {r}} = 0.02$ m (the total number of domain particles is 1600); horizon size $\delta = 0.06$ m; and time step 10 s. In this problem, the mirror-type boundary conditions on all domain boundaries are realised by introducing three layers of boundary particles, where temperatures are defined explicitly at every time step. The simulated physical time is 10 h.

Fig. 8

Heating of a hollowed cylinder. $L = 0.2$ m. The thick blue line is the zero heat flux boundary condition. The thick red line is the constant temperature boundary condition

The same problem is solved by the finite element method. The temperature distributions in the domain after 5 h obtained by PD and by FEM simulations are presented in Fig. 9. Differences in the results between the two numerical methods cannot be observed.

Fig. 9

Temperature distribution in the cross-section of the hollowed cylinder at the time instance 5 h. The left panel is the result of PD simulation; the right panel is the result of FEM simulation

Comparisons of temperature distribution in the domain along the line $z = 0.4$ m obtained by PD and FEM simulations for several time instances are presented in Fig. 10a and the change of the enthalpy $\Delta H = H - H^0$, where $H^0$ is the initial enthalpy, is presented in 10b. The comparisons show excellent agreement between the two methods.

Fig. 10

Temperature distribution in the hollowed cylinder at the height of 0.4 m (a) and the change of enthalpy in this body (b). The solid line, FEM results. The dotted lines, PD results

To assess the accuracy of the developed PD implementation, the global numerical error, $\varepsilon$, can be calculated using the following relation [41]:

$$\begin{aligned} { \varepsilon = \frac{1}{\vert {T^{ref}}\vert _{\mathrm {max}}} \sqrt{\frac{1}{K} \sum _{n=1}^{K}\left( T_n^{ref}-T_n^{PD}\right) ^2}} \end{aligned}$$

(27)

where $T_n^{PD}$ and $T_n^{ref}$ are the temperatures at the position of particle n, obtained by the PD and by a reference method, respectively, K is the total number of PD particles, and $\vert {T^{ref}}\vert _{\mathrm {max}}$ is the maximum absolute reference temperature. Since an analytical solution for this problem does not exist, a FEM solution with a sufficiently fine mesh is used as a reference.

The solution accuracy depends on the particle size (or the distance between the neighbouring particles) $\Delta _\mathbf {r}$, and the relative size of the PD horizon $m=\delta /\Delta _\mathbf {r}$. To demonstrate the influence of these factors, results obtained with $\Delta _\mathbf {r} = 0.04, 0.02, 0.01 \mathrm {m}$ and with m from 2 to 8 are presented in Fig. 11a for the time instance 5 h. The area of the PD horizon decreases with decreasing particle size $\Delta _\mathbf {r}$; hence, it is appropriate to present this dataset as a function of the ratio between the area of a single particle horizon, $\pi \delta ^2$, and the area of the simulation domain, $A_{dom}$. The results are presented in Fig. 11b in a logarithmic scale.

Fig. 11

The global numerical error of the model for the time instance 5 h

As can be seen from Fig. 11, the PD solution converges to the reference solution with decreasing particle size $\Delta _\mathbf {r}$. Notably, the solution for $\Delta _\mathbf {r} = 0.01 \text{ m}$ and $m=2$ breaks the general convergence trend. This solution underestimates the temperature at the PD particles, whereas all other cases overestimate it. It can be inferred that for the case with $\Delta _\mathbf {r} = 0.01 \text{ m}$, the numerical error $\varepsilon$ becomes zero somewhere between $m=2$ and $m=3$, and then increases. This is supported by the results for $m=2.5$ shown in Fig. 11. The overall accuracy of the PD approach is acceptable, as the global error does not exceed 0.5 % for reasonable values of $\Delta _\mathbf {r}$ and m. It should be noted that the increase of m decreases the computational speed, which is crucial for problems with phase change because the matrix of coefficients is updated after every time step.

3.5 Axisymmetric Heat Transfer with Phase Change: a 2D Problem

In this subsection, the problem of cooling a ’glass’ of boiling water with a cylindrical piece of ice is analysed. The geometry of the problem is presented in Fig. 12. The ’glass’ is a cylinder with a diameter 0.2 m and a height 0.2 m. A cylindrical block of ice with a diameter 0.8 m and a height 0.8 m is placed in hot water. The material properties are the same as those given in Sect. 3.3. The initial temperature of the water is 373.15 K, and the temperature of the ice block is 268.15 K. In this problem, the effect of convection is neglected; therefore, heat transfer is only by conduction.

Fig. 12

Cooling of water by a block of ice. The red line is the simulated domain

As the body has two types of symmetry—axial symmetry and plane symmetry—only a square domain with size 0.1 $\times$ 0.1 shown in the figure is simulated. The simulations are performed with particle size ${\Delta _\mathbf {r}= 0.01}$ m (the total number of domain particles is 10,000); horizon size $\delta = {0.03}$ m; and time time step 2 s. The simulated physical time is 120 min.

The mirror-type boundary conditions on the ’bottom’, ’top’, and ’right’ domain boundaries are realized by introducing three layers of boundary particles, where temperatures are defined explicitly at every time step. At the ’left’ domain boundary along the z-axis there is no space to add layers of additional boundary particles; it is assumed that the absence of bonds that cross this axis represents a zero flux boundary condition. Therefore, a slight decrease in the accuracy of the results should be expected, which was previously discussed in [11].

The temperature fields after 20 and 60 min obtained by PD and FEM simulations are presented in Fig. 13. The results demonstrate that the agreement between the two numerical solutions is very good.

Fig. 13

Temperature distribution in the ’glass’ of water with a piece of ice at the time instances 20 and 60 min. The left panels are the results of PD simulation; the right panels are the results of FEM simulation. The red line is the interface between liquid and solid phases with temperature $T_f = 273.15 K$

The temperature distribution along the horizontal plane of symmetry is presented in Fig. 14a, and along the vertical axis of symmetry in Fig. 14b. The agreement between PD and FEM results is very good. It should be noted that for the particle adjacent to the z-axis, their horizons include only half of the number of necessary particles, which can lead to some inaccuracy if the particle size is large. However, for the presented case, no tangible difference is observed.

Fig. 14

Temperature distribution in the ’glass’ of water with ’a piece of ice’: a the temperature at the horizontal plane of symmetry ($z = 0.002$m); b the temperature at the vertical axe of symmetry ($r = 0.001$m). The solid line, FEM results. The dotted lines, PD results

4 Conclusions

A bond-based peridynamic formulation of heat transfer with phase change in domains with axial and spherical symmetry is presented and verified. This allows for simplifying real 3D problems to corresponding 2D problems in the case of axial symmetry, and to 1D problems in the case of spherical symmetry with corresponding significant reductions of PD particles used and computational times required. The presented non-local formulation can be used to develop mathematical models of physical problems that involve strong non-linearities and discontinuities. Such problems are challenging for the classical local formulations that use partial differential equations to describe the fundamental conservation laws in physics.

The formulation uses enthalpy as the main variable and discontinuous functions for properties of bonds. The specialisation to symmetrical problems is achieved by introducing a shape factor $\Omega (\mathbf {r}', \mathbf {r})$ that modifies the heat transfer in bonds depending on the position of their adjacent particles with respect to the axis or centre of symmetry.

The accuracy of the proposed formulation and its numerical implementation is verified using several test problems that either have analytical solutions or could be solved with a well-established and trusted numerical method. The verification examples demonstrated that the proposed PD approach for problems with axial and spherical symmetry can simulate accurately heat transfer and heat transfer with phase change problems.