Static condensation of peridynamic heat conduction model

Abstract

A model order reduction methodology for reducing the order of the peridynamic transient heat model is proposed. This methodology is based on the static condensation procedure. To set the development of the model reduction procedure on a sound mathematical setting, a nonlocal vector calculus was employed in the formulation of the heat transport problem. The model order reduction framework proposed in this study provides a technique to reduce the dimensionality of a peridynamic transport model while still maintaining accurate prediction of the model response. Moreover, the methodology can be adaptively applied to accommodate different resolution requirements for different sections of the model. Using numerical experiments, the proposed methodology is shown to be capable of accurately reproducing results of the full peridynamic transient heat transport problem.

1. Introduction

The peridynamic model for solid mechanics was proposed by Silling [1] to provide a suitable alternative to the classical continuum model where the later fails to provide meaningful results, either because the classical model fall short of capturing the underlying physics such as the existence of long-range or nonlocal interaction or where the partial derivatives in the classical model become unevaluable due to the presence of singularities such as cracks, shockwaves, and interfaces. The existence of long-range interaction in dynamic fracture was suggested in Ramulu et al. [2] in which a crack curving criterion was proposed as a sufficient condition for crack branching. This criterion is predicated on the postulate that the direction of crack propagation is dictated by the microcracks ahead of the crack tip. Under this postulate, a crack will deflect if a microcrack within a critical distance $r_c$ to the crack tip is actuated by a critical crack tip stress field. The distance $r_c$ thus introduces a notion of a long-range or nonlocal interaction. Similarity between the $r_c$ and the peridynamic horizon has been suggested in Bobaru and Hu [3].

Because differential operators are the building block of the mathematical framework on which the classical continuum theory is based, modelling of systems with singular surfaces where large changes occur in some field properties of the system becomes challenging or impossible. This is because, differential operators as composed of derivatives are undefined on such singularities. To make models based on the classical continuum theory applicable in situations where such singularities exist, various remedial methodologies have been proposed within the framework of the classical continuum theory. These include the Linear Elastic Fracture Mechanics [4–6] and methodologies based on the cohesive zone model [7–10].

Classical continuum-based methods have been used to model heat transport in bodies with singularities. In Sih [11], complex variable method was used to find solution to the problem of heat conduction in a medium with line discontinuities. The heat flux obtained from the classical continuum model exhibited the characteristic inverse square root singularity as is observed with stress near singularities in linear elastic fracture mechanics. The finite element method has been applied in modelling heat transport in fractured medium [12–16]. A common characteristic shared by all these methodologies is the requirement to either modify the underlying equations or to provide additional criteria for special treatment of the fractured surface. Another problem associated with application of finite element method in modelling heat transport in medium containing multiple random cracks is non-triviality of mesh generation [17].

Despite the advancement recorded in the use of methods based on the classical continuum theory in modelling physical phenomena involving singularities such as cracks, these methodologies are still unable to accurately predict a wide range of problems [18]. Peridynamics (PD) provide a framework that admits discontinuity and long-range forces. A key feature of the mathematical framework of peridynamic theory that gives it this robustness is the use of integral operators instead of spatial differential operators. Owing to this mathematical framework and its theoretical structure, the peridynamic theory accounts for nonlocality and allows for the modelling of a wide range of physical phenomena in the realm of mechanics and heat transport in materials with discontinuous or evolving topology.

Peridynamics has been successfully applied in studying a wide spectrum of physical phenomena in the realm of solid and fluid mechanics. These include studies in dynamic crack propagation [19–35], fatigue analysis [36–40], vibration and wave dispersion [41–46]. In the field of heat transport, considerable research efforts have been expended in developing and applying peridynamic framework to model heat transport phenomena. The first attempt at applying peridynamic framework to model transport phenomena appeared in Gerstle et al. [47]. The developed framework is a multi-physics framework that allowed for the coupling of thermo-mechanical, electrical potential distribution and atomic diffusion mechanisms in a single model. Further development of the peridynamic electromigration appeared in previous works [48,49]. A Bond-based peridynamic transient thermal conduction model based on Fourier heat conduction law was developed in Bobaru and Duangpanya [50] and was later extended for application on bodies with evolving topology in Bobaru and Duangpanya [51].

A state-based thermal diffusion model was proposed in Oterkus et al. [52]. The formulation utilized the Lagrange’s equation to derive an ordinary state-based peridynamic heat diffusion equation. The formulation allows for the utilization of the heat flux from the classical heat conduction model within the nonlocal peridynamic framework. A deviation from deriving the peridynamic heat conduction equation based on the Fourier law appeared in Wang et al. [53] where the nonlocal peridynamic heat conduction equation was derived based on the Dual-Phase-Lag model [54].

However, the increased accuracy and modelling capability offered by peridynamics is partially offset by several numerical and modelling challenges: some of which include increased computational cost due to the nonlocality of the peridynamic model [55] and nontrivial prescription of boundary conditions [56–58]. Another challenge associated with peridynamics is selection of appropriate model parameters such as the so-called peridynamic horizon [3]. Our interest in this paper is focused on developing a framework that will reduce the computational cost of simulating the nonlocal peridynamic transient heat conduction model. Efforts have been made to implement peridynamic model with fewer degrees of freedom to reduce computational cost. A coarsening method was proposed in Silling [59] which allows accurate solution of the linearized peridynamic model to be obtained using fewer degrees of freedom. This framework was later extended for two-dimensional implementation in Galadima et al. [60]. The coarsening methodology has the limitation of not being suitable for reducing dynamic models. A model order reduction for peridynamics based on static condensation was proposed in Galadima et al. [61] like the coarsening method but with extended capability of reducing the order of dynamic models. All the frameworks above have been applied only to solid mechanics problems. This research work aim at developing a model order reduction framework for peridynamic transport model.

The model order reduction technique proposed here is achieved by applying static condensation to the peridynamic heat transport model. Although the focus in this work is on heat conduction, the framework developed here can easily be applied in modelling other transport phenomena such as mass transport, owing to the common mathematical framework shared by these phenomena. Because our goal is to achieve the model order reduction of a nonlocal peridynamic heat transport model, we will first introduce some element of nonlocal calculus that will provide a concise and systematic formulation of the heat transport problem. This review will be presented in section 1.1 and is based on and will largely follow the exposition in Du et al. [62]. This presentation will review only results that are deemed relevant to this work. The key result to be presented are the nonlocal divergence and gradient operators. These operators in conjunction with appropriate thermal conduction coefficients serve as the building blocks in the development of the nonlocal thermal conduction model in section 1.2. In section 2, a framework for static condensation of the transient peridynamic heat conduction model will be proposed. The framework developed will be applied to numerical problems in section 3 to be followed by discussion of results and concluding remarks.

1.1. Nonlocal calculus

The exposition about nonlocal vector calculus will start with the definition of some key concepts. We will begin with the definition of functions and operators. Let $k,m$ and $n$ be positive integers and let $\mathit{x}$ and $\mathit{y}$ be points in ${\mathbb{R}}^n$ . For a given domain $\mathrm{\Omega }\subseteq {\mathbb{R}}^n$ , functions or operators that maps $\mathrm{\Omega }$ into ${\mathbb{R}}^{m\times n}$ or ${\mathbb{R}}^n$ or $\mathbb{R}$ are called point functions or operators, respectively. On the other hand, functions, or operators from $\mathrm{\Omega }\times \mathrm{\Omega }$ into ${\mathbb{R}}^{m\times n}$ or ${\mathbb{R}}^n$ or $\mathbb{R}$ are called two-point functions or operators, respectively. Point functions and two-point functions could be scalar, vector, or tensor valued functions. Since one of the key components of a balance law particularly as it relates to transport phenomena is the concept of flux, it is therefore necessary to introduce next in this discussion the notion of a nonlocal flux. Given a scalar two-point function $\psi :{\mathbb{R}}^n\times {\mathbb{R}}^n\to \mathbb{R}$ , then the definition

$$\mathcal{F}({\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2;q):=\underset{{\mathrm{\Omega }}_1}{\int }\underset{{\mathrm{\Omega }}_2}{\int }\psi (\mathit{x},\mathit{y})d\mathit{y}d\mathit{x}$$

(1)

is the nonlocal flux of $q$ from ${\mathrm{\Omega }}_1$ into ${\mathrm{\Omega }}_2$ where $\underset{{\mathrm{\Omega }}_2}{\int }\psi (\mathit{x},\mathit{y})d\mathit{y}$ is identified as the nonlocal flux density into the region ${\mathrm{\Omega }}_2$ from point $\mathit{x}\in {\mathrm{\Omega }}_1$ . The flux is nonlocal because the operator may be nonzero even if ${\overline{\mathrm{\Omega }}}_1$ and ${\overline{\mathrm{\Omega }}}_2$ have an empty intersection. If $\psi (x,y)$ is assumed to be antisymmetric, then the following statements are true:

• 1. There is no self-interaction, i.e.,

$$\underset{\mathrm{\Omega }}{\int }\underset{\mathrm{\Omega }}{\int }\psi (\mathit{x},\mathit{y})d\mathit{y}d\mathit{x}=0.$$

(2)

• 2. The nonlocal action-reaction principle holds for ${\mathrm{\Omega }}_1$ , ${\mathrm{\Omega }}_2\subset \mathrm{\Omega }$

$$\underset{{\mathrm{\Omega }}_1}{\int }\underset{{\mathrm{\Omega }}_2}{\int }\psi (\mathit{x},\mathit{y})d\mathit{y}d\mathit{x}+\underset{{\mathrm{\Omega }}_2}{\int }\underset{{\mathrm{\Omega }}_1}{\int }\psi (\mathit{x},\mathit{y})d\mathit{y}d\mathit{x}=0.$$

(3)

The action–reaction principle given by equation (3) simply states that the flux from ${\mathrm{\Omega }}_1$ into ${\mathrm{\Omega }}_2$ is equal and opposite to the flux that exits ${\mathrm{\Omega }}_2$ into ${\mathrm{\Omega }}_1$

1.1.1. Nonlocal divergence and gradient operators and their adjoint

Given a vector two-point function $\mathit{v}:{\mathbb{R}}^n\times {\mathbb{R}}^n\to {\mathbb{R}}^k$ , and a scalar point function $u\left(\mathit{x}\right):{\mathbb{R}}^n\to \mathbb{R}$ . Let $\mathbf{\alpha }(\mathit{x},\mathit{y}):{\mathbb{R}}^n\times {\mathbb{R}}^n\to {\mathbb{R}}^k$ be an antisymmetric vector two-point function. The action of nonlocal divergence operator $\mathcal{D}$ and its adjoint ${\mathcal{D}}^{*}$ on $\mathit{v}$ and $u$ , respectively, are defined as

$$\mathcal{D}\left(\mathit{v}\right)\left(\mathit{x}\right):=\underset{{\mathbb{R}}^n}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\alpha (\mathit{x},\mathit{y}))d\mathit{y}\forall \mathit{x}\in {\mathbb{R}}^n$$

(4)

and

$${\mathcal{D}}^{*}\left(\eta \right)(\mathit{x},\mathit{y}):=-(u\left(\mathit{y}\right)-u\left(\mathit{x}\right))\alpha (\mathit{x},\mathit{y})\forall \mathit{x}\in {\mathbb{R}}^n,$$

(5)

where $\mathcal{D}\left(\mathit{v}\right)\left(\mathit{x}\right):{\mathbb{R}}^n\to \mathbb{R}$ and ${\mathcal{D}}^{*}\left(\eta \right)(\mathit{x},\mathit{y}):{\mathbb{R}}^n\times {\mathbb{R}}^n\to {\mathbb{R}}^k$ .

Given a scalar two-point function $\eta :{\mathbb{R}}^n\times {\mathbb{R}}^n\to \mathbb{R}$ and a vector point function $\mathbf{v}\left(\mathit{x}\right):{\mathbb{R}}^n\to {\mathbb{R}}^k$ . Let $\mathit{\beta }(\mathit{x},\mathit{y}):{\mathbb{R}}^n\times {\mathbb{R}}^n\to {\mathbb{R}}^k$ be an antisymmetric vector two-point function. Then the action of the nonlocal gradient operator $\mathcal{G}$ and its adjoint ${\mathcal{G}}^{*}$ on $\eta$ and $\mathbf{v}$ , respectively, are defined as

$$\mathcal{G}\left(\eta \right)\left(\mathit{x}\right):=\underset{{\mathbb{R}}^n}{\int }(\eta (\mathit{x},\mathit{y})+\eta (\mathit{y},\mathit{x})\mathit{\beta }(\mathit{x},\mathit{y}))d\mathit{y}\forall \mathit{x},\mathit{y}\in {\mathbb{R}}^n$$

(6)

and

$${\mathcal{G}}^{*}\left(\eta \right)(\mathit{x},\mathit{y})=-(\mathbf{v}\left(\mathit{y}\right)-\mathbf{v}\left(\mathit{x}\right))·\mathit{\beta }(\mathit{x},\mathit{y})\forall \mathit{x},\mathit{y}\in {\mathbb{R}}^n.$$

(7)

Observe that, unlike in local calculus which deals with point functions only, nonlocal calculus involves two kinds of functions: point and two-point functions. This therefore necessitates the definition of alternative forms of the nonlocal operators defined in equations (4)–(7). The alternative forms of the nonlocal divergence and gradient operators were given in Galadima et al. [61] to be the pairs $\mathcal{D},-{\mathcal{G}}^{*}$ and $\mathcal{G},-{\mathcal{D}}^{*}$ . For example, in the alternative pair $\mathcal{D},-{\mathcal{G}}^{*}$ , while $\mathcal{D}$ operates on tensor two-point functions, $-{\mathcal{G}}^{*}$ operates on tensor point functions. Similarly, for the pair $\mathcal{G},-{\mathcal{D}}^{*}$ , $\mathcal{G}$ operates on vector two-point functions while $-{\mathcal{D}}^{*}$ operate on vector point functions.

1.1.2. Interaction kernels and domains

In equations (4)–(7), the two-point vector functions $\mathbf{\alpha }(\mathit{x},\mathit{y})$ and $\mathbf{\beta }(\mathit{x},\mathit{y})$ are also known as the interaction kernels. In the broad context of Peridynamics, these interaction kernels are assumed to have a finite domain that does not map to zero.

Given two points $\mathit{x},\mathit{y}\in {\mathbb{R}}^n$ , and $\delta \in {\mathbb{R}}^{+}$ . Let ${\mathcal{B}}_{\delta }(\mathit{x})$ be a ball or radius $\delta$ centred at $\mathit{x}$ , then, for example, the interaction kernel $\alpha (\mathit{x},\mathit{y})$ is nonzero only if $\mathit{y}\in {\mathcal{B}}_{\delta }(\mathit{x})$ , i.e.

$$\alpha (\mathit{x},\mathit{y})\{\begin{array}{c}\ne 0\forall y\in {\mathcal{B}}_{\delta }(\mathit{x})\hfill \\ \multicolumn{1}{c}{=0\forall y\notin {\mathcal{B}}_{\delta }(\mathit{x})}\end{array}$$

(8)

where $\delta$ is the interaction radius also known as the horizon in the context of PD. Interaction kernels that satisfy equation (8) are called truncated kernels [63] or localized kernels [64]. Another key concept that is connected to the notion of truncated kernels is the interaction domain. Let ${\mathrm{\Omega }}_{\mathrm{s}}\subset {\mathbb{R}}^n$ be a bounded open set. The interaction domain ${\mathrm{\Omega }}_I$ consists of points outside of ${\mathrm{\Omega }}_{\mathrm{s}}$ that interact with points in ${\mathrm{\Omega }}_{\mathrm{s}}$ . Given the truncated interaction kernel $\mathbf{\alpha }(\mathit{x},\mathit{y})$ , an interaction domain ${\mathrm{\Omega }}_I$ is defined as

$${\mathrm{\Omega }}_I=\{\mathit{y}\in {\mathbb{R}}^n\backslash {\mathrm{\Omega }}_s:\exists \mathit{x}\in {\mathrm{\Omega }}_s:\mathbf{\alpha }(\mathit{x},\mathit{y})\ne 0\}.$$

(9)

The interaction subdomain contains all the points $\mathit{y}$ in the complement domain ${\mathbb{R}}^n\backslash {\mathrm{\Omega }}_s$ that interact with points $\mathit{x}$ in ${\mathrm{\Omega }}_{\mathrm{s}}$ . Many geometrical relationships exist between ${\mathrm{\Omega }}_{\mathrm{s}}$ and ${\mathrm{\Omega }}_I$ [62]. A typical such relationship is shown in Figure 1.

Figure 1.

A typical relationship between $\mathrm{\Omega }$ and ${\mathrm{\Omega }}_I$ .

1.1.3. Nonlocal interaction operators

Given a domain ${\mathrm{\Omega }}_{\mathrm{s}}$ , let ${\mathrm{\Omega }}_I$ be the interaction domain associated with ${\mathrm{\Omega }}_{\mathrm{s}}$ as defined in section 1. Corresponding to the nonlocal divergence operator $\mathcal{D}\left(\mathit{v}\right)\left(\mathit{x}\right)$ , a point interaction operator $\mathcal{N}\left(\mathit{v}\right):{\mathrm{\Omega }}_I\to \mathbb{R}$ is defined as

$$\mathcal{N}\left(\mathit{v}\right)\left(\mathit{x}\right):=-\underset{{\mathrm{\Omega }}_s\cup {\mathrm{\Omega }}_I}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}\forall \mathit{x}\in {\mathrm{\Omega }}_I.$$

(10)

Corresponding to the nonlocal gradient operator $\mathcal{G}\left(\eta \right)\left(\mathit{x}\right)$ , a point interaction operator $\mathcal{S}\left(\eta \right)\left(\mathit{x}\right):{\mathrm{\Omega }}_I\to {\mathbb{R}}^k$ is defined as

$$\mathcal{S}\left(\eta \right)\left(\mathit{x}\right):=-\underset{{\mathrm{\Omega }}_s\cup {\mathrm{\Omega }}_I}{\int }(\mathbf{\eta }(\mathit{x},\mathit{y})+\mathbf{\eta }(\mathit{y},\mathit{x})\mathbf{\beta }(\mathit{x},\mathit{y}))d\mathit{y}\forall \mathit{x}\in {\mathrm{\Omega }}_I.$$

(11)

1.1.4. Nonlocal integral theorem

A very important outcome of the nonlocal operators developed in the proceeding subsections is statement of the nonlocal Gauss theorem. Recall from equation (2) that

$$\begin{array}{c}\underset{\mathrm{\Omega }}{\int }\underset{\mathrm{\Omega }}{\int }\psi (\mathit{x},\mathit{y})d\mathit{y}d\mathit{x}=0\hfill \\ \underset{{\mathrm{\Omega }}_{\mathrm{s}}\cup {\mathrm{\Omega }}_I}{\int }\underset{{\mathrm{\Omega }}_{\mathrm{s}}\cup {\mathrm{\Omega }}_I}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}=0\hfill \\ \underset{{\mathrm{\Omega }}_{\mathrm{s}}}{\int }\underset{{\mathrm{\Omega }}_{\mathrm{s}}\cup {\mathrm{\Omega }}_I}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}+\underset{{\mathrm{\Omega }}_I}{\int }\underset{{\mathrm{\Omega }}_{\mathrm{s}}\cup {\mathrm{\Omega }}_I}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}=0\hfill \\ \underset{{\mathrm{\Omega }}_s}{\int }\mathcal{D}\left(\mathit{v}\right)\left(\mathit{x}\right)d\mathit{x}=-\underset{{\mathrm{\Omega }}_I}{\int }\underset{{\mathrm{\Omega }}_s\cup {\mathrm{\Omega }}_I}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}.\hfill \end{array}$$

(12)

$$\underset{{\mathrm{\Omega }}_s}{\int }\mathcal{D}\left(\mathit{v}\right)\left(\mathit{x}\right)d\mathit{x}=\underset{{\mathrm{\Omega }}_I}{\int }\mathcal{N}\left(\mathit{v}\right)(\mathit{x})d\mathit{x}.$$

(13)

Consider the right-hand side of equation (12). If the interaction kernel function $\mathbf{\alpha }(\mathit{x},\mathit{y})$ is antisymmetric and equation (2) holds, then

$$\begin{array}{c}-\underset{{\mathrm{\Omega }}_I}{\int }\underset{\mathrm{\Omega }\cup {\mathrm{\Omega }}_I}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}=-\underset{{\mathrm{\Omega }}_I}{\int }\underset{\mathrm{\Omega }}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}\hfill \\ =\underset{\mathrm{\Omega }}{\int }\underset{{\mathrm{\Omega }}_I}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}.\hfill \end{array}$$

(14)

From equation (14) and considering equation (13), it can be deduced that $\underset{{\mathrm{\Omega }}_I}{\int }\mathcal{N}\left(\mathit{v}\right)(\mathit{x})d\mathit{x}$ in equation (13) represents the flux from ${\mathrm{\Omega }}_{\mathrm{s}}$ into ${\mathrm{\Omega }}_I$ . Thus, equation (13) is the mathematical statement of the nonlocal Gauss theorem which postulate that the integral of the nonlocal divergence of $\mathit{v}$ over ${\mathrm{\Omega }}_s$ is equal to the total flux exiting ${\mathrm{\Omega }}_s$ into ${\mathrm{\Omega }}_I$ .

1.1.5. Nonlocal flux

From the statement of the divergence theorem which can be stated as the flux of a vector $\mathit{v}\in {\mathbb{R}}^n$ out of a region ${\mathrm{\Omega }}_1\subset {\mathbb{R}}^n$ is equal to the volume integral of the divergence of $\mathit{v}$ over ${\mathrm{\Omega }}_1$ . Thus, from equations (1) and (4),

$$\mathcal{F}({\mathrm{\Omega }}_1,{\mathrm{\Omega }}_2;q)=\underset{{\mathrm{\Omega }}_1}{\int }\mathcal{D}\left(\mathit{v}\right)\left(\mathit{x}\right)d\mathit{x}=\underset{{\mathrm{\Omega }}_1}{\int }\underset{{\mathrm{\Omega }}_2}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}.$$

(15)

1.2. Peridynamic model

The goal of the peridynamic theory is to provide a nonlocal alternative to the classical continuum theory that can unify the treatment of both continuous and discontinuous response of bodies in a single mathematical framework. Our goal in this section is to derive the nonlocal peridynamic transient heat conduction model using the principle of the nonlocal calculus presented in section 1.1. To proceed with our objective, it is necessary to restate some important concept in peridynamic theory. Given a bounded open domain ${\mathrm{\Omega }}_{\mathrm{s}}\in {\mathbb{R}}^n$ , in PD, a continuum point $\mathit{x}\in \mathrm{\Omega }$ interacts with infinitely many other points located within its domain of influence. If this domain of influence is assumed to be a ball ${\mathcal{B}}_{\delta }(\mathit{x})$ of radius $\delta >0$ centred at $\mathit{x}$ , then $\delta$ is called the horizon of $\mathit{x}$ , such that

$${\mathcal{B}}_{\delta }\left(\mathit{x}\right)=\{{\mathit{x}}^{\prime }\in \mathcal{R}:|{\mathit{x}}^{\prime }-\mathit{x}|<\delta \},$$

(16)

where ${\mathcal{B}}_{\delta }(\mathit{x})$ is called the family of $\mathit{x}$ . Interaction between two points $\mathit{x}$ and ${\mathit{x}}^{\prime }$ is called a bond and the distance $\xi ={\mathit{x}}^{\prime }-\mathit{x}$ in the undeformed reference configuration is called the bond vector.

In the original development [1] that came to be known as the Bond-Based Peridynamics, interaction within a bond is independent of all other bonds. This lack of dependence imposes certain limitation to the peridynamic model. For example, in the realm of solid mechanics, this restricted the value of Poisson’s ratio to 1/3 for 2D and 1/4 for 3D isotropic solids. To circumvent this limitation and obtain a more general material model, State Based Peridynamics was developed. To pursue the development of the state based peridynamics, mathematical objects called states which are functions defined on bonds in ${\mathcal{B}}_{\delta }\left(\mathit{x}\right)$ were introduced in [65]. To define the domain $\mathcal{H}$ of a state, let $\delta >0$ be the horizon of a point $\mathit{x}$ in a body $\mathfrak{B}$ . Then,

$$\mathcal{H}=\{\mathbf{\xi }\in (\mathbb{R}\backslash 0)|(\mathbf{\xi }+\mathit{x})\in ({\mathcal{B}}_{\delta }\left(\mathit{x}\right)\cap \mathfrak{B})\}$$

(17)

is the family of bonds for the point $\mathit{x}$ .

1.2.1. The nonlocal peridynamic balance law

The governing equation of motion in PD can be formulated using the nonlocal vector calculus presented in section 1.1 based on a statement of balance law which postulate the dependence of the rate of change in the content of an extensive quantity over a given domain on the rate at which the quantity is produced within the domain and a flux through the boundary of the domain. Let the region occupied by a body $\mathcal{B}$ be given by the open domain $\mathrm{\Omega }\subseteq {\mathbb{R}}^n$ . Let $\tilde{\mathrm{\Omega }}\subseteq \mathrm{\Omega }$ be an open subregion, then a quantitative statement of a balance law for $\mathrm{\Omega }$ can be stated as

$$\mathcal{A}(\tilde{\mathrm{\Omega }};q)=\mathcal{P}\left(\tilde{\mathrm{\Omega }}\right)-\mathcal{F}(\tilde{\mathrm{\Omega }},{\mathbb{R}}^n\backslash \tilde{\mathrm{\Omega }};q),$$

(18)

where equation (18) postulate that $\mathcal{A}(\tilde{\mathrm{\Omega }},q)$ (the time rate of change of the intensive quantity $q$ ) is equal to $\mathcal{P}\left(\tilde{\mathrm{\Omega }}\right)$ (the rate at which the quantity is produced within the subdomain by sources, minus $\mathcal{F}(\tilde{\mathrm{\Omega }},{\mathbb{R}}^n\backslash \tilde{\mathrm{\Omega }};q)$ (the rate at which the intensive quantity exits the subdomain). The rate of change and production of $q$ within $\tilde{\mathrm{\Omega }}$ can be written as

$$\mathcal{A}(\tilde{\mathrm{\Omega }};q)=\frac{\partial }{\partial t}\underset{\tilde{\mathrm{\Omega }}}{\int }q(\mathit{x},t)d\mathit{x}\mathrm{and}\mathcal{P}\left(\tilde{\mathrm{\Omega }}\right)=\underset{\tilde{\mathrm{\Omega }}}{\int }b(\mathit{x},t)d\mathit{x}.$$

(19)

If $\tilde{\mathrm{\Omega }}\subset \mathrm{\Omega }$ is assumed to interact only with the open subregion ${\tilde{\mathrm{\Omega }}}_I\subset {\mathrm{\Omega }}_I$ , then equation (18) can be written as

$$\mathcal{A}(\tilde{\mathrm{\Omega }};q)=\mathcal{P}\left(\tilde{\mathrm{\Omega }}\right)-\mathcal{F}(\tilde{\mathrm{\Omega }},{\tilde{\mathrm{\Omega }}}_I;q).$$

(20)

If we write ${\mathrm{\Omega }}_1=\tilde{\mathrm{\Omega }}$ and ${\mathrm{\Omega }}_2={\tilde{\mathrm{\Omega }}}_I$ , then equation (15) becomes

$$\begin{array}{c}\mathcal{F}(\tilde{\mathrm{\Omega }},{\tilde{\mathrm{\Omega }}}_I;q)=\underset{\tilde{\mathrm{\Omega }}}{\int }\underset{{\tilde{\mathrm{\Omega }}}_I}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}\hfill \\ =-\underset{{\tilde{\mathrm{\Omega }}}_I}{\int }\underset{\tilde{\mathrm{\Omega }}}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}\hfill \\ =-\underset{{\tilde{\mathrm{\Omega }}}_I}{\int }\underset{\tilde{\mathrm{\Omega }}\cup {\tilde{\mathrm{\Omega }}}_I}{\int }(\mathit{v}(\mathit{x},\mathit{y})+\mathit{v}(\mathit{y},\mathit{x})·\mathbf{\alpha }(\mathit{x},\mathit{y}))d\mathit{y}d\mathit{x}\hfill \\ \mathcal{F}(\tilde{\mathrm{\Omega }},{\tilde{\mathrm{\Omega }}}_I;q)=\underset{{\tilde{\mathrm{\Omega }}}_I}{\int }\mathcal{N}\left(\mathit{v}\right)\left(\mathit{x}\right)d\mathit{x}\forall \tilde{\mathrm{\Omega }}\subset \mathrm{\Omega },\hfill \end{array}$$

(21)

where the third equality follows from consideration of equation (2) and the fourth is due to equation (10). The balance law equation (18) can now be written as

$$\frac{\partial }{\partial t}\underset{\tilde{\mathrm{\Omega }}}{\int }q(\mathit{x},t)d\mathit{x}+\underset{{\tilde{\mathrm{\Omega }}}_I}{\int }\mathcal{N}\left(\mathit{v}\right)\left(\mathit{x}\right)d\mathit{x}=\underset{\tilde{\mathrm{\Omega }}}{\int }b(\mathit{x},t)d\mathit{x}\forall \tilde{\mathrm{\Omega }}\subset \mathrm{\Omega }.$$

(22)

Applying the nonlocal divergence theorem equation (13) yields

$$\underset{\tilde{\mathrm{\Omega }}}{\int }(\frac{\partial q(\mathit{x},t)}{\partial t}+\mathcal{D}\left(\mathit{v}\right)-b(\mathit{x},t))d\mathit{x}=0\forall \tilde{\mathrm{\Omega }}\subset \mathrm{\Omega },$$

(23)

which owing to the arbitrariness of $\tilde{\mathrm{\Omega }}$ , equation (23) localizes to the field equation

$$\frac{\partial q(\mathit{x},t)}{\partial t}+\mathcal{D}\left(\mathit{v}\right)-b(\mathit{x},t)=0\forall \mathit{x}\in \mathrm{\Omega }.$$

(24)

To uniquely determine $q$ , equation (24) is required to satisfy some constraint conditions. Let the interaction domain ${\mathrm{\Omega }}_I$ be split into two disjoint subdomains ${\mathrm{\Omega }}_{\mathrm{Id}}$ and ${\mathrm{\Omega }}_{\mathrm{In}}$ such that ${\mathrm{\Omega }}_{\mathrm{Id}}\cap {\mathrm{\Omega }}_{\mathrm{In}}=Ø$ and either of ${\mathrm{\Omega }}_{\mathrm{Id}}$ and ${\mathrm{\Omega }}_{\mathrm{In}}$ could be an empty set. ${\mathrm{\Omega }}_{\mathrm{Id}}$ is the subdomain where Dirichlet boundary condition is applied and ${\mathrm{\Omega }}_{\mathrm{In}}$ is the subdomain where Neumann boundary condition is applied. Analogous to the boundary value problem of the classical local theory, constraints on the solution $q$ of equation (24) over $\mathrm{\Omega }$ is applied as follows: a given function value ${\mathit{g}}_d$ is prescribed on the solution over ${\mathrm{\Omega }}_{\mathrm{Id}}$ such that

$$q\left(\mathit{x}\right)=g_d\forall \mathit{x}\in {\mathrm{\Omega }}_{\mathrm{Id}}.$$

(25)

To prescribe the Neumann type constraint, recall that in the classical boundary value problem, this involves prescribing a flux density $\mathit{q}·\mathit{n}$ over the flux boundary. From equation (14) and the discussion that follows, the nonlocal flux density over ${\mathrm{\Omega }}_{\mathrm{In}}$ is given by $\underset{{\mathrm{\Omega }}_I}{\int }\mathcal{N}\left(q\right)\left(\mathit{x}\right)d\mathit{x}$ . Let $g_n$ be a given function value of the flux density over ${\mathrm{\Omega }}_{\mathrm{In}}$ . The Neumann constraint can be stated as

$$\underset{{\mathrm{\Omega }}_{\mathrm{In}}}{\int }\mathcal{N}\left(q\right)\left(\mathit{x}\right)d\mathit{x}=g_{\mathit{n}}\forall \mathit{x}\in {\mathrm{\Omega }}_{\mathrm{In}}.$$

(26)

The presence of the first order time derivative of the solution $q$ in equation (24) means in addition to the boundary constraints equation (25), equation (26), initial conditions also need to be specified. The initial condition involves prescribing the initial values of the solution and its first derivative. Let $q_I$ and ${\stackrel{·}{q}}_I$ be the initial values of $q(\mathit{x})$ and $\stackrel{·}{q}\left(\mathit{x}\right)$ , respectively, then

$$q(\mathit{x},0)=q_I\forall \mathit{x}\in \mathrm{\Omega },\mathrm{for}t=0$$

(27)

and

$$\stackrel{·}{q}(\mathit{x},0)={\stackrel{·}{q}}_I\forall \mathit{x}\in \mathrm{\Omega },\mathrm{for}t=0$$

(28)

are the initial conditions. So that equations (24)–(28) give the complete definition of the nonlocal problem which can be stated as

$$\{\begin{array}{cc}\mathrm{\rho }\stackrel{·}{q}=-\mathcal{D}\left(q\right)\left(\mathit{x}\right)+b(\mathit{x},t)\hfill & \forall \mathit{x}\in {\mathrm{\Omega }}_s\hfill \\ \multicolumn{1}{c}{q\left(\mathit{x}\right)=g_d}& \forall \mathit{x}\in {\mathrm{\Omega }}_{\mathrm{Id}}\hfill \\ \multicolumn{1}{c}{\underset{{\mathrm{\Omega }}_{\mathrm{In}}}{\int }\mathcal{N}\left(q\right)\left(\mathit{x}\right)d\mathit{x}=g_{\mathit{n}}}& \forall \mathit{x}\in {\mathrm{\Omega }}_{\mathrm{Id}}\hfill \\ \multicolumn{1}{c}{q(\mathit{x},0)=q_I}& \forall \mathit{x}\in \mathrm{\Omega },\mathrm{for}t=0\hfill \\ \multicolumn{1}{c}{\stackrel{·}{q}(\mathit{x},0)={\stackrel{·}{q}}_I}& \forall \mathit{x}\in \mathrm{\Omega },\mathrm{for}t=0\hfill \end{array}.$$

(29)

For the specific case of heat conduction, $q$ is the stored internal energy density expressed as

$$q=\mathrm{\rho }{C}_v(\mathrm{\tau }(\mathit{x},t)-{\mathrm{\tau }}^{*}\left(\mathit{x}\right)),$$

(30)

where $\mathrm{\rho }$ and $C_v$ are the density and specific heat capacity of the material of the body, respectively. $\mathrm{\tau }$ is the temperature of the current state and ${\mathrm{\tau }}^{*}$ is the temperature of a reference state. Then the rate of stored energy is given by

$$\frac{\partial q(\mathit{x},t)}{\partial t}=\mathrm{\rho }{C}_v\frac{\partial \mathrm{\tau }(\mathit{x},t)}{\partial t}.$$

(31)

The stored internal energy density $q$ is related to the interaction vector $\mathit{v}$ through a constitutive law. Let the constitutive law be such that

$$\mathit{v}=\frac{1}{2}\mathcal{D}(\mathbf{\Theta }·{\mathcal{D}}^{*}\left(\mathrm{\tau }\right)),$$

(32)

where ${\mathcal{D}}^{*}\left(\mathrm{\tau }\right)=\mathcal{G}\left(\mathrm{\tau }\right)$ is the nonlocal gradient of the temperature field (see section 1.1.1). $\mathbf{\Theta }$ is a second order tensor given by $\mathbf{\Theta }=k\mathbf{I}$ , where $k$ is a constant and $\mathbf{I}$ is the identity tensor. Substituting equations (31) and (32) into equation (24) and dropping the time and spatial dependence of variables for the sake of brevity and write $\mathrm{\tau }(\mathit{y},t)={\mathrm{\tau }}^{\prime }$ and $\mathrm{\tau }(\mathit{x},t)=\mathrm{\tau }$ , we obtain

$$\begin{array}{c}\frac{\partial \mathrm{\tau }}{\partial t}=\frac{1}{2\mathrm{\rho }{C}_v}[-\mathcal{D}(\mathrm{\Theta }·{\mathcal{D}}^{*}\left(\mathrm{\tau }\right))+b]\forall \mathit{x}\in \mathrm{\Omega }\hfill \\ =\underset{\mathrm{\Omega }}{\int }\frac{k}{\mathrm{\rho }{C}_v}(({\mathrm{\tau }}^{\prime }+\mathrm{\tau })\mathit{B}·\mathbf{\alpha }-({\mathrm{\tau }}^{\prime }+\mathrm{\tau })\mathbf{\alpha }·{\mathbf{\alpha }}^{\prime })\mathrm{dx}+\frac{k}{\mathrm{\rho }{C}_v}b\hfill \\ =\underset{\mathrm{\Omega }}{\int }\frac{k}{\mathrm{\rho }{C}_v}(({\mathrm{\tau }}^{\prime }+\mathrm{\tau })\mathbf{\alpha }·\mathbf{\alpha }+({\mathrm{\tau }}^{\prime }+\mathrm{\tau })\mathbf{\alpha }·\mathbf{\alpha })\mathrm{dx}+\frac{k}{\mathrm{\rho }{C}_v}b\hfill \\ \frac{\partial \mathrm{\tau }}{\partial t}=\underset{\mathrm{\Omega }}{\int }\frac{k}{\mathrm{\rho }{C}_v}\mathbf{\alpha }·\mathbf{\alpha }({\mathrm{\tau }}^{\prime }-\mathrm{\tau })d\mathit{x}+\frac{k}{\mathrm{\rho }{C}_v}b\forall \mathit{x}\in \mathrm{\Omega }.\hfill \end{array}$$

(33)

If we write $\mathbf{\gamma }=\mathbf{\alpha }·\mathbf{\alpha }$ , we obtain

$$\frac{\partial \mathrm{\tau }}{\partial t}=\underset{\mathrm{\Omega }}{\int }\frac{k}{\mathrm{\rho }{C}_v}\gamma ({\mathrm{\tau }}^{\prime }-\mathrm{\tau })d\mathit{x}+\frac{k}{\mathrm{\rho }{C}_v}b\forall \mathit{x}\in \mathrm{\Omega }.$$

(34)

The integrand of the integral in equations (33) and (34) without the factor $1/(\mathrm{\rho }{C}_v)$ is the thermal response function, and $k$ is a peridynamic parameter called the micro-conductivity of bonds. The general form [47, 66] of the micro-conductivity is given by

$$k(\xi ,\delta )=k(0,\delta )g(\xi ,\delta ),$$

(35)

where $k(0,\delta )=k_0$ is a constant and $g(\xi ,\delta )$ is a kernel function and gives the spatial distribution of the intensity of long-range thermal interaction between material points in a material. The explicit form of the micro-conductivity will be given in section 1.2.2. Analogous to thermal diffusivity of the classical heat conduction model, we define a micro-diffusivity function as

$$K=\frac{k}{\mathrm{\rho }{C}_v}\gamma .$$

(36)

Just like the micro-conductivity parameter, the micro-diffusivity function is also a parameter intrinsic to the material. Thus equation (34) can be written as

$$\frac{\partial \mathrm{\tau }}{\partial t}=\underset{\mathrm{\Omega }}{\int }K({\mathrm{\tau }}^{\prime }-\mathrm{\tau })d\mathit{x}+h\forall \mathit{x}\in \mathrm{\Omega },$$

(37)

where $h=\frac{b}{\mathrm{\rho }{C}_v}$ is a heat source density. The explicit form of the micro-diffusivity function and hence the thermal response function depends on the choice of the interaction kernel function. In what follows, some candidate kernel functions will be utilized to derive the explicit forms of the thermal response functions that appeared in the literature. It is important to state that the analysis and classification or systematic justification of the interaction kernels is not an objective of this work, rather, the objective is to demonstrate the robustness of the nonlocal operators by showing how for appropriate choice of the interaction kernel, the nonlocal operators are able to yield specific form of the Peridynamic heat equation. Our approach is to a posteriori select those kernels that will result in the form of thermal response functions that are of interest to this study.

If we define the kernel function as $\mathbf{\alpha }(\mathit{x},\mathit{y})=\frac{\mathbf{\xi }}{{\left|\mathbf{\xi }\right|}^2}$ , then diffusivity modulus (equation (36)) becomes

$$K=\frac{k}{\mathrm{\rho }{C}_v{\left|\mathbf{\xi }\right|}^2},$$

(38)

which leads to a form of the response function that appeared in Bobaru and Duangpanya [50] herein denoted as $f_1$ given by

$$f_1=\frac{k}{{\left|\mathbf{\xi }\right|}^2}({\mathrm{\tau }}^{\prime }-\mathrm{\tau }).$$

(39)

If we, however, define $\mathbf{\alpha }(\mathit{x},\mathit{y})=(\mathbf{\xi }/\left|\mathbf{\xi }\right|)$ , then from equations (33) and (37), equation (36) becomes

$$K=\frac{k}{\mathrm{\rho }{C}_v},$$

(40)

resulting in a response function that was proposed in Gerstle et al. [47] herein denoted as $f_2$ and is expressed as

$$f_2=k({\mathrm{\tau }}^{\prime }-\mathrm{\tau }).$$

(41)

On the other hand, if we define the kernel in equation (34) as $\gamma =\frac{1}{|y-x|}$ , we obtain from equation (36)

$$K=\frac{k}{\begin{array}{c}\mathrm{\rho }{C}_v\left|\xi \right|\hfill \end{array}},$$

(42)

which corresponds to the form of the response function in previous works [62,67] given as

$$f_3=\frac{k}{\left|\mathbf{\xi }\right|}({\mathrm{\tau }}^{\prime }-\mathrm{\tau }).$$

(43)

1.2.2. Relationship between micro-conductivity and thermal conductivity

The peridynamic micro-conductivity $k$ is related to the thermal conductivity $\kappa$ from the classical heat conduction model using the principle of constitutive correspondence [50,52]. In Oterkus et al. [52], this relationship is established by assuming equivalence of thermal potential between the peridynamic and classical continuum models. To show the explicit form of this relationship, the following definitions are made based on Agwai [67]: The peridynamic thermal potential is given as

$$Z(\mathit{x},t)=\frac{1}{2}\underset{\mathcal{H}}{\int }z(\mathit{y},\mathit{x},t)d{V}_{\mathit{y}},$$

(44)

where $z$ is the microthermal potential. The microthermal potential is related to the normalized response function through

$$f=\frac{\partial z}{\partial \theta },$$

(45)

where $\theta ={\mathrm{\tau }}^{\prime }-\mathrm{\tau }$ . The form of the microthermal potential corresponding to different form of the response functions are obtained by substituting either of equations (39), (41), or (43) into equation (45) to, respectively, obtain

$$z_1=k\frac{{\mathrm{\tau }}^2}{2{\left|\mathbf{\xi }\right|}^2},z_2=k\frac{{\mathrm{\tau }}^2}{2},z_3=k\frac{{\mathrm{\tau }}^2}{2\left|\mathbf{\xi }\right|}.$$

(46)

The thermal potential from the classical model is given as

$$\hat{Z}=\frac{1}{2}(\nabla \mathrm{\tau }·\kappa \nabla \mathrm{\tau }),$$

(47)

where $\kappa$ is the thermal conductivity of the medium.

1.2.2.1. One-dimensional micro-conductivity

If a linear temperature field of the form $\mathrm{\tau }\left(\mathit{x}\right)=c\left|\mathit{x}\right|$ is applied such that $c$ is a constant and $\left|\mathit{x}\right|$ is the length of the coordinate vector $\mathit{x}$ , then equation (46) becomes

$$z_1=k\frac{{\mathbf{\xi }}^2}{2{\left|\mathbf{\xi }\right|}^2},z_2=k\frac{{\mathbf{\xi }}^2}{2},z_3=k\frac{{\mathbf{\xi }}^2}{2\left|\mathbf{\xi }\right|}$$

(48)

and equation (47) becomes

$$\hat{Z}=\frac{1}{2}\kappa c^2.$$

(49)

Let $\delta$ be the horizon of material points in the medium, substituting equation (48) into equation (44) assuming one-dimensional problem, yields

$$Z_1(\mathit{x},t)=\frac{1}{2}\underset{-\delta }{\overset{\delta }{\int }}k\frac{c^2{\mathbf{\xi }}^2}{2{\left|\mathbf{\xi }\right|}^2}\mathrm{Ad}\xi ,Z_2(\mathit{x},t)=\frac{1}{2}\underset{-\delta }{\overset{\delta }{\int }}k\frac{c^2{\mathbf{\xi }}^2}{2}\mathrm{Ad}\xi ,Z_3(\mathit{x},t)=\frac{1}{2}\underset{-\delta }{\overset{\delta }{\int }}k\frac{c^2{\mathbf{\xi }}^2}{2\left|\mathbf{\xi }\right|}\mathrm{Ad}\xi ,$$

(50)

where $A$ is the cross-sectional area associated with the point $x$ . If we neglect the effect of distance between material points on their interaction and assume a micro-conductivity function of the form $k(x,y)=k_0$ , then equation (50), respectively, evaluates to

$$\begin{array}{c}{Z}_1(\mathit{x},t)=\frac{1}{2}\underset{-\delta }{\overset{\delta }{\int }}{k}_0\frac{c^2{\mathbf{\xi }}^2}{2{\left|\mathbf{\xi }\right|}^2}\mathrm{Ad}\xi =\frac{k_0{c}^2\delta A}{2}\hfill \\ \multicolumn{1}{c}{Z_2(\mathit{x},t)=\frac{1}{2}\underset{-\delta }{\overset{\delta }{\int }}{k}_0\frac{c^2{\mathbf{\xi }}^2}{2}\mathrm{Ad}\xi =\frac{k_0{c}^2{\delta }^{3}A}{6}}\\ \multicolumn{1}{c}{Z_3(\mathit{x},t)=\frac{1}{2}\underset{-\delta }{\overset{\delta }{\int }}{k}_0\frac{c^2{\mathbf{\xi }}^2}{2\left|\mathbf{\xi }\right|}\mathrm{Ad}\xi =\frac{k_0{c}^2{\delta }^{2}A}{4}.}\end{array}$$

(51)

Equating (49) and (51), respectively, yields

$$k_1=\frac{\kappa }{A\delta },k_2=\frac{3\kappa }{A{\delta }^3},k_3=\frac{2\kappa }{A{\delta }^2}.$$

(52)

In the case where the distance between points is assumed to have considerable effect on the long-range thermal interaction between them, then the kernel function chosen for the micro-conductivity function must be able to capture the distribution of intensity of interaction. Typical kernel functions used in this case are the linear and quadratic kernel functions given, respectively, as $g(\left|\mathbf{\xi }\right|,\delta )=(1-(\left|\mathbf{\xi }\right|/\delta ))$ and $g(\left|\mathbf{\xi }\right|,\delta )=(1-(\mathbf{\xi }/\delta {)}^2{)}^2$ so that the micro-conductivity functions takes the form $k(\mathit{x},\mathit{y})=k_0(1-(\left|\xi \right|/\delta ))$ and $k(\mathit{x},\mathit{y})=k_0(1-(\mathbf{\xi }/\delta {)}^2{)}^2$ , respectively. For example, substituting these into the first of equation (50), working out the resulting integral and equating thermal potentials from peridynamics with that from the classical theory allow us to obtain expressions for the micro-conductivity corresponding to the different forms of the kernel functions, respectively, as

$$k=\frac{2\kappa }{A\delta }(1-\frac{\left|\xi \right|}{\delta })$$

(53)

and

$$k=\frac{15\kappa }{8A\delta }{(1-{\left(\frac{\xi }{\delta }\right)}^2)}^2.$$

(54)

The micro-conductivity functions (equations (53)–(54)) all have compact support within the horizon. In other words, the interaction between a given point and any other point in a body cease if the distance between them is greater the horizon $(\delta )$ . The typical forms of different micro-conductivity functions are presented in Figure 2 which represents the graph of $k_3$ in equation (52) as well as equations (53) and (54).

Figure 2.

Comparison of micro-conductivity functions.

Although the different forms of the thermal response functions and their associated micro-conductivity functions have been established in the literature as cited, results of numerical experiments [67] have shown that the response function (equation (43)) provides the best approximation to the classical solution and would henceforth be used in the development of the proposed model reduction methodology. Also, the three different forms of the micro-conductivity functions presented in this section will be statically condensed in section 2.2. However, the choice of the micro-conductivity used in section 3 is limited to the constant function because of its simplicity of implementation and because the condensation of the peridynamic heat model with the other linear and quadratic micro-conductivity functions follows the same procedure.

1.2.2.2. Two-dimensional micro-conductivity function

Following the same procedure laid out in section 1.2.2.1, the explicit form of the constant and triangular micro-conductivity as well as the corresponding diffusivity modulus function can be obtained. A constant micro modulus function with its corresponding diffusivity function are given as

$$k=\frac{6\kappa }{\pi h{\delta }^3},K=\frac{6\kappa }{\pi h{\delta }^3\left|\mathbf{\xi }\right|}.$$

(55)

Interested readers are referred to Agwai [67] for discussion on the various forms of the response functions and their corresponding micro-conductivity functions.

1.2.3. The discrete heat conduction equation

Different discretization schemes have been proposed for the numerical approximation of balance laws arising from the peridynamic theory, such as the meshfree method [68,69], the collocation methods [70,71], and methods based on finite elements mesh [55,72]. Due to its simple implementation algorithm and relatively low computational cost, the meshfree method suggested in Silling and Askari [68] is the most widely used [73] and is the preferred method in this work for these same reasons. In this approximation method, the discrete form of equation (37) is

$${\stackrel{·}{\mathrm{\tau }}}_i=\sum _i^{N}K({\mathrm{\tau }}_j-{\mathrm{\tau }}_i)V_j+h_i,$$

(56)

where ${\mathrm{\rho }}_i=\mathrm{\rho }({\mathit{x}}_i)$ , $\stackrel{·}{\mathrm{\tau }}=\frac{\partial {\mathrm{\tau }}_i}{\partial t}$ , and $N$ is the number of nodes in the neighbourhood of node $i$ .

2. Static condensation of the peridynamic head conduction model

The assembled peridynamic transient heat conduction equations for the body in matrix form is given by

$$\mathit{I}\stackrel{·}{\mathrm{\tau }}+\mathit{K}\mathrm{\tau }=h,$$

(57)

where $\mathit{I},\mathit{K}\in {\mathbb{R}}^{n\times n}$ are the identity and diffusivity modulus matrices, respectively, and $n$ is the total number of degrees of freedom in the system. In this context, equation (57) is referred to as the full order model of the system or simply full model. The objective of the model reduction is to replace the full model having $n$ degree of freedom with a reduced order system having $m$ degree of freedom, such that $m<<n$ . The reduced model is expected to preserve the heat conduction characteristics of the full model. To proceed with the model condensation, the degrees of freedom of the system are separated into retained and truncated degrees of freedom. The retained degrees of freedom are those to be preserved while truncated degrees of freedom are those to be condensed out in the reduced model. Let $r$ denote the retained degrees of freedom and let $t$ denote the truncated degrees of freedom, then equation (57) can be partitioned as follows:

$$\left[\begin{array}{cc}{\mathit{I}}_{\mathrm{rr}}\hfill & {\mathit{I}}_{\mathrm{rt}}\hfill \\ \multicolumn{1}{c}{{\mathit{I}}_{\mathrm{tr}}}& {\mathit{I}}_{\mathrm{tt}}\hfill \end{array}\right]\left\{\begin{array}{c}{\stackrel{·}{\mathrm{\tau }}}_r\hfill \\ \multicolumn{1}{c}{{\stackrel{·}{\mathrm{\tau }}}_t}\end{array}\right\}+\left[\begin{array}{cc}{\mathit{K}}_{\mathrm{rr}}\hfill & {\mathit{K}}_{\mathrm{rt}}\hfill \\ \multicolumn{1}{c}{{\mathit{K}}_{\mathrm{tr}}}& {\mathit{K}}_{\mathrm{tt}}\hfill \end{array}\right]\left\{\begin{array}{c}{\mathrm{\tau }}_r\hfill \\ \multicolumn{1}{c}{{\mathrm{\tau }}_t}\end{array}\right\}=\left\{\begin{array}{c}{h}_r\hfill \\ \multicolumn{1}{c}{h_t}\end{array}\right\}.$$

(58)

Neglecting the transient term and assuming there is no heat source at the truncated degrees of freedom, then the solution of the second submatrix equation for ${\mathrm{\tau }}_t$ yields

$${\mathrm{\tau }}_t={\mathit{R}}_G{\mathrm{\tau }}_r,$$

(59)

where ${\mathit{R}}_G\in {\mathbb{R}}^{t\times r}$ is the condensation matrix, defined as

$${\mathit{R}}_G=-{\mathit{K}}_{\mathrm{tt}}^{-1}{\mathit{K}}_{\mathrm{tr}}$$

(60)

The condensation matrix relates the retained degrees of freedom and truncated degrees of freedom and is load independent because it is assumed there is no heat source at the truncated degrees of freedom. Also note that in arriving at equation (59), the transient term in equation (58) has been neglected; hence, this method is called static condensation method. The temperature state of the full model may be expressed as

$$\mathrm{\tau }=\left\{\begin{array}{c}{\mathrm{\tau }}_r\hfill \\ \multicolumn{1}{c}{{\mathrm{\tau }}_t}\end{array}\right\}={\mathit{T}}_G{\mathrm{\tau }}_r,$$

(61)

where ${\mathit{T}}_G\in {\mathbb{R}}^{n\times r}$ is a transformation matrix given by

$${\mathit{T}}_G=\left[\begin{array}{c}\mathit{I}\hfill \\ \multicolumn{1}{c}{{\mathit{R}}_G}\end{array}\right],$$

(62)

where $\mathit{I}$ is $r\times r$ identity matrix. If the transformation matrix ${\mathit{T}}_G$ is assumed to be independent of time, differentiating equation (60) with respect to time gives

$$\stackrel{·}{\mathrm{\tau }}={\mathit{T}}_G{\stackrel{·}{\mathrm{\tau }}}_r.$$

(63)

Substituting equations (62) and (63) into equation (57) and pre-multiplying by ${\mathit{T}}_G^T$ yields

$${\mathit{I}}_G{\stackrel{·}{\mathrm{\tau }}}_r+{\mathit{K}}_G{\mathrm{\tau }}_r=h_G,$$

(64)

where ${\mathit{I}}_G$ , ${\mathit{K}}_G$ , and $h_G$ are the reduced identity matrix, reduced diffusivity matrix, and reduced load vector, respectively, associated with the reduced model, defined as

$${\mathit{I}}_G={\mathit{T}}_G^T\mathit{I}{\mathit{T}}_G,{\mathit{K}}_G={\mathit{T}}_G^T\mathit{K}{\mathit{T}}_G,h_G={\mathit{T}}_G^{T}h.$$

(65)

2.1. Handling boundary conditions

It is necessary to provide some notes on how the Dirichlet and Neumann boundary conditions are handled in this model order reduction scheme. When implementing the static condensation procedure, the region constrained by a Dirichlet type boundary condition must be retained in the reduced model. This is because in Dirichlet boundary condition, a specified value of the temperature is imposed on the boundary volume. As can be seen from equation (8), there is no transformation on the temperature field. Hence, to preserve the dynamics of the system, any boundary with Dirichlet condition imposed must be retained. This is not the case with Neumann boundary condition. This is because in Neuman type boundary condition, instead of temperature, a specified value of heat is applied to the boundary volume. From equation (8) and the definitions in equation (9), in going from the full model to a reduced model, there is a transformation relation between the heat source term in the full model and its counterpart in the reduced model (see the third definition in equation (9)). Therefore, irrespective of what region is retained in the reduced model, the characteristics of the full model are preserved. This affords more flexibility on the final configuration of the reduced model compared with the Dirichlet boundary condition scenario.

2.2. Condensation of parameters of peridynamic heat conduction model

The condensation of the micro-conductivity and micro-diffusivity functions will be illustrated using an example of a bar as shown in Figure 3. The goal is to demonstrate the behaviour of the peridynamic parameters when subjected to static condensation. Let the bar be of length $1$ unit. Three micro-conductivity functions and their corresponding micro-diffusivity functions will be studied: the constant ( $k_3$ in equation (52)), linear (equation (53)), and quadratic (equation (54)) micro-conductivity functions.

Figure 3.

A discretized bar to illustrate condensation of the micro-conductivity functions.

The horizon is $\delta =0.4$ and $\mathrm{\rho },C_v=1$ . The bar is discretized into nodes with spacing $\mathrm{dx}=0.01$ . Condensation is carried out by retaining every fourth node in the full model as shown in Figure 3. To proceed with the condensation, matrices of the micro-conductivity $\mathit{k}$ and the diffusivity $\mathit{K}$ of the bar is computed. These matrices are computed using the expression for $k_3$ in equation (52) as well as equations (53) and (54) in the case of micro-conductivity while the diffusivity matrices corresponding to these micro-conductivity functions are obtained by substituting the micro-conductivity functions into equation (43). Condensation is achieved by introducing $\mathit{k}$ or $\mathit{K}$ into the second expression of equation (65) to obtain the condensed micro-conductivity matrix ${\mathit{k}}_G$ or the condensed diffusivity matrix ${\mathit{K}}_G$ .

The condensed curves of the micro-conductivity functions and the micro-diffusivity functions are characterized by sharp spikes as can be seen from Figures 4 and 5, respectively. This behaviour is expected since the condensed model is defined only at the retained degrees of freedom.

Figure 4.

Static condensation of micro-conductivity functions of a homogeneous bar: (a) constant micro-conductivity, (b) linear micro-conductivity, and (c) quadratic micro-conductivity.

Figure 5.

Static condensation of micro-diffusivity functions of a homogeneous bar: (a) corresponding to constant micro-conductivity, (b) corresponding to linear micro-conductivity, and (c) corresponding to quadratic micro-conductivity.

3. Numerical verification

The static condensation scheme described in section 2 will be tested on numerical problems to verify and demonstrate its capabilities in reducing the order of a peridynamic transient heat conduction model. To achieve this goal, the temperature field in a one-dimensional bar and two-dimensional plate is predicted using the proposed static condensation scheme and the result is compared with prediction using the full model.

3.1. A homogeneous bar with constant temperature applied at both ends

A homogeneous bar initially at zero temperature is subjected to boundary temperature of 1°C at both ends as shown in Figure 6. The bar has a length and thickness of 1 m and 0.01 m, respectively. The specific heat capacity, thermal conductivity, and density of the bar are specified as $C_v=64\mathrm{J}/\mathrm{kgK}$ , $\kappa =233\mathrm{W}/\mathrm{mK}$ , and $\mathrm{\rho }=260\mathrm{kg}/{\mathrm{m}}^3$ , respectively. The material of the bar is assumed to have a micro-conductivity function of the form given in equation (48).

Figure 6.

A homogeneous bar subjected to constant temperature at both ends.

The bar is discretized into $100$ nodes; thus, the spacing between nodes is $\Delta x=0.01$ . Three instances of the model reduction algorithm would be considered in this example. In the first instance, model order reduction of the bar is achieved by retaining every second node of the full model. In the second instance, reduction is achieved by retaining every fifth node in the full model, and finally in the third instance reduction is achieved by retaining every 10th node in the full model. A time step size of $\Delta t={10}^{-2}\mathrm{s}$ is used and the simulation was run for $3000$ time steps.

Temperature prediction for both full model and reduced model for the three instances of reduction are presented in Figures 7–9. The reported results are temperature distribution across the bar corresponding to time steps $t=10\mathrm{s},20\mathrm{s},\mathrm{and}30\mathrm{s}$ . These results show a good correlation between the full model and reduce model irrespective of the reduction algorithm employed, thus demonstrating the capability of the model order reduction scheme in preserving the dynamics of a one-dimensional transient peridynamic heat conduction model in-spite of using fewer degrees of freedom.

Figure 7.

Temperature distribution in a homogeneous bar subjected to constant temperature at both ends corresponding to retaining every second node of the full model and response function (a) $f_1$ , (b) $f_2$ , and (c) $f_3$ .

Figure 8.

Temperature distribution in a homogeneous bar subjected to constant temperature at both ends corresponding to retaining every fifth node of the full model and response function (a) $f_1$ , (b) $f_2$ , and (c) $f_3$ .

Figure 9.

Temperature distribution in a homogeneous bar subjected to constant temperature at both ends corresponding to retaining every 10th node of the full model and response function (a) $f_1$ , (b) $f_2$ , and (c) $f_3$ .

3.2. Numerical study of heat conduction in plate with a pre-existing crack

The static condensation scheme will be used to study temperature distribution in a plate with a pre-existing crack. The geometry of the plate is specified as length $L=0.1\mathrm{m}$ , width $W=0.1\mathrm{m}$ , thickness $H=3.3\times {10}^{-4}\mathrm{m}$ , and the crack length $a=0.0016\mathrm{m}$ as shown in Figure 10. The thermal conductivity, specific heat capacity, and density are specified, respectively, as: $\kappa =233\mathrm{W}/\mathrm{mK}$ , $C_v=64\mathrm{J}/\mathrm{kgK}$ , and $\mathrm{\rho }=260\mathrm{kg}/{\mathrm{m}}^3$ .

Figure 10.

Example 2 problem setup: A plate with pre-existing crack.

The numerical implementation of the peridynamic transient heat conduction model proceed by discretizing the plate into $90$ nodes along each edge, so that we have $8100$ nodes. The distance between nodes is $\Delta x=\Delta y=1.1\times {10}^{-4}\mathrm{m}$ . The horizon of the full model is $\delta =3.015\Delta x$ . Model order reduction on the plate is carried out as shown schematically in Figure 11. The condensation algorithm adopted is to retain every second node of every second row. The total simulation time is $t={10}^{-3}s$ with a time step of $\Delta t={10}^{-6}\mathrm{s}$ . The plate would be studied for two case scenarios: Dirichlet and Neumann boundary conditions.

Figure 11.

Schematic representation of static condensation.

3.2.1. Example 1. Case I: Dirichlet boundary condition

In this case, the plate is subjected to the following initial and boundary conditions, respectively:

$$T((x,y),0)=0^{o}C,-L/2\le x\le \frac{L}{2},-W/2\le y\le W/2,$$

(66)

and

$$T((x,y),t)=1^{o}C,-\frac{L}{2}\le x\le \frac{L}{2},y=-W/2.$$

(67)

Note that in this case, although the condensation algorithm may exclude the boundary nodes, we are constrained to still retain them. Temperature distribution in time across the plate for both full model and the reduced model are shown Figure 12. A profile of the temperature along a grid of nodes parallel to the vertical axis of the plate at $x=-0.000055\mathrm{m}$ is also shown in Figure 13.

Figure 12.

Dirichlet boundary condition: Temperature distribution across the plate: (a) Full model corresponding to $t=7\times {10}^{-5}\mathrm{s},$ (b) Full model corresponding to $t=7\times {10}^{-4}\mathrm{s}$ , (c) Full model corresponding to $t=7\times {10}^{-3}\mathrm{s}$ , (d) Reduced model corresponding to $t=7\times {10}^{-5}\mathrm{s},$ (e) Reduced model corresponding to $t=7\times {10}^{-4}\mathrm{s},$ , (f) Reduced model corresponding to $t=7\times {10}^{-3}\mathrm{s}$ ,.

Figure 13.

Temperature profile along the grid at $x=-0.000055\mathrm{m}$ parallel to the y-axis of the plate with Dirichlet boundary condition.

The results presented in Figures 12 and 13 show a very good match between predictions from the full model and those from the reduced model. The difference between the full model and reduce model curves close to the crack location in Figure 13 reflects wider spacing between nodes in the reduced model due to fewer nodes. Note that all boundary nodes subjected to the Dirichlet boundary condition were retained in the reduced model.

3.2.2. Example 2. Case II: Neumann boundary condition

The goal in this example is to demonstrate the performance of the model order reduction in reducing a 2D model with a Neumann type boundary condition. In this case, in addition to equation (66), the plate is subjected to the following conditions:

$$h((\mathit{x},\mathit{y}),t)={10}^7,0\le x\le L,y=0.$$

(68)

Figure 14(a)–(c) illustrates the full model prediction of the distribution of temperature over the plate, while Figure 14(d)–(f) are the predicted temperature distribution in the plate from the reduced model corresponding to simulation times $t=0.00007\mathrm{s}$ , $t=0.0007\mathrm{s}$ and $t=0.007\mathrm{s}$ , respectively. Figure 15 presents the curves of temperature profile along a grid parallel to the vertical axis of the plate for both full and reduced model corresponding to simulation times $t=0.0007\mathrm{s}$ , $t=0.0035\mathrm{s},$ and $t=0.007\mathrm{s}$ , respectively.

Figure 14.

Neumann boundary condition: Temperature distribution across the plate: (a) Full model corresponding to $t=7\times {10}^{-5}\mathrm{s},$ (b) Full model corresponding to $t=7\times {10}^{-4}\mathrm{s}$ (c) Full model corresponding to $t=7\times {10}^{-3}\mathrm{s}$ (d) Reduced model corresponding to $t=7\times {10}^{-5}\mathrm{s},$ , (e) Reduced model corresponding to $t=7\times {10}^{-4}\mathrm{s},$ , (f) Reduced model corresponding to $t=7\times {10}^{-3}\mathrm{s}$ ,.

Figure 15.

Temperature profile along the grid at $x=-0.000055\mathrm{m}$ parallel to the y-axis of the plate with Neumann boundary condition.

The results presented show a near identical match between predictions from the full model and those from the reduced model. It is worthy to note that in the case of the Neumann type boundary condition, there was no constraint to the choice of nodes to retain and those to be condensed out. This potentially allows for the use of fewer degrees of freedom than in the scenario where a Dirichlet type boundary condition is imposed.

3.2.3. Example 3. Neumann boundary condition with retained nodes selected randomly

In this example, the plate of example 2 will be analysed using the same boundary conditions and plate properties. However, in coarsening the plate, instead of following a predetermined pattern as done in the previous examples, the retained nodes in this case will be selected randomly. The goal is to demonstrate the robustness of the condensation methodology in handling any kind of selection pattern during model reduction. In this condensation scheme, 2,025 nodes were randomly selected from a total of 8,190 nodes. Result of the temperature distribution over the surface of the plate is shown in Figure 16, while the temperature profile of nodes that falls along a grid located at $x=-0.000055\mathrm{m}$ .

Figure 16.

Neumann boundary condition: Temperature distribution across the plate: (a) Full model corresponding tot = 7×10^(-5) s, (b) Full model corresponding to t = 7×10^(-4) s (c) Full model corresponding to t = 7×10^(-3) s,(d) Reduced model corresponding.

The results presented in Figures 16 and 17 show good agreement between prediction from the full model and the condensed model. This is a demonstration of the robustness of the condensation methodology in predicting accurate result irrespective of the algorithm or pattern employed in the condensation process.

Figure 17.

Temperature profile along a grid of the plate with Neumann boundary condition and condensation achieved through random selection of nodes.

3.2.4. Example 3. Nonuniform condensation

It is sometimes necessary to have high resolution at certain region of the model while such level of detail may not be required in other regions. This model order reduction can be selectively applied to accommodate this different resolution requirement. For example, high resolution may be required around regions close to the crack of the plate in section 3.2.2 to allow for more detailed information than is required in regions further away from the cracks. To achieve this, a nonuniform condensation algorithm will be implemented in this example where all nodes around the crack are retained and regions further away from the crack will be condensed. The region around the crack is defined by a square of dimension $0.0033\times 0.0033$ centred around the crack. Model reduction of regions away from the crack is achieved by retaining every second node of every second row.

Prediction of temperature distribution in time using the adaptive condensation scheme is presented alongside prediction from the full model in Figure 18. Temperature profile along a grid of nodes parallel to the vertical axis is also presented in Figure 19. In addition to the near identical prediction from both models as can be seen from Figures 18 and 19, the temperature profile near the crack for both curves are also near identical. This contrasts with the case in sections 3.2.1 and 3.2.2 where there is a slight difference between the curves from the full and reduced models.

Figure 18.

Adaptive condensation: Temperature distribution across the plate: (a) Full model corresponding to $t=10\mathrm{s},$ (b) Full model corresponding to t = 20 s, (c) Full model corresponding to t = 30 s, (d) Reduced model corresponding to t = 10 s, (e) Reduced model corresponding to t = 20s, (f) Reduced model corresponding to t = 30s.

Figure 19.

Adaptive condensation. Temperature profile along a grid of the plate with Neumann boundary condition.

4. Conclusion

In this paper, we have been able to utilize a nonlocal vector calculus to derive the general form of the nonlocal peridynamic transient heat conduction equation. Using specific kernel functions, we show how the different explicit forms of the peridynamic heat transport equation developed in the literature can be recovered from the general form.

A model order reduction of the nonlocal peridynamic transient heat conduction model based on static condensation algorithm was implemented. The results from predictions based on the full order model and reduced order model show a near perfect match, thus demonstrating the capability of the model reduction scheme in preserving the characteristics of the system. The model reduction scheme also allows for adaptive implementation of the condensation algorithm so that a more detailed model is implemented in regions where higher resolution results are needed for greater insight into the numerical predictions.

The development in this study only treats systems with time invariant physical and geometrical characteristics. For example, model parameters such as the thermal conductivity or peridynamic bonds existing between material points must be time invariant. Because of this constraint, treatment of thermal conduction problems in the presence of propagating cracks is beyond the capabilities of the present scheme. In a future endeavour, the authors would like to extend the capabilities of this model order reduction scheme to problems with time dependent properties.