Continuum-kinematics-inspired peridynamics. Mechanical problems

The main objective of this contribution is to develop a novel continuum-kinematics- inspired approach for peridynamics (PD), and to revisit PD’s thermodynamic foundations. We distinguish between three types of interactions, namely, one-neighbour interactions, two-neighbour interactions and three-neighbour interactions. While one-neighbour interactions are equivalent to the bond-based interactions of the original PD formalism, two- and three-neighbour interactions are fundamentally different to state-based interactions in that the basic elements of continuum kinematics are preserved exactly. In addition, we propose that an externally prescribed traction on the boundary of the continuum body emerges naturally and need not vanish. This is in contrast to, but does not necessarily violate, standard PD. We investigate the consequences of the angular momentum balance and provide a set of appropriate arguments for the interactions accordingly. Further- more, we elaborate on thermodynamic restrictions on the interaction energies and derive thermodynamically-consistent constitutive laws through a Coleman–Noll-like procedure. an open access article


Introduction
Peridynamics (PD) is an alternative approach to formulate continuum mechanics ( Silling,20 0 0 ) the roots of which can be traced back to the pioneering works of Piola 2016; which prepared the foundations for nonlocal continuum mechanics and peridynamics. PD has experienced prolific growth as an area of research, with a significant number of contributions in multiple disciplines. PD is a non-local continuum mechanics formulation. However, it is fundamentally different from common non-local elasticity (e.g. Eringen, 2002 ) in that the concepts of stress and strain are not present. As a non-local theory, the behaviour of each material point in PD is dictated by its interactions with other material points in its vicinity. Furthermore, in contrast to classical continuum mechanics, the governing equations of PD are integro-differential equations appropriate for problems involving discontinuities such as cracks and interfaces.
While the discretized format of PD bears a similarity to discrete mechanics formulations such as molecular dynamics (MD), it is still a continuum formulation and only takes advantage of basic MD concepts such as the cutoff radius and Fig. 1. Schematic illustration and comparison between the standard PD formulation (left) and the proposed continuum-kinematics-inspired alternative (right). One-neighbour interactions in our framework are identical to bond-based interactions in the PD formulation of Silling (20 0 0) . Two and threeneighbour interactions corresponding to Eq. (4) and Eq. (5) , respectively, are alternatives to state-based interactions. The difference between the bondbased, ordinary state-based, and non-ordinary state-based PD formulations lies in the magnitude and direction of the interaction forces (green arrows) between the materials points. In our approach, the difference between the one-, two-and three-neighbour interactions lies in their kinematic descriptions. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) point-wise interactions. For further connections and differences between PD theory, continuum mechanics and particle systems see the fundamental contributions by Fried (2010) , Murdoch (2012) , Fosdick (2013) , and Podio-Guidugli (2017) , among others. PD inherently accounts for geometrical discontinuities, hence it is readily employed in fracture mechanics and related problems. However, the applications of PD extend far beyond fracture and damage. For an extensive study of the balance laws, applications, and implementations, see Madenci and Oterkus (2014) , and for a brief description of PD together with a review of its applications and related studies in different fields to date, see Javili et al. (2018) . Table 1 categorises various PD applications and the associated key contributions in the literature. It is clear that the range of PD applications is broad and not limited to fracture mechanics.
The original PD theory of Silling (20 0 0) was restricted to bond-based interactions. This limited its applicability for material modelling, including the inability to account for Poisson's ratio other than 1/4 for isotropic materials. This shortcoming was addressed in various contributions and finally rectified by Silling et al. (2007) via the introduction of the notion of state and categorising the interactions as bond-based, ordinary state-based and non-ordinary state-based as schematically illustrated in Fig. 1 (left). Despite the large amount of research on PD, its thermodynamic foundations have not been fully investigated. Fundamental works on PD are limited in number but include those of Silling and Lehoucq (2010) , Ostoja-Starzewski et al. (2013) , and Oterkus et al. (2014a) . The starting point of these contributions is the PD theory and constitutive formulation of Silling et al. (2007) . The goal here is to adopt a continuum-kinematics-inspired approach and thereby bridge the gap between classical continuum thermodynamics and PD. More precisely, we propose an alternative PD formulation whose underlying concepts are reminiscent of classical continuum mechanics. In particular, we firstly propose to decompose the interaction potentials into three parts corresponding to one-neighbour interactions , two-neighbour interactions and three-neighbour interactions within the horizon, as illustrated in Fig. 1 (right). Note, one-neighbour interactions are identical to bond-based interactions in the PD formulation of Piola  and Silling (20 0 0) . Secondly, we derive the balance of linear and angular momentum corresponding to our interaction potentials and identify the fundamental properties of these potentials such that angular momentum balance is a priori fulfilled. Finally, we derive the dissipation inequality and propose thermodynamically-consistent constitutive laws. Crucially, we postulate the virtual power equivalence as the key requirement of our approach and build our entire framework solely on this variational assumption.
Remark. Before proceeding, we revisit the notions of a "localization procedure" and a "point-wise equation" since in the current context they serve a broader purpose than they usually do in classical continuum mechanics. Localization refers to the process of deriving a point-wise relation from an integral form over a domain. The resulting point-wise relation itself may or may not be an integral form. Applying the localization procedure on global forms in CCM renders point-wise relations at each X that are not integrals and thus are local. On the contrary, point-wise equations at each X in CPD include integrals over the horizon and are hence non-local. It is possible to apply a localization procedure on these non-local forms to derive neighbour-wise equations that are point-wise forms at each neighbouring particle's location X | . Henceforth, we use the term "local form" exclusively to indicate the point-wise quantities and equations of CCM. The term "non-local form" on the other hand refers to point-wise integral forms associated with CPD. Finally, the term "neighbour-wise form" refers to non-integral quantities and relations in CPD obtained via localization of their non-local forms.
The manuscript is organized as follows. Section 2 introduces the notation, elaborates on the kinematics of the problem and presents the geometrical aspects of the proposed framework. Here the novelty is to introduce two-and three-neighbour interactions inspired by basic elements of classical continuum kinematics. Firstly, as a motivation, we derive the governing equations using the Dirichlet principle in Section 3 via minimizing the total energy functional, for the special case of a quasistatic, conservative problem. Next, for the general case, thermodynamic balance laws are discussed in Section 4 . In particular, we detail the kinetic energy, energy and entropy balance equations. Afterwards, through a Coleman-Noll-like procedure based on the dissipation inequality, we provide thermodynamically-consistent constitutive laws. Section 5 concludes this work.

Kinematics
Consider a continuum body that occupies the material configuration B 0 ∈ R 3 at time t = 0 and that is mapped to the spatial configuration B t ∈ R 3 via the nonlinear deformation map y as in which X and x identify the points in the material and spatial configurations, respectively illustrated in Fig. 2 . Central to the PD theory, and in contrast to standard local continuum mechanics, is the non-locality assumption that any point X in the material configuration can interact with points within its finite neighbourhood H 0 ( X ) . The neighbourhood H 0 is referred to as the horizon in the material configuration. The measure of the horizon in the material configuration is denoted δ 0 and is generally the radius of a spherical neighbourhood at X . The spatial horizon H t is the image of the material horizon H 0 under the deformation map y and its measure is denoted δ t , that is H t = y (H 0 , t ) with δ 0 := meas (H 0 ) and δ t := meas (H t ) = y (δ 0 ) . Note that the horizon H 0 coincides with the point X in the limit of an infinitesimal neighbourhood and therefore recovering the kinematics of the local continuum mechanics formalism.
To be more precise and to better distinguish the PD formalism from conventional continuum mechanics, we identify the points (neighbours) within the horizon by a superscript. For instance the point X | ∈ H 0 ( X ) denotes a neighbour of point X Fig. 2. Motion of a continuum body. Illustration of classical continuum mechanics formalism (left) and the peridynamics formulation (right). The continuum body that occupies the material configuration B 0 ∈ R 3 at time t = 0 is mapped to the spatial configuration B t ∈ R 3 via the nonlinear deformation map y .
in the material configuration. The point x | within the horizon of x is the spatial counterpart of the point X | defined through the nonlinear deformation map y as (1) For our proposed framework, we identify the neighbour set of point X as These neighbours of X denoted X | , X || , X ||| are mapped onto x | , x || , x ||| , respectively. The relative positions, i.e. the finite line elements, in the material and spatial configurations are denoted as { • } and ξ { • } , respectively, where the superscript { • } identifies the neighbour, that is In addition, we define the conventional infinitesimal line elements, by a limit operation, as d X | := lim In order to overcome the bond-based restrictions of early PD formulations, and in the spirit of classical constitutive modelling, we first recall the three local kinematic measures of relative deformation, namely the deformation gradient F , its cofactor K and its determinant J , where F := Grad y and K := Cof F and J := Det F .
We now introduce three non-local PD kinematic measures of relative deformation chosen to resemble the local measures (3) .
(i) The first relative deformation measure ξ | mimics the linear map F from the infinitesimal line element d X | in the material configuration to its spatial counterpart d x | . The infinitesimal spatial line element d x | is related to its material counterpart d X | via a Taylor expansion at X as where G is the second gradient of the deformation map y . In view of our proposed PD formalism, the relative deformation measure x | − x is the main ingredient to describe one-neighbour interactions , see Fig. 3 .
(ii) Similar to finite line elements, we introduce finite area elements constructed from two finite line elements. For instance, the vectorial area element A | / || in the material configuration corresponds to the vector product of the line elements | and || as A | / || := | × || with its counterpart in the spatial configuration denoted as a  The second relative deformation measure a | / || mimics the linear map from the infinitesimal (vectorial) area element d A in the material configuration to its spatial counterpart d a | / || . An infinitesimal area element is constructed from three points within the horizon in the limit of infinitesimal horizon measure as This is essentially the Nanson's formula frequently used in conventional continuum kinematics. In our proposed framework, In a similar fashion to finite line elements and area elements, we define finite volume elements formed by three finite line elements. Let V | / || / ||| denote the finite volume element in the material configuration with its spatial counterpart being v The volume elements V | / || / ||| and v | / || / ||| are obtained by a scalar triple product, also referred to as a mixed product, of their edges as The third and last deformation measure v | / || / ||| mimics the linear map J from the infinitesimal volume element d V | / || / ||| in the material configuration to its spatial counterpart d v | / || / ||| . However unlike J that must be strictly positive, the volume elements v | / || / ||| and V | / || / ||| can be positive or negative as long as they are consistent in the sense that v | / || / ||| /V | / || / ||| > 0 must hold. The infinitesimal volume elements are formed from four points within the horizon in the limit of infinitesimal horizon measure as ] is the main ingredient to describe three-neighbour interactions , see Fig. 5 .

Dirichlet principle setting
To gain insight into the thermodynamic balance laws before investigating the general case in Section 4 , we begin with the special case of a quasi-static conservative problem. Thus, in order to set the stage and to motivate the structure of the governing equations for the important problem of a conservative system that is equipped with a total potential energy functional, we consider the Dirichlet principle. More precisely, we obtain the governing equations by minimizing the corresponding total potential energy functional via setting its first variation to zero. The total potential energy functional consists of internal and external contributions, denoted as int and ext , respectively, and is given by The internal and external contributions are detailed in Sections 3.1 and 3.2 , respectively. In Sections 3.3 the governing equations are derived and their connection to classical (local) Cauchy continuum mechanics is highlighted. The discussion on the variational setting in this section is entirely restricted to non-dissipative processes. As outlined by dell'Isola and Placidi (2011) , however, this variational setting can be extended to more generic dissipative cases using the Hamilton-Rayleigh variational principle, as will be explored in a separate contribution.

Internal potential energy
The internal potential energy of the system int is assumed without loss of generality to be separable, i.e. to be composed of the internal potential energy due to one-neighbour interactions 1 int , two-neighbour interactions 2 int and threeneighbour interactions 3 int , that is where the number in the subscript indicates the type of interaction. These contributions to the internal potential energy are now explored.

One-neighbour interactions
To proceed, we define the one-neighbour interaction energy density per volume squared in the material configuration w 1 | as a function of the relative position ξ | between two points, that is where the semi-colon delineates arguments of a function from its parametrisation. Furthermore, we define the more familiar energy density per volume as half of the integral of w 1 over the horizon H 0 , that is wherein the factor one-half is introduced to prevent double counting since we visit each point twice due to the resulting double-integration in the next step. Consequently, the internal potential energy due to one-neighbour interactions 1 int is defined by The last step holds since at any point X one-neighbour interactions with points outside the horizon vanish. Next, the variation of 1 int can be expressed as in which the previously introduced factor one-half disappears due to the variation rules on multiple integrals. Motivated by the structure of Eq. (7) , we define the force density per volume squared due to one-neighbour interactions by and therefore the variation of 1 int , using δξ | = δy | − δy from Eqs. (1) and (2) , reads We identify the internal force density per volume in the material configuration due to one-neighbour interactions b int 01 as Note, we recognize the right-hand side of Eq. (10) as an internal force density since it is the virtual power conjugated quantity to δy according to Eq. (9) . Finally, the variation of the internal potential energy due to one-neighbour interactions 1 int reads

Two-neighbour interactions
Next, we define the two-neighbour interaction energy density per volume cubed in the material configuration w 2 | / || as a function of the area element a | / || between three points, that is Furthermore, we define the more familiar energy density per volume as one third of the double integral of w 2 over the horizon H 0 , that is The factor one-third is introduced to prevent triple counting due to the resulting triple-integrals that come next. Note that the sequence of integration may be exchanged. The internal potential energy due to two-neighbour interactions denoted 2 int is defined by Again, the last step holds since at any point X two-neighbour interactions with points outside the horizon vanish. Next, the variation of 2 int can be written as in which the previously introduced factor one-third disappears and the factor one-half is introduced for convenience. The double force density per volume cubed is defined by m Importantly, 1 m is assumed to be homogeneous of degree one in a | / || so that Using the relation δa | / || = δξ | × ξ || + ξ | × δξ || from Eq. (4) , the variation of 2 int reads To proceed, we change the order of integration for the second term and relabel the quantities, which yields Motivated by the structure of Eq. (12) , we define the force density per volume squared due to two-neighbour interactions by This result should be compared with the force density per volume squared due to one-neighbour interactions (8) . The variation of 2 int with δξı = δy ı − δy reads where we identify the internal force density per volume in the material configuration due to two-neighbour interactions b int 0 2 as b int 0 2 := Again, we recognize the right-hand side of Eq. (15) as an internal force density since it is the virtual power conjugated quantity to δy according to Eq. (14) . Finally, the variation of the internal potential energy due to two-neighbour interactions 2 int reads

Three-neighbour interactions
The three-neighbour interaction energy density per volume to the fourth power in the material configuration w 3 | / || / ||| is a function of the volume element v | / || / ||| between four points and reads We define the more familiar energy density per volume as one quarter of the triple integral of w 3 over the horizon H 0 by with the factor one-fourth preventing quadruple counting due to the following quadruple interchangeable integrals. Consequently the internal potential energy due to three-neighbour interactions denoted 3 int reads Next, the variation of 3 int can be written as wherein the previously introduced factor one-fourth disappears due to the variation rules on multiple integrals and the factor one-third on the last term is introduced for convenience. The triple force density per volume to the fourth power p | / || / ||| is defined by We note that p is invariant with respect to even permutations in ξ | , ξ || and ξ ||| since We emphasize that m was assumed to be homogeneous of degree one such that the property m | / || = −m || / | holds. However, p is invariant with respect to even permutations by definition . Using the relation δv in which in the second step we changed the order of integration and relabelled the quantities. Motivated by the structure of Eq. (17) , we define the force density per volume squared due to three-neighbour interactions as This should be compared with the force density per volume squared due to one-neighbour interactions (8) and the force density per volume squared due to two-neighbour interactions (13) . The variation of 3 int with δξ | = δy | − δy reads in which we identify the internal force density per volume in the material configuration due to three-neighbour interactions b 0 3 int as b int The right-hand side of Eq. (19) is again an internal force density since it is the virtual power conjugated quantity to δy according to Eq. (18) . Finally, the variation of the internal potential energy due to three-neighbour interactions 3 int reads

External potential energy
Let ext be the external potential energy functional consisting of the contributions from both the externally prescribed forces within the bulk and tractions on the surface of the body. Note that, in contrast to Silling (20 0 0) , we allow for the externally prescribed tractions exclusively acting on the boundary of the body, i.e. they are not considered as the result of a cut-out volume within the body. We emphasize that our assumption is in contrast to, but not necessarily in violation of the standard PD since the tractions could be embedded within the internal force densities. The external potential energy ext can thus be expressed as where b ext 0 denotes the external force density per volume in the material configuration, with units N/m 3 , and t ext 0 is the external traction on the boundary in the material configuration, with units N/m 2 . Note, this format of the external potential energy is a particular sub-case of a more general case applicable to higher gradient and non-local continua as elaborated by Auffray et al. (2015) ; Javili et al. (2013a) , among others.

Governing equations
The total potential energy functional that we seek to minimize with respect to all admissible (spatial) variations δy at fixed material placement is composed of the internal and external contributions according to Eq. (6) , that is in which From the structure of the variational form (20) we can readily extract the governing equation as The internal body force density here corresponds to the stress divergence in the classical continuum mechanics formalism where b int 0 = Div P and P is the Piola stress. The variational governing Eq. (22) should be compared to its counterpart in classical continuum mechanics where or in its more familiar local form The above should, in turn, be compared to the relation subject to the boundary conditions that can be extracted form Eq. (22) 2 and will be clarified shortly. We emphasize that the requirement for the virtual power equivalence , is an underlying postulate of our framework and a key feature of this contribution. This relation allows one to introduce and prescribe external tractions on the boundary. If the external boundary is traction free or if only displacement type boundary conditions are prescribed, the right-hand side of the requirement (24) vanishes and it reduces to This is the more familiar requirement in classical PD.
Remark. Central to the state-based PD is the notion of "correspondence" which essentially states that a peridynamic constitutive model can "correspond" to a classical constitutive continuum model for homogeneous deformations. Correspondence allows one to calibrate a PD material model such that it furnishes the same result as the corresponding classical continuum model for a given homogeneous deformation. However, the constitutive correspondence framework of PD can lead to nonphysical deformation modes including material collapse, matter inter-penetration at discontinuities, and may suffer from zero-energy mode instability. These have been addressed by Tupek and Radovitzky (2014) and Silling (2017) among others. The first and most significant step in deriving a correspondence model ( Silling et al., 2007 ) is to approximate the deformation gradient F . Since the kinematics in our approach coincides with that of classical continuum mechanics, the deformation gradient need not be approximated and can be obtained exactly. Therefore, it seems reasonable that the continuumkinematics-inspired approach can alleviate some of the aforementioned issues with the correspondence framework of PD. Nonetheless, this task is beyond the scope of the current manuscript and shall be elaborated in a separate contribution.
Remark While one-neighbour interactions are the more familiar type in mechanics, multiple-neighbour interactions are commonly employed for atomistic modeling and molecular dynamics simulations. Such multiple-neighbour interactions are often described in terms of angles and bond-length instead of area and volume. We adopt a continuum-kinematics-inspired approach; our deformation measures are common continuum measures, namely length, area and volume. So equipped, various interaction energy densities can be proposed and their respective coefficients calculated via a suitable parameter identification procedure. The generic forms of interaction energy densities will be given in Section 4.2 . To aid understanding, elastic one-neighbour interactions can be viewed as the resistance against the change of length between a point and its neighbours, reminiscent of the elastic modulus in classical continuum mechanics. Elastic two-neighbour interactions can be interpreted as the resistance against the change of the area of the triangle formed by a point and a pair of neighbours, analogous to Poisson-like effects of classical continuum mechanics in two dimensions. Finally, elastic three-neighbour interactions are essentially the resistance against the change of the volume of the tetrahedron formed by each point and its triplet of neighbours, similar to Poisson-like effects of classical continuum mechanics in three dimensions.
Remark Our formulation can impose both plane-strain and plane-stress assumptions via satisfying respective boundary conditions on a three-dimensional domain. Our formulation in 2D, however, corresponds to a purely two-dimensional case similar to the surface elasticity theory of Gurtin and Murdoch (1975) , see also ( Javili et al., 2013b ). Obviously, threeneighbour interactions do not contribute in the 2D case. Note that neither "stress" nor "strain" is present in the peridynamic formulation and they can only be computed through post-processing. Therefore, the notions of "plane strain" or "plane stress" become naturally less relevant as they correspond to a local view on continuum mechanics.
Remark Our formulation is inherently non-local and, similar to classical non-local theories, can capture a frequencydependent wave speed. Thus, dispersion of elastic waves occurs naturally. Dispersion would certainly be a feature of our model and motivates further investigation, in a separate contribution, in the spirit of the analyses provided by Bazant et al.(2016 ) and Butt et al.(2017 ). The current model may inherit some of the pathological behaviours reported in this context and thus this potential issue shall be explored. An in-depth analysis is left for a future contribution.

Thermodynamic balance laws
Equipped with the virtual power equivalence (24) for a quasi-static conservative case, we can proceed to derive the thermodynamic balance laws for more general cases. Note the virtual power equivalence (24) must hold for any arbitrary δy . Among all admissible motions, we select rigid translation and rigid rotation in what follows. For a rigid translation of the body δy = δy | = const . and therefore the virtual power equivalence reduces to This can be understood as traction equivalence and serves as the boundary condition for balance Eq. (23) . Its more familiar counterpart in classical continuum mechanics, obtained using the Gauss theorem, and the traction boundary conditions are given by For a rigid rotation of the body δy = ω δ × y and δy | = ω δ × y | with ω δ = const . being the variational analogue to the angular velocity vector and the virtual power equivalence reduces to which can be understood as torque equivalence . This should be compared with its more familiar counterpart in classical continuum mechanics given by

Momentum balances
To derive the momentum balance equations of a dynamic and possibly non-conservative problem, we follow the standard procedure of classical continuum mechanics. In doing so, we begin with the global form of the force or moment balance in their integral forms and identify terms using the traction equivalence (25) and the torque equivalence (26) relations. Obviously, this process must be carried out for both the linear momentum balance and the angular momentum balance separately.
Let v denote the velocity of the material point X and ρ 0 the mass density per volume in the material configuration. The global form of the linear momentum balance reads The integral of external traction t ext 0 is now replaced by the traction equivalence (25) and the definition of the internal body force density (21) 2 is employed to obtain which yields the non-local form of the linear momentum balance via localization as Its counterpart in classical continuum mechanics is given by To derive the angular momentum balance, we start from the global form of the moment balance The integral of external traction moment y × t ext 0 is now replaced by the torque equivalence (26) to yield Using the linear momentum balance (27) , this reduces to the global form of the angular momentum balance and upon localization yields the non-local form of the angular momentum balance Its counterpart in classical continuum mechanics is given by with ε the third-order permutation tensor .

Consequences of balance of angular momentum on elastic interaction forces
Next, we explore the consequences of the angular momentum balance (28) on the interactions and the possible restrictions it imposes on interaction potentials. In particular, we investigate the conditions required for the nature of interactions such that the angular momentum balance is a priori fulfilled. The force density per volume squared p ı is additively composed of force densities per volume squared due to one-neighbour, two-neighbour and three-neighbour interactions. Thus, the expanded version of the angular momentum balance reads Each of the three integrals must vanish identically in order to sufficiently satisfy the angular momentum balance. Accordingly, we require Therefore, we next investigate separately the one-neighbour, two-neighbour and three-neighbour interactions.

One-neighbour interactions
Recall the one-neighbour interaction energy density per volume squared in the material configuration w 1 in its most generic form is a function of the relative position ξ | , i.e. the finite line element, between two points with the force density per volume squared denoted as p 1 | , that is Inserting p 1 | from Eq. (30) into Eq. (29) 1 yields the condition required to satisfy the angular momentum balance due to one-neighbour interactions with ζ 1 = ζ 1 ( ξ | ) being an arbitrary function of ξ | .
Example If the one-neighbour interaction energy density w 1 takes the form w 1 | = w 1 ( l ) with l being a square function l := | ξ | | 2 / 2 , it sufficiently fulfils the angular momentum balance condition (31) . That is with the function ζ 1 = ζ 1 ( ξ | ) defined as This is a generic example of central interaction forces corresponding to the original bond-based model of Silling (20 0 0) .

Two-neighbour interactions
The two-neighbour interaction energy density per volume cubed in the material configuration w 2 in its most generic form is a function of the finite vectorial area element a | / || among three points with the force density per volume squared denoted as p 2 | , that is Inserting p 2 | from Eq. (32) into Eq. (29) 2 yields the condition Using the identity ] for arbitrary vectors a , b and c , the condition (33) can be rewritten as The first double-integral on the right-hand side can be expressed equivalently by changing the labels as in which the last step holds since m | / || = −m | / || according to assumption (11) . The relation (35) indicates that the first double-integral on the right-hand side of Eq. (34) is equal to its negative and thus, it can only be zero. In order to guarantee that the second term on the right-hand side of Eq. (34) vanishes, we further require m with ζ 2 being an arbitrary function of a | / || holding the property ζ 2 ( a | / || ) = ζ 2 ( a || / | ) so that m | / || = −m | / || according to assumption (11) . This requirement enforces ξ | to be orthogonal to m | / || and therefore, ξ | · m | / || vanishes identically and thus the condition (34) is fulfilled a priori.
Example If the two-neighbour interaction energy density w 2 takes the form w 2 | = w 2 ( a ) , with a being a square function such as a := | a | / || | 2 / 2 , it sufficiently fulfils the angular momentum balance condition (33) . That is with the function ζ 2 = ζ 2 ( a | / || ) defined as ζ 2 := 2 ∂w 2 ∂ a which satisfies the property ζ 2 ( a | / || ) = ζ 2 ( a || / | ) so that m | / || = −m | / || and finally in which the first term can be expressed via the virtual power equivalence (24) particularized to δy = v as the power equivalence given by with P ext := Here, P int denotes the internal mechanical power due to the interaction forces and P ext the external mechanical power due to the externally prescribed forces and tractions.

Balance of energy
Let u 0 denote the internal energy density in the material configuration. The integral of u 0 over B 0 renders the global internal energy U with its rate denoted U, that is We also allow for (external) thermal power denoted by Q ext . Note that Q ext is the thermal power due to the externally prescribed heat within the continuum body and externally prescribed heat flux on its boundary. Thus, the sum of the external mechanical power P ext and the thermal power Q ext equals the sum of the rate of internal energy U and the rate of the kinetic energy K as K + U = P ext + Q ext .
Alternatively, by substituting K from Eq. (40) , the sum of the internal mechanical power P int and the thermal power Q ext equals the rate of internal energy U, that is In other words, the internal energy balance states that the internal mechanical power and (external) thermal power cause a change in the internal energy. The thermal power Q ext is composed of thermal power Q ext B within the body and thermal power Q ext ∂B on the boundary as The external thermal power term Q ext B is the integral of the heat source density R ext 0 in the material configuration. The (external) heat source density R ext 0 should be compared with the externally prescribed body force density b ext 0 for the mechanical problem. In a similar fashion, the external heat flux density Q ext 0 is reminiscent of the externally prescribed traction t ext 0 for the mechanical problem. Similar to the virtual power equivalence (24) , or its direct consequence (23) 2 for the mechanical problem, a thermal power equivalence can be established as with q | being the heat flux density per volume squared in the material configuration with units N.m/s.m 6 . The thermal power equivalence (45) should be compared with its more familiar counterpart in classical continuum mechanics, obtained via the Gauss theorem and the flux boundary condition, and given by with Q being the heat flux vector in the material configuration. Inserting the thermal power equivalence (45) and the definition (44) into the thermal power (43) results in To proceed, we substitute the thermal power Q ext from Eq. (46) and the internal (mechanical) power P int from Eq. (41) 2 into the internal energy balance (42) which yields the important global relation The non-local form of the energy balance follows via localization of the global form as This is also referred to the first law of thermodynamics. The non-local form of the energy balance (47) should be compared with its local form in classical continuum mechanics ˙ u 0 = P : ˙ F + R ext 0 − Div Q . Table 3 Governing equations of classical continuum mechanics (CCM) and continuum-kinematics-inspired peridynamics (CPD).