Open system peridynamics

We propose, for the first time, a thermodynamically consistent formulation for open system (continuum-kinematics-inspired) peridynamics. In contrast to closed system mechanics, in open system mechanics mass can no longer be considered a conservative property. In this contribution, we enhance the balance of mass by a (nonlocal) mass source. To elaborate a thermodynamically consistent formulation, the balances of momentum, energy and entropy need to be reconsidered as they are influenced by the additional mass source. Due to the nonlocal continuum formulation, we distinguish between local and nonlocal balance equations. We obtain the dissipation inequality via a Legendre transformation and derive the structure and constraints of the constitutive expressions based on the Coleman–Noll procedure. For the sake of demonstration, we present an example for a nonlocal mass source that can model the complex process of bone remodelling in peridynamics. In addition, we provide a numerical example to highlight the influence of nonlocality on the material density evolution.

in two-dimensional and 1/4 in three-dimensional problems [14], Silling extended the original PD theory to (ordinary and non-ordinary) state-based PD [15]. Alternatively, Javili et al. [9] recently introduced continuumkinematics-inspired peridynamics (CPD), where the interactions that govern the nonlocal material behaviour are referred to as one-, two-and three-neighbour interactions dependent on the contributing number of neighbours. Since the interactions in CPD are governed by the change of pair-length, triplet-area and tetrad-volume, the CPD framework represents a kinematically exact formulation. Since one-neighbour interactions are equivalent to bond-based interactions, the Poisson effect is captured by considering two-and three-neighbour interactions [16].
Although the range of applications for PD has grown tremendously in recent years, the number of theoretical contributions is rather limited. To name a few, besides the already mentioned theoretical contributions of Silling et. al [1,15] and Javili et al. [9], additional theoretical contributions involve the one of Silling and Lehoucq [17], where the first and second law of thermodynamics are exploited based on peridynamic quantities. Ostoja-Starzewski et al. [18] discuss the restrictions imposed by the second law of thermodynamics in bond-based and state-based PD. Further theoretical contributions, e.g. the ones of Kilic et al. [19], Bobaru and Duangpanya [20] and Oterkus et al. [21], address thermomechanical problems. Additionally, Javili et al. extended recently the CPD formulation to account for thermomechanical coupling [22] and elasto-plastic material behaviour [23].
However, to the author's knowledge, only formulations for closed system peridynamics have been introduced in recent years. Therefore, the aim of this contribution is to propose a formulation for open system (continuum-kinematics-inspired) peridynamics for the first time. In contrast to closed systems, where the system can exchange energy in terms of heat and work with its surrounding, open systems can additionally exchange matter with its surrounding and generate mass locally. Thus, closed systems represent special cases of an open system, where the conservation of mass has to be satisfied. In the literature, mainly closed systems are considered in continuum mechanics, although it is a simplification for most mechanical problems. Particularly for chemomechanical and biomechanical processes, the concept of open systems holds. Since additional mass sources influence not only the balance of mass, but also the balances of momentum, energy and entropy, all balance equations have to be reconsidered. Thus, in terms of open system mechanics it is appropriate to start with the balance of mass and continue with the balances of momentum, energy and entropy. For open system CCM, Kuhl and Steinmann [24] showed that with the help of the balance of mass all other balance equations can be reformulated to reduced balance equations, which are comparable to the ones valid in closed system CCM.
To propose a formulation for open system (continuum-kinematics-inspired) peridynamics, the paper is structured as follows. In Sect. 2, we revisit the kinematic description relevant for PD and state general forms for the global, point-wise and neighbour-wise balance equations. In Sect. 3, the governing equations valid for open system peridynamics are stated. First, the balance of mass is introduced. Subsequently, the balance of linear and angular momentum are discussed, where the balance of mass is incorporated to obtain a reduced form for both balance equations. Following this, the kinetic energy theorem is examined to identify the external and internal power that are helpful to establish the balance of internal energy. Lastly, the balance of entropy is stated that is required to obtain the dissipation inequality, discussed in Sect. 4. In Sect. 5, we obtain the structure and constraints of the constitutive expressions by applying the Coleman-Noll procedure. For the sake of demonstration, we present in Sect. 6 an example for a nonlocal mass source and constitutive expressions for modelling the process of bone remodelling in a peridynamic sense. Section 7 concludes the main findings of the paper and provides further outlook.

Kinematics in peridynamics
As illustrated in Fig. 1, we consider a continuum body in its material configuration B 0 ⊂ R 3 at time t 0 ⊂ R + with its boundary ∂B 0 . The material coordinates X that identify points of B 0 are mapped via the nonlinear deformation map y to their spatial coordinates as x = y(X, t) where B t ⊂ R 3 at t > t 0 ⊂ R + denotes the continuum body in its spatial configuration.
In contrast to local continuum mechanics, the peridynamic theory is based on the fundamental idea that every continuum point is in interaction with continuum points in its finite neighbourhood, the so-called horizon H 0 ⊂ B 0 . The interaction region is typically a spherical neighbourhood measured via its radius δ in the material configuration. The horizon H t in the spatial configuration is obtained via the deformation map by H t = y(H 0 ). The horizon size δ is considered a material parameter of PD as it influences the degree of nonlocality at X.
applies for nonlocal quantities, where ρ | 0 is the material density of the respective neighbour. Both material densities, ρ 0 and ρ | 0 , are local continuum quantities and considered homogenised versions of the underlying heterogeneous material. Note that {•} | 0 is a neighbour-wise quantity per volume squared and {•} | consequently per volume and per mass. Next, we establish general forms for balance equations in our peridynamic formulation and compare them to their counterparts in local classical continuum theory. To proceed, we first address the difference between the PD and CCM theory in terms of flux or flux-like quantities. Therefore, we consider a continuum body B 0 := B 1 ∪ B 2 exposed to external loads on its boundary ∂B 0 . To analyse the current state of a continuum point X, we fictitiously divide the body B 0 into two sub-domains, B 1 and B 2 , schematically illustrated in Fig. 2. The resultant force measures acting on the continuum point X α in CPD theory that result only from one-neighbour interactions (left) are compared to the resultant force measures acting on X in CCM theory (right). In CCM, the resultant force is measured via the Cauchy theorem by the traction vector t = P · N, as depicted in Fig. 2 (right). Therein, P denotes the Piola stress and N the outward-pointing unit normal vector on the area element dA. In PD, however, we cannot apply the Cauchy theorem due to the nonlocal characteristic of PD. As depicted in Fig. 2 (left), the force densities p | iα with i = {β, γ } originate from interactions of the continuum point X α with the continuum points X | β and X | γ , which exemplarily represent the neighbourhood of X α in the sub-domain B 2 . The resultant force is consequently computed by an integral of the force density per volume squared over the entire horizon H 0 . Note that Fig. 2 corresponds only to one-neighbour interactions in CPD or bond-based PD, wherein the interactions are pairwise by nature, unlike two-and three-neighbour interactions in CPD [9].
In continuum mechanics, the rate of a balanced quantity of a continuum body B 0 is in general equal to the sum of source terms and flux or flux-like terms [24]. Generalising the above example to "flux" densities { } | 0 as nonlocal quantities, assuming sufficient smoothness regarding the balanced quantity { * } 0 , the "flux" density { } | 0 and the source quantity {•} 0 and integrating over the continuum body B 0 , the general global balance equation in PD is given by where { } | 0 is in particular a "flux" density per volume squared that additionally needs to be integrated over the horizon H 0 .
Remark 1 In PD, the term flux has to be used very carefully, since a flux quantity is a local quantity that defines the flux through a surface. Due to the nonlocal characteristic of the PD theory the term flux density is not the physically correct description for the neighbour-wise density { } | 0 . However, we maintain the term flux and put it in quotation marks to highlight that { } | 0 is the nonlocal counterpart to the flux quantity { } per area in CCM. Consequently, we rather speak of a "flux" or flux-like quantity in PD.
After localisation, the point-wise and yet nonlocal balance equation in PD reads Note that the balanced and source quantity themselves can be local or nonlocal quantities. In case of the latter, we define these quantities also in integral form over the horizon H 0 as respectively. For balance equations with only nonlocal quantities the point-wise form reads that can be further localised to a neighbour-wise and also local form This consequently provides a more restrictive balance equation than the point-wise form (4). In contrast to Eq. (1), the general balance equation in global form in CCM is stated by integrating a flux term { } · N over the boundary ∂B 0 of the continuum body and is given by where N is the outward-pointing unit normal vector to ∂B 0 and { } is the flux quantity per area [24]. After localisation and via the Gauss theorem, the point-wise balance equation in CCM reads Remark 2 Balance equations can also be formulated in the spatial configuration. The conversion of the material and spatial configuration can be accomplished in CCM via the deformation gradient F := Grad y, its Cofactor K := CofF and its determinant J := DetF. Similar holds for CPD [9] as it provides a geometrically exact framework.

Governing equations
Based on the general forms of the balance equation, we will elaborate in the following the mechanical and thermodynamic governing equations that are valid for open system peridynamics. The respective governing equations for closed systems can be found in [9,22] and are recovered in the present formulation by setting the temporal rate of mass to zero. Table 1 provides a comparison of the governing equations in the peridynamic and classic continuum formulation.

Balance of mass
In the case of open system mechanics, the mass contained in B 0 can no longer be considered a conservative property, since in-or out-"flux" of mass or creation of mass influence the temporal rate of mass. In this contribution, we only consider a mass source R 0 for the sake of simplicity. The balance of mass in its global form is consequently given by Energy balance Dissipation inequality for the PD and CCM formulation, where the material density ρ 0 at the continuum point is the balanced quantity. For both formulations, the point-wise form of Eq. (8) reads Although the balance of mass (9) is a point-wise balance equation without any nonlocal contribution due to the omitted "flux" term, the mass source R 0 of our peridynamic formulation can be a nonlocal quantity by its constitutive expression, providing the balance of mass with a nonlocal characteristic.

Remark 3
In case of an in-or out-"flux" of mass, the balance of mass (9) needs to be enhanced by an additional "flux" term. In PD, we would introduce a mass term H 0 r | 0 dV | with the mass "flux" density r | 0 per volume squared that governs the mass "flux" in a peridynamic sense. The balance of mass is established following the general point-wise form (2) with ρ 0 being the balanced quantity and R 0 being the source quantity. In CCM, on the other hand, we would introduce a mass flux R per area as flux quantity resulting in the mass flux term DivR.

Balance of linear and angular momentum
We begin with the balance of linear and angular momentum in their global form. For both balance equations the momentum density π 0 = ρ 0 v, defined as the velocity v of the continuum point X weighted with the continuum point density ρ 0 , is or contributes to the balanced quantity.
The balance of linear momentum embodies Newton's second law and states that the temporal rate in π 0 integrated over the body B 0 is in equilibrium with the internal body force density b int 0 and the external body force density b ext 0 integrated over B 0 as Therein, the internal body force density b int 0 of a continuum point in PD is an integral over the horizon H 0 of the neighbour-wise force density p | 0 between the continuum point X and its neighbours X | as see [9]. In CPD, the neighbour-wise force density p | 0 can be further subdivided into where p 0 | 1 , p 0 | 2 and p 0 | 3 denote the neighbour-wise force densities according to one-, two-and three-neighbour interactions [9], respectively. The amount of force density invoked by the mass source R 0 is included in the external body force density by b ext 0 :=b ext 0 + vR 0 , wherebyb ext 0 is accordingly called the reduced external body force density.
Localising Eq. (10) and incorporating the expressions for the internal and external body force densities result in the point-wise balance of linear momentum In comparison, the balance of linear momentum in CCM reads Incorporating the balance of mass (9) into Eqs. (13) and (14), results in the reduced balance of linear momentum in PD and in CCM, respectively.

The balance of angular momentum in PD is derived by starting with the global form of the momentum balance
With D t y × π 0 = 0, y | = ξ | + y and the internal body force density (11) at hand, the global balance of angular momentum can be rewritten to The point-wise balance of angular momentum is obtained via localisation and reads Incorporating the balance of mass (9) into Eq. (19) yields the reduced point-wise balance of angular momentum Given the balance of linear momentum (15), the point-wise balance of angular momentum (20) can be further reduced to Javili et al. [9] discuss the required conditions for the interaction potentials of one-, two-and threeneighbour interactions such that the balance of angular momentum is satisfied a priori. With the balance of linear momentum (16) in CCM, the balance of angular momentum results in the following statement with the third-order permutation tensor ε.

Kinetic energy theorem
Next, we consider the balance equations for internal energy and entropy and establish the dissipation inequality. First, however, we examine the kinetic energy in order to identify the external and internal mechanical power. Therefore, we introduce the global kinetic energy K in the material configuration as being the kinetic energy density per volume. The mass-specific counterpart of the global kinetic energy reads where k is the kinetic energy density per mass. The temporal rate of the volume-specific kinetic energy is consequently given by with the identity Given the balance of linear momentum (10) and balance of mass (8), we obtain the kinetic energy theorem The definition of the internal body force density (11), the external body force density b ext 0 :=b ext 0 + vR 0 and the relations v · v = 2k and v = v | − D t ξ | yield the more detailed form of the kinetic energy theorem Incorporating the virtual power equivalence proposed in [9], with particularising δ y = v results in where we can identify the external mechanical power due to the mass source and externally prescribed body forces and tractions. In addition, we identify the internal mechanical power due to interaction forces as The kinetic energy theorem can eventually be summarised in terms of the externally and internally generated mechanical power as

Balance of internal energy
Let U denote the global internal energy that is obtained by the integral of the internal energy density u 0 over B 0 as and its temporal rate We consider that the internal energy density u 0 per volume of a continuum point is a nonlocal quantity, which is why it is given by an integral over H 0 of the internal energy density u | 0 per volume squared as To proceed, we briefly address the balance of total energy E = K + U . In closed and open systems the temporal rate of the total energy E = K + U is equal to the external power of the system that is composed of a mechanical and non-mechanical contribution. Denoting the external thermal power as Q ext yields the balance of total energy Note that in contrast to closed systems, we have additional external power contributions in P ext and Q ext due to the mass source R 0 .
With the balance of kinetic energy (33) and the balance of total energy (37) at hand, it follows that the temporal rate of internal energy is equal to the sum of the internal mechanical power P int and the external thermal power Q ext , that is The external thermal power Q ext summarises the thermal power induced by externally prescribed heat within the body B 0 and heat "flux" on the boundary ∂B 0 as with Therein, R ext 0 is the external heat source that is composed of the reduced external heat sourceR ext 0 and a nonlocal source term caused by the mass source R | 0 and the internal energy density u | of the neighbour as The external heat "flux" density Q ext 0 per area is equivalently expressed in the peridynamic formulation as a nonlocal quantity by with the heat "flux" density q | 0 per volume squared. In CCM, the external heat source R ext 0 :=R ext 0 + uR 0 of Eq. (40) is a local quantity, dependent on the internal energy u and mass source R 0 of a continuum point [24]. Furthermore, the external heat flux density Q ext 0 is expressed by the heat flux vector Q per area as Q ext 0 = Q · N that is comparable to t ext 0 = P · N for the mechanical problem and yields the counterpart to Eq. (42) Equipped with the detailed expressions of the external thermal power, we obtain via localisation the pointwise form in PD and its counterpart in CCM Next, we incorporate the balance of mass into both balance equations, (44) and (45). However, for the peridynamic balance equation (44) we use the balance of mass of the neighbours instead of the continuum point as the internal energy u 0 is considered a nonlocal quantity. Consequently, we obtain the reduced balance of internal energy in PD and its local counterpart in CCM

Balance of entropy
Based on the second law of thermodynamics, we introduce a new extensive quantity, the entropy of the system, which dictates the direction of a thermodynamical process. For this purpose, let S denote the global entropy that is obtained by the integral of the entropy density s 0 per volume over B 0 as and its temporal rate We consider that the entropy density s 0 of a continuum point is a nonlocal quantity, which is obtained by an integral over H 0 of the entropy density s | 0 per volume squared, that is The entropy of an open system is balanced by the input of entropy as where H ext is the external entropy contribution and H prd ≥ 0 is the positive entropy production. Similar to the external heat source Q ext , the external entropy input H ext is composed of a source term H ext B within the body B 0 and a flux term H ext ∂B on the boundary ∂B 0 as Next, we establish a relation between the external entropy source on the boundary H ext ∂B and the external heat "flux" density q | 0 in dependency on the absolute neighbour-wise temperature T | > 0, that reads The external entropy source within the body H ext B is expressed in terms of an external entropy density H ext 0 as Therein, the external entropy density H ext 0 is composed of the reduced external heat sourceR ext 0 divided by the absolute temperature T , a nonlocal source term resulting from the mass source R Introducing the dissipation density D 0 ≥ 0 of a continuum point, we can express the positive entropy production as Thus, we obtain the point-wise balance of entropy Since the entropy density s 0 is a nonlocal quantity, we further reduce the balance of entropy with the help of the balance of mass of the neighbours, resulting in Under isothermal conditions with T = T | , we eventually obtain that should be compared to its local counterpart in CCM

Dissipation inequality
We first introduce the Helmholtz energy density in the material configuration that can be expressed by the internal energy u 0 and the entropy s 0 multiplied with the absolute temperature T , i.e. via a Legendre transformation as with its time derivative Under isothermal conditions Eq. (61) reduces to since D t T = 0. Note that we only state the quantities per unit volume for the sake of brevity.
To continue, we first need to express nonlocal quantities in integral form. In PD, we consider the Helmholtz energy density 0 per volume as a nonlocal quantity that is given by for isothermal conditions with T = T 5 Coleman-Noll procedure and further reduced neighbour-wise dissipation inequality 1 ≤ n ≤ 3.5. In CCM, numerical stability and uniqueness of solutions is ensured for n < m, see [25] and [32], respectively. Accordingly, the Helmholtz energy density per volume that enters the mass source function R 0 reads For the underlying purely mechanical material behaviour, we utilise the harmonic pairwise energy density function that is typically used in PD literature [9,15], since an elastic energy density is typically employed to model open-pored hard tissue [33]. In Eq. (87), L = | | | and l = |ξ | | are the bond lengths in the material and spatial configuration, respectively. The peridynamic material parameter C indicates the resistance against the change of bond length. The more detailed expression of the mass source function in our peridynamic model for bone remodelling, clearly emphasises the nonlocal characteristic of the mass source R 0 . From the neighbour-wise dissipation inequality (81) the entropy source s | 0 follows as

Conclusion
This contribution proposes a framework for the thermodynamics of open system (continuum-kinematicsinspired) peridynamics. We introduced the balance of mass as a local balance equation that allows the system to gain or lose mass dependent on the (nonlocal) mass source. Equivalent to classical continuum open system mechanics, we assume that the change in mass affects the balances of linear momentum, internal energy and entropy. For the peridynamic formulation we distinguish between local and nonlocal balance equations. By incorporating the balance of mass into the volume-specific balance equations, we obtained the corresponding reduced balance equations, whose structures are comparable to the ones valid in closed system peridynamics. We have presented an example for a nonlocal mass source, where we exploit the constitutive equations of bone remodelling processes. In order to demonstrate the influence of nonlocality on the governing equations, a numerical example of a biaxial deformation test is shown in Fig. 3. The relative density evolution for different horizon sizes δ that results from a stepwise applied deformation function highlights the influence of nonlocality on the temporal change in relative density. Further numerical examples and information on the computational implementation can be found in [12].
In summary, this contribution proposes a thermodynamically consistent framework for open system (continuum-kinematics-inspired) peridynamics for the first time. We believe that this framework can be used to study and better understand nonlocal material behaviour in biomechanical and chemomechanical processes where the material density changes in time due to a (nonlocal) mass source.
In future research, our open system (continuum-kinematics-inspired) peridynamic framework will be exploited to model and study nonlocal material behaviour occurring in implant-bone-interfaces between nonliving and living material. Furthermore, we will use this framework to account for material density changes during bone healing processes, revisiting the initial scope of PD in fracture mechanics. Regarding the example of a nonlocal mass source, this contribution has focused, for the sake of demonstration, on one-neighbour interactions of the CPD framework of Javili et al. [9]. In future, we will also study the influence of two-and three-neighbour interactions on the density evolution. In a next step, the balance of mass can be enhanced by a mass "flux" density, allowing in-or out-"flux" of mass. The relative density evolution (bottom) for different horizon sizes δ, resulting from a stepwise biaxial deformation u =ū with a stepwise incremental strain of ε = 0.015 in horizontal and vertical direction applied on a unit square, shows the influence of nonlocality. For all simulations, we uniformly set δ/ = 3.01, ρ * 0 = 1.0, c = 1000, m = 3, n = 2, * 0 = 0.001 and E = 1.0. The peridynamic parameter C in Eq. (87) for one-neighbour interactions is computed by the relation C = 9E/δ 3 π with the Young's modulus E [34] Funding Open Access funding enabled and organized by Projekt DEAL. Ali Javili gratefully acknowledges the support provided by the Scientific and Technological Research Council of Turkey (TÜBITAK) Career Development Program, grant number 218M700.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.