Anisotropic Continuum-Molecular Models: A Unified Framework Based on Pair Potentials for Elasticity, Fracture and Diffusion-Type Problems

This paper presents a unified framework for continuum-molecular modeling of anisotropic elasticity, fracture and diffusion-based problems within a generalized two-dimensional peridynamic theory. A variational procedure is proposed to derive the governing equations of the model, that postulates oriented material points interacting through pair potentials from which pairwise generalized actions are computed as energy conjugates to properly defined pairwise measures of primary field variables. While mass is considered as continuous function of volume, we define constitutive laws for long-range interactions such that the overall anisotropic behavior of the material is the result of the assigned elastic, conductive and failure micro-interaction properties. The non-central force assumption in elasticity, together with the definition of specific orientation-dependent micromoduli functions respecting material symmetries, allow to obtain a fully anisotropic non-local continuum using a purely pairwise description of deformation and constitutive properties. A general and consistent micro-macro moduli correspondence principle is also established, based on the formal analogy with the classic elastic and conductivity tensors. The main concepts presented in this work can be used for further developments of anisotropic continuum-molecular formulations to include other mechanical behaviors and coupled phenomena involving different physics.


Introduction
Continuum-molecular models in solid mechanics assume generalized distance forces as ultimate actions, and describe macroscopic behavior of bodies from a microscopic level, making use at the same time of some concepts and tools proper of mathematical field theories. The origin of this perspective lies in the classical molecular models of Nineteenth century [1,2], developed by Cauchy, Poisson and Navier [3][4][5][6] with the aim of providing an explanation per causas of elasticity [7]. They derived the equilibrium equations of elastic bodies considered as composed of molecules [8] interacting with elastic forces linearly variable with their mutual displacements. Indeed, these molecules were regarded as simple centres of forces endowed with the property of mass, their displacement being defined by a vector field in the continuum of euclidean geometry associated with the elastic body [4,8]. Within this mechanistic view, macroscopic constitutive relations derive from intermolecular properties, and the notion of stress is introduced basing it upon an hypothesis concerning intermolecular forces [8]. Two main features characterize such models: the assumption of central forces between pair of molecules, and the existence of a region of molecular activity beyond which these pairwise actions are not supposed to extend. As an inevitable consequence of the assumed structure theory [8], stress-strain relations have to respect internal constrains affecting the number of independent elastic constant, namely the Cauchy relations [8,9]. Therefore, equations of motion of anisotropic bodies contain 15 independent moduli, which reduce to a single elastic constant in case of isotropic symmetry [3,4,8,9]. In the development of the theory, it also became evident to its founders that Cauchy relations do not depend upon considering matter as continuous or discontinuous.
Poisson proposed later molecules as polar objects, capable of rotation and displacement [8]. This suggestion was worked out in detail later by Voigt [10,11] who pointed out that the drawback of the classical molecular models lies in the too restrictive hypothesis of central actions depending only on the actual distance between molecules [3,8]. According to Voigt's model for elasticity, molecules interact in pairs via a system of actions reducible to a force and a couple [3]. He obtained that the equations of a linearized elasticity theory based on molecular assumptions, do not necessarily verify Cauchy relations, the elastic constants being 21 in the most general case. Influenced by Green's energetic approach [8,12], Voigt showed also that in nondissipative processes pairwise actions can be derived from a quadratic potential which depends, in general, on the actual position and orientation of pairs of molecules. Another different model for elasticity based on corpuscular assumptions was developed by Poincaré who postulated the existence of a force function (related to a multi-body potential) depending not only on the actual position of a pair of molecules but also on the relative position of other molecules [13]. Born developed later a molecular-type model accounting for thermal and other properties of solid bodies as well as elastic properties [8]. It should be underlined that in this context the term "molecular" is strictly related to the assumption about the ultimate nature of the internal actions [14], which in this case occur at a distance. The models for elasticity here mentioned, each one with its own peculiarities, provide a description of phenomena from a perspective of a given underlying structure, making use at the same time of the powerful tools of Calculus.
After a long time, a fundamental contribution to the development of mechanical models of this kind was given by Silling [15], who proposed at the beginning of the twentyfirst century a theory of solid mechanics, namely Peridynamics, based on actions at a distance with the aim of modeling discontinuities and non-linear material response adopting a set of integro-differential equations without partial derivatives in space [15]. Such theory postulates material points in a continuum interacting in pairs within a cut-off distance, the horizon, which represents in certain sense a generalization of the concept of radius of the sphere of molecular activity characterizing molecular models of the XIX century. Indeed, since the horizon is, in general, a finite measure, the model results to be non-local, while constitutive laws for long-range interactions are defined such that the macroscopic physical behavior of the material is the result of the assigned micro-interactions properties. Even if a nonlinearized kinematics is assumed, pairwise forces can be derived from elastic central pair potentials depending on the actual distance between pairs of material points. Hence, as for classical molecular models based on central-forces, Cauchy relations do hold [15,16]. In any case, the distinctive key point of Silling's continuum-molecuar model lies in its non-local character, and in the fact that spatial derivatives do not appear in the peridynamic governing equations, with the consequence that they remain equally valid at points or surfaces of discontinuity [15]. This aspect lends itself to the description of modern problems in mechanics involving spontaneous formation of cracks [17][18][19][20][21][22][23][24][25]. It should be noted that some peridynamic fundamental concepts can also be found in other previous works on quasi-continuum models and non-local elasticity [26][27][28]. However, Silling's work introduced a general mathematical theory including elasticity, fracture, damage as well as other mechanical behaviors.
In order to overcome constitutive restrictions such as those imposed by Cauchy relations, the original peridynamic formulation (referred to as bond-based) was extended by introducing the assumption that the interaction forces between two material points depend on the displacement of all the material points within their neighborhood (e.g. the horizon regions of the two material points) [29,30]. These enriched models, named state-based and non-ordinary statebased, still consider the concept of distance action, requiring however the definition of point-wise defined deformation measures, namely the deformation states [29,[31][32][33]. This aspect represents an important breaking point with the original theory, since the simplicity given by the paradigm of purely pairwise formalism is lost [15]. Indeed, as in Poincaré theory for elasticity, due to the new assumption concerning the nature of the interactions, energy conserving state-based models are associated with multi-body potentials [34]. Other peridynamic models which require, directly or not, the definition of multi-body potentials are detailed in [35,36].
A different perspective was given instead by Gerstle [37], who proposed to maintain an analytical formulation based on pair potentials, while assuming an enriched kinematics to describe isotropic elasticity by two independent material micromoduli. Actually, as already pointed out by Voigt, material points rotation (as additional kinematic descriptor) and the non-central forces assumption in a mechanism-based molecular model for elasticity, lead to resulting governing equations which do not necessarily verify Cauchy relations [4,8]. However, the kinematic of pair interactions assumed by Gerstle et al. stems from that of an ideal Euler-Bernoulli beam model. As a consequence, the pair potential function consists of terms involving the definition of a bending micromodulus that, in general, cannot be neglected. This introduces a constitutive internal length into the model which is proper of materials homogenized as polar continua, leading to size dependent results in elasticity under general inhomogeneous deformation [38]. Diana and Casolo then proposeda polar peridynamic continuum-molecular model based on a generalized kinematic employing the definition of a noncentral pair potential consisting of three independent terms, quadratic functions of equal number of pairwise deformation measures [39]. The implication is twofold. On one hand an explicit definition of a pairwise shear deformation and shear micromodulus can be obtained, and on the other, the micromodulus associated with the mutual rotation of oriented material points' (e.g. the bending micromodulus), is independent of the other constitutive parameters and can be neglected if necessary (consistently with standard elasticity) [40]. Other polar formulations for classical elasticity have been proposed afterwards [41][42][43]. Most of them, however, still consider a kinematic of interacting material points inspired by beam models with non-vanishing bending micromoduli, these being function of the other micromoduli and/or abstract geometrical parameters.
The universal nature of integral equations used in peridynamics allows the theory to be extended also to other physical fields. This is motivated, for instance, by the need of describing transport phenomena in bodies with (possibly evolving) discontinuities, non-local diffusion in microstructured materials etc [44][45][46]. Such extensions can be obtained by adapting the integrand function expressing the ultimate action density function to the specific features of the physical phenomenon under consideration. Gerstle et al. [47] proposed a continuum-molecular (peridynamic) model for one-dimensional diffusion-type problems along with a framework which allows four coupled physical processes to be modeled simultaneously: mechanical deformation, heat transfer, electric potential distribution, and vacancy diffusion. Bobaru and Duangpanya proposed a one dimensional peridynamic heat transfer integro-differential equation using a different form of the kernel, and derived an analytical formulation for transient heat conduction by defining the peridynamic heat-flux while following a bond-based approach [44,45,48]. A model for transient advection-diffusion problems was instead presented by Zhao et al. [49]. Prakash and Seidel [50,51] proposed a coupled formulation for the study of the electrical and piezoresistive response of nanocomposite materials. A peridynamic diffusion model for damage from pitting corrosion was proposed by Chen and Bobaru [52], whereas Diana and Carvelli derived a polar electro-mechanical analytical formulation accounting for damage evolution in solids [53]. Other models of this kind can be found in [54][55][56]. All the here mentioned formulations for diffusiontype problems and coupled phenomena fall in the category of bond-based models, which means in this case that the generic pairwise action between two material points depends on the actual value of their primary field variables (e.g. temperature, electric potential, concentration etc.). As stated before, models of this kind can be associated to quadratic pair potentials. As for state-based models instead, analytical formulations applied to diffusion were presented, among others, by Oterkus et al. [57,58] and by Katiyar et al. [59].
Most of the research carried out on continuum-molecular formulations is devoted to isotropic materials, whose stiffness, conductivity and strength properties are not directiondependent, while relatively few studies on anisotropic materials can be found in literature. In the context of elasticity and fracture, mathematical formulations of peridynamic anisotropic continua, in which Cauchy relations do not apply, were proposed later by disregarding the paradigm of purely pairwise interactions [60,61], the available continuum bond-based models being characterized by the aforementioned restriction in the number of independent material constants [18,[62][63][64][65][66][67][68]. In particular, Ghajari et al. [63] proposed a peridynamic formulation based on continuous trigonometric functions for elasticity and fracture in orthotropic materials using an analytical approach similar to that detailed in [69]. Due to the central force assumption and the directional dependency laws adopted, two of the four elastic moduli which define orthotropic in-plane elasticity have to be preset, whereas square-symmetric materials [70] cannot be modeled. A mechanical continuum-molecular formulation based on pair potentials without restrictions in the number of material constants was proposed later by Diana et al. [40,71,72].
Regarding anisotropic diffusion-type models instead, theoretical studies were carried out by Seleson [46] and Mikata [73]. The latter analyzed the problem of anisotropic heat transfer focusing on the micro-conductivity function definition without making particular assumptions about the nature of the interactions. Recently, Diana and Carvelli have derived an analytical peridynamic formulation for anisotropic materials accounting for their anisotropic overall conductivity and fracture properties, together with a general procedure for the identification of the model parameters [72]. For the best of the author's knowledge, other anisotropic models based on pair potentials with similar features have not been proposed so far.
In a computational context, pair potential based continuum-molecular models lead to mathematical formulations resulting in an easier implementation and higher computational efficiency compared to models based on multi-body potentials as state-based models [74]. Actually, their use is preferred, when possible, as the large number of literature works on bond-based type models demonstrates. Moreover a purely pairwise formalism provides an intuitive mechanism-based description of macroscopic physical properties of materials, while being at the same time closely related to the mechanics of full-discrete models. This aspect may lead to the possibility to design real lattice-like systems whose physical and mechanical effective properties may be tailoring designed by controlling those assigned to the microstructure (i.e. the pairwise interactions) [75].
This paper provides a unified theoretical and computational scheme for pair potentials based continuum-molecular modelling of anisotropic elasticity, fracture and diffusiontype problems within the framework of a revised bond-based peridynamic theory with oriented material points. As important contribution of this work, governing equations of both mechanical and diffusion-type models are derived using a unified approach based on a variational formalism. An implicit meshfree implementation strategy is also detailed together with an analytical expression of the equivalent stiffness operators. Particular attention has been given to establish a consistent micro-macro moduli correspondence between material parameters of anisotropic peridynamics and classical continuum physics. The proposed method is validated against analytical solutions and experimental data, when available. By summarizing the work that we have been doing on this topic [40,71,72,76], the manucript introduces a general approach to anisotropic problems that can be used for further developments of continuum-molecular formulations to include other mechanical behaviors and coupled phenomena involving different physics. The proposed model couples the intuitive simplicity of the concept of mutual interaction in molecular models based on pair potentials (removing Cauchy relations restrictions), with the (modern) mathematical formalism of a peridynamic continuum formulation for anisotropic materials.
The paper is organized as follows: In Sect. 2 the anisotropic model for elasticity is described. The energy-based fracture model accounting for directional-dependent toughness in solids is discussed in Sect. 3, where several numerical experiments are also proposed. Finally, in Sect. 4, the analytical model for anisotropic diffusion-based problems is derived, with a focus on heat-transfer and steady-state electrical conduction.

Elasticity
A two-dimensional continuum body Ω , composed of oriented material points , is considered.
Material points and ′ within a finite distance (e.g. the horizon [15]), are assumed to interact with each other through non-central pairwise actions depending on pairwise constitutive parameters and pairwise deformation variables (Fig. 1). The set of all material points � ∈ Ω such that ‖ � − ‖ ≤ is denoted by H , which is the horizon region of . The generic vector = � − is called bond [15] or virtual fiber. A theoretical (peridynamic) model of this kind is here referred to as (polar) continuum-molecular (CM) model.
Given the reference orthonormal basis 1 , 2 , the body configuration at time t is described by the displacement field ( , t) = u i i and rotation field ( , t) defined over Ω . The applied body forces and couples are denoted by b( , t) and c( , t) , respectively. The functions ( ) = and ( ) = denote instead the densities of mass and mass moment of inertia, and are assumed constant. Three linearized pairwise deformation measures, functions of the relative position = � − , relative displacement = ( � , t) − ( , t) = � − and rotations ( , t) = and ( � , t) = � of the generic pair of interacting points, are defined. The deformation in the direction of the material points' joining line = ∕‖ ‖ is the pairwise stretch s, where n = u � n − u n is the component of the relative displacement vector along the unit vector . The pairwise angular or shear deformation is defined as the angle difference where, since t = u � t − u t is the component of along the unit vector ∶ ⊥ (Fig. 1), the shear deformation can be interpreted as the difference between the linearized rotation angle of the virtual fiber t ∕‖ ‖ and the average rotation ̄ of the interacting oriented material points and ′ . The third pairwise deformation variable depends instead on the relative rotation of two oriented material points according to the dimensional ratio which can be interpreted as the average curvature of the virtual fiber 1 .
Assuming that the material is conservative, we define the elastic pairwise potential or micropotential function such that the scalar-valued mutual actions between pairs of oriented material points can be obtained as where w( , � , t) = w( , � , t)∕‖ ‖ is the pairwise elastic potential energy function per unit distance ‖ ‖ . From Eqs. 5-7 we obtain where k n = k n ( , � ) = k n ( ) , k t = k t ( , � ) = k t ( ) and k b = k b ( , � ) = k b ( ) are the micromoduli functions or pairwise constitutive parameters depending, in the general case of anisotropic materials, on the spatial orientation = Arg ⋅ 1 + ⋅ 2 of the virtual fiber. In the above, Λ = Λ(‖ ‖) is instead the influence or attenuation function, that weights the nonlocal interactions within the spatial domain H with respect to ‖ ‖ [77,78].
At this point, it should be remarked that focusing on standard (Cauchy) solids rather than on heterogeneous materials homogenized as polar continua [79,80], the definition of an elastic pair potential term related to mutual rotations, namely w ( , � , t) , is not strictly required [39].
The general form of the nonlocal elastic energy density at , namely the macroelastic energy density [18], is obtained by integrating the pairwise potential w( , � , t) over the horizon region H of radius whereas the total macroelastic energy (e.g. the elastic potential energy of the body) is given by The Hamiltonian action integral is defined as where the density of kinetic energy is while the potential energy related to prescribed body forces and couples, in the first instance assumed conservative, is By substituting Eqs. 11, 14 and 15 in Eq. 13, considering Eqs. 1-3, and imposing the stationarity of the functional H , we obtain where δ denotes the mathematical symbol for variation [76]. Equation 16 allows deriving the field equations at ∈ Ω, whose general form reads where ( � , , � , , � , ) = ( , � , , � , ) is the force density vector function given by It is worth considering that, alternatively, could be taken directly as a pairwise deformation parameter, as reported in [39,71]. and ( � , , � , , � , ) = ( , � , , � , ) is the moment density function given by As stated previously, when considering the continuummolecular model for equivalent Cauchy-elastic materials [81], the micro-constitutive parameter associated to the deformation measure defined in Eq. 3 is not required, and is here set equal to zero (e.g. k b = 0 ). In this case, using Eqs. 8-9, the force and moment density functions can be expressed as As a consequence, both linear and angular momentuum conservation laws in Eq. 18 can be written exclusively in terms of micro-forces, with constitutive prescriptions being then needed for f n and f t only. Eq. 21 denotes the active part of the internal action, whereas micro-moments reduce to mere constraint reactions t that can be determined via pairwise balance conditions. Pairwise (mutual) or active couples b are not present if k b = 0 is assumed. A similar feature can be found in Voigt's molecular model [3,11], however here it is not a consequence of a kinematic constrain imposed to rotations within H , but it comes from the direct constitutive assumption k b = 0.

Discretized Equations: Elastodynamics
A meshfree discretization approach [18] is considered which requires the material domain Ω be divided into a set of N sub-domains Δ p , each of which associated to a particle p of coordinates x 1 p , x 2 p . Hence, a proximity search algorithm identifies particles q belonging to the p-centered horizon region H p according to the one-point quadrature scheme proposed by Hu et al. [82,83] which accounts for partial neighbor intersections 2 . The displacements p = u p 1 1 + u p 2 2 , q = u q 1 1 + u q 2 2 and rotations p , q of the material particles p and q can be collected column-wise in vector form as (20) such that, for each pair of interacting particles p and q ∈ H p , we can collect column-wise p and q as where pq is a properly defined transformation matrix [38], and the displacements and rotations collected column-wise in vector form as are aligned with the local reference basis ̂ 1 ,̂ 2 , where the unit vectors ̂ 1 ≡ and ̂ 2 ≡ . The pairwise compatibility equation relating the pairwise deformation variables collected in the vector pq = s pq pq pq ⊤ , to the interacting particles generalized displacements can be written in a compact matrix form as where ‖ ‖ pq being the distance between the generic particles p and q. The pairwise constitutive equation is instead where defines the specific elastic property of each interaction. It relates the scalar-valued mutual actions defined in Eqs. 8-10 and collected in the vector pq , to the pairwise deformations measures defined by Eqs. 1-3. The non-dimensional factor Λ pq = Λ pq (‖ ‖ pq ) controls the radial dependence of the non-local interaction between two particles, as specified in Sect. 2.
As for the balance of linear and angular momentuum of the continuum, the discrete algebraic system of governing equations in elastodynamics, and the analytical expression of the stiffness operator, can be derived from Hamilton's variational principle referred to a closed discretized domain. The variation of the first term of Eq. 13, i.e. δ[∫ can be written in discrete form as where integration-by-parts is used between the first and the second steps, and The variation of the second term of Eq. 13, i.e. δ[∫ , can be treated in discrete form as c p } T is the vector of body forces and body couples applied at particle p in the global reference system of unit vectors 1 , 2 . Finally, considering Eq. 4, the last term of Eq. 13, i.e. δ[∫ , can be written instead for a discretized body as where H denotes the number of particles q within the p-centered horizon , while pq is the partial factor for the subdomains Δ q , related to the specific quadrature rule adopted [83]. Considering Eq. 13 together with Eqs. 30, 32 and 33, the stationary condition δH = 0 gives where pq = 3 3 ⊤ is a specific topology incidence matrix for non-pairwise defined matrices and vectors [76]. The assembly operator replaces the double sum symbol, and is introduced so that each algebraic object is added to the appropriate location in properly defined global matrices and vectors [85]. In this way, Eq. 34 can be then rewritten in compact form as where the global mass matrix is given by

Anisotropic Elastic Pair Potentials and Material Micromoduli
Given the reference orthogonal basis 1 , 2 defined in Sect. 2, the classical constitutive stress-strain relation ij = C ijhk hk of a two-dimensional homogeneous linearly elastic anisotropic Cauchy continuum can be written in Voigt notation as which in component form is given by  respectively. The angle is positive if the basis rotation is anticlockwise, and the tensor ̂ in Eq. 42 is defined as According to Eq. 45, the off-axis axial Ĉ 1111 and shear Ĉ 1212 elastic moduli can be written as circular functions of the angle by where C 1111 , C 1122 , C 2222 , C 1112 , C 2212 and C 1212 are the six given elastic constants defining in-plane fully anisotropic Cauchy elasticity. The mechanistic nature of the continuummolecular model requires the definition of pairwise constitutive laws for microstructure (i.e. the virtual fibers and then the pairwise interactions) such that the overall elastic behavior of the material is the result of the assigned micro-interactions properties [76]. Besides, to assign microscopic properties on the basis of macroscopic ones, the identification of a representative cell is of the essence [76,86]. In non-local models of this kind, the representative cell is the horizon region H over the integral in Eq. 11 is defined, while the micro-macro moduli correspondence principle is based on a micro-macro energy equivalence, which allows macroscopic properties to be directly associated to microscale-defined constitutive laws [72]. Since in-plane Cauchy anisotropic linear elasticity is defined by six independent elastic moduli, it follows that k n ( ) and k t ( ) must be functions of at least six independent microelastic constants.
In order to preserve material symmetries, and in analogy with the classical continuum, we can assume 3 without loss of generality that the pairwise axial and shear stiffness k n ( ) sin 2 cos sin sin 2 cos 2 − cos sin −2 cos sin 2 cos sin cos 2 and k t ( ) , exhibit a directional dependency as Ĉ 1111 ( ) and Ĉ 1212 ( ) described by Eqs. 46 and 47, respectively where K 1111 and K 2222 are the microelastic axial moduli of virtual fibers parallel to the unit vectors 1 and 2 , respectively, whereas K 1212 is the microelastic shear modulus corresponding to the two directions defined by the aforementioned unit vectors. Differently, K 1122 , K 1112 and K 2212 are the microelastic moduli related to the axial-axial and axial-shear couplings in anisotropic Cauchy elasticity. It is worth to note that the coupling between macro shear and axial deformation is here obtained without introducing additional terms in the elastic micropotential function expressed by Eq. 4, hence without modifying the formal structure of the pairwise constitutive law characterizing the generic virtual fiber. In the following subsection, we establish analytical relations between the microelastic moduli K 1111 , K 2222 , K 1122 , K 1212 , K 1112 and K 2212 , and the six elastic constants of Cauchy anisotropic two-dimensional continua, adopting a general and consistent approach that does not require the definition of specific deformation fields.

Micro-Macro Moduli Correspondence
Let us consider a general time-independent two-dimensional homogeneous deformation field. The strain energy density of a linear elastic fully anisotropic Cauchy continuum at a generic position is then with 11 , 22 and 12 macro strain components (see Eq. 41) in the reference basis 1 , 2 .
According to Eq. 11, the corresponding quantity in the continuum-molecular model is instead given by the general integral (48) k n ( ) = K 1111 cos 4 + K 2222 sin 4 + 2K 1122 sin 2 cos 2 + 4K 1212 sin 2 cos 2 + 4K 1112 cos 3 sin + 4K 2212 cos sin 3 (49) k t ( ) = K 1111 sin 2 cos 2 + K 2222 sin 2 cos 2 − 2K 1122 sin 2 cos 2 + where, under the assumed conditions, the pairwise deformations s( ) and ( ) of the generic virtual fiber can be expressed as linear functions of the homogeneous macro strain components ij and quadratic functions of the direction cosines of the orthonormal unit vectors = n i i and = t i i (as previously defined), by which are obtained through the Cauchy-Born rule [87,88], particularized to our case 4 , and according to ̂= , where is given by Eq. 44 [72,76]. Hence, one obtains Substituting Eq. 54 in Eq. 51, and assuming k n ( ) and k t ( ) by Eqs. 48, 49, and 51 gives as general solution where N = N( ) is a scalar-valued function of the horizon , depending on the specific influence function considered. In case of dimensionless attenuation functions, N can be expressed as In what follows, we assume without loss of generality that Λ(‖ ‖) = 1 , thus Comparing Eqs. 50 and 51, and collecting the terms that multiply the same strain components ij , six independent  In the case of general in-plane orthotropy with principal material exes defined by the unit vectors ̌ 1 ≡ 1 , ̌ 2 ≡ 2 , C 1112 = C 2212 = 0 , thus K 1112 = K 2212 = 0 . Moreover, twodimensional materials of square symmetry or square-symmetric materials [70,89], called also two-dimensional cubic materials or materials with two-dimensional cubic symmetry [90][91][92] (symmetry rotations ∕2 3 ), have the further condition C 1111 = C 2222 , hence K 1111 = K 2222 . In the special case of elastic isotropy, C 1111 = C 2222 , and C 1122 = C 1111 − 2C 1212 . Therefore, Eqs. 65-70 provide two independent micromoduli depending on C 1111 and C 1212 macromoduli only, by from which it follows that, in this case, Eqs. 65-70 lead to angle invariant pairwise axial and shear stiffness. Moreover, it can be noted that while k n in Eq. 71 is always positive definite (see Fig. 3), k t given by Eq. 72 is non-negative for C 1212 ≤ C 1111 ∕3 , thus, for Poisson's ratios ≤ 1∕3 and ≤ 1∕4 in plane stress and plane strain, respectively. An important consideration is that this analytical evidence is independent of the attenuation function considered and is rather a well-known intrinsic trait of mechanistic models associated with elastic pair potentials, 5 and characterized by uniform and overall isotropic properties [14,39,[94][95][96].
Besides, this feature is proper of energy conserving fulldiscrete models such as periodic mass-spring lattices and architected materials featuring repetitive unit cells composed of rigid elements connected by elastic interfaces [14,[96][97][98][99], whose constitutive laws involve pairwise deformation measures, and postulate the existence of elastic mutual potentials allowing inter-molecular actions to be derived from centroidal kinematics. For further readings one may also refer to [100][101][102][103]. Equations 48 and 49 together with Eqs. 65-70 determine the axial and shear microelastic stiffness associated with each bond, thus with each orientation , in such a way to reproduce fully anisotropic Cauchy elasticity. Under general inhomogeneous deformations, the micropotential function w( , and C 2212 = C 1111 ∕5 . Macromoduli and micromoduli are normalized with respect to the maximum values of Ĉ 1111 ( ) and k n ( ) , respectively

Fig. 3
Normalized macro-moduli C 1j1j , j = 1, 2 of Cauchy isotropic continuum, and corresponding micromoduli K 1j1j , j = 1, 2 of the continuum-molecular model as function of the Poisson's ratio (in the studied range −1 < < 1∕2 ), for a given Young modulus E totally defines the linear elastic macroscopic behavior of an equivalent classical continuum [76], while providing a mechanism-based description of material anisotropy. In the special case of isotropic symmetry, Eqs. 71-72 hold, and Eqs. 73-74 reduce to Moreover, it is self-apparent that in the limit case of C 1111 = 3C 1212 , the conceived elastic continuum-molecular model reduces to the well-known Silling's non-local continuum with central actions [15]. It should be mentioned that arbitrary material rotations can be easily taken into account in the proposed model. Consider an anisotropic Cauchy solid, and ̌ 1 ,̌ 2 be an orthonormal basis embedded in the material such that, initially, we have ̌ 1 ≡ 1 , ̌ 2 ≡ 2 . For arbitrary rotations of the material, Eqs. 48 and 49 modify substituting the bond angle with the angle difference − , where = Arg ̌ 1 ⋅ 1 +̌ 2 ⋅ 2 is positive if denotes an anticlockwise rotations of ̌ i with respect to i . Furthermore, in case of orthotropic elasticity with principal material axes in the direction of the unit vectors ̌ i Eqs. 48 to be defined for given values of the angle such that the directional dependent constitutive parameters are assigned in a finite number of bond directions. In case of a regular grid with spacing Δx and given quadrature rule [83], the number of bond directions to be associated with k n ( ) and k t ( ) , depends on the density parameter m = ∕Δx [104] expressing the ratio between the horizon and the grid spacing itself. Higher values of the density parameter correspond to finer angular discretizations of the the continuous circular functions in Eqs. 48 and 49. Hence, m should be chosen in such a way that the effective anisotropic elastic behavior of the material is correctly described by the discretized model [76]. Given a a homogeneous macrodeformation field defined by a deformation gradient of components F ij with i, j = 1, 2 in the reference basis 1 , 2 , the macroelastic energy of the discretized continuum-molecular model should equal the classical continuum counterpart (i.e. its strain energy density) for any arbitrary rotation of the material [76]. As representative cases, we consider a uniaxial extension and a pure shear macro-deformation field. As Fig. 4 shows, to correctly represent the variation of the elastic energy density as function of the rigid rotation of the basis ̌ 1 ,̌ 2 embedded in the material, m > 3 should be adopted in discretized continuum-molecular models for anisotropic elasticity. It is apparent that higher values of the density parameter lead to higher computational costs since larger is the number of non-zero elements in the elastic stiffness operator described by Eq. 37 [105]. In what follows, unless otherwise specified, m = 5 is used (see Fig. 4).
At this point it should be noted that the analytical identification procedure (namely the micro-macro correspondence scheme) developed in this section is derived by implicitly assuming an unbounded domain, or referring to material points located in the bulk [15]. In fact, material points in Ω that have a distance D < from the nearest point on the boundary do not have a full horizon region H [15], thus their material properties result to be slightly different from those of points in the bulk. In order to take into account this important aspect in numerical simulations, a specific surface correction algorithm is required. In this work we refer for simplicity to the volume method by Lee and Bobaru [106], particularized to our case.

Crack-Tip Problem in Anisotropic Bodies
Consider the plane elastostatic problem of a horizontal semi-infinite crack under far-field loading in a homogeneous anisotropic body. The reference basis is 1 , 2 , with unit vector 1 aligned with the axis of the crack, and a local polar coordinate system (r, ) is defined at the crack tip (Fig. 5). In the Lekhnitskii formalism [107], the problems of two-dimensional anisotropic elasticity are formulated in terms of two independent analytic functions of complex variables, whose determination for given boundary conditions, require the solution of the characteristic equation in which S ijhk are the components of = −1 with defined in Eq. 40). Due to the positive-definiteness of the elastic energy, the roots of Eq. 79 are always complex or purely imaginary and occur in conjugate pairs as 1 , 1 and 2 , 2 [108,109].
The displacement field in a small region surrounding the crack tip are analytically related to the mode-I and mode-II stress intensity factors by [109,110]  When examining the near-crack front region, the universal character of the asymptotic fields at the crack-tip eliminates the need to model finite geometry specimens and particular loading conditions [40]. The dominance of the asymptotic solution allows for a boundary layer analysis involving only the near-front region, to which the effects of loading and specimen geometry are transmitted by prescribing to its boundary the displacement field from this elastic solution and the associated stress intensity factors [111,112]. This "boundary layer" analysis considerably reduces the size of the problem and allows the results to be general and applicable to finite geometry configurations characterized by their own (known) stress intensity factors [40,111]. As representative cases, we consider far-field loadings defined by mixity mode factors K = tan −1 (K ∞ II ∕K ∞ I ) corresponding to Mode I ( K = 0 ) and mode II ( K = ∕2 ). The performed analysis setup is shown schematically in Fig. 5, where a circular shape of diameter 2a, unit thickness, and crack of length a (note that this model represents a crack whose length is infinite compared to the size of the modeled region), is subjected to displacement boundary conditions (80) along its entire boundary corresponding to Eqs. 80 and 81. The spatial domain is discretized using a regular grid with spacing Δx = a∕120 , resulting in a model composed of 45344 particles. The horizon results to be < a∕25 whereas m = 5 , which lead to negligible non-local effects in the immediate vicinity of the crack-tip [105].
The computed displacement fields are reported in Fig. 6. Figure 7 shows that the displacements obtained by the anisotropic CM formulation are in excellent agreement with the analytical solution [109], demonstrating that the model is capable to describe correctly general inhomogeneous deformation fields in classical fully-anisotropic materials.
Besides the concept of stress is not required here since the deformation measures defined are work conjugates to actions at a distance and not to contact actions, the stress field 6 can be obtained from pairwise forces by particularizing to our case the original stress definition by Saint Venant and Cauchy [4,8] given in the context of classical molecular theory of elasticity. Considering for simplicity a discrete lattice of particles, let define an arbitrary plane Π of unit normal * , passing through the centroid of particle p and dividing the family region H p into two pieces denoted as H + p and H − p , respectively. The intersection between the plane Π and the sub-domain Δ p is denoted by A p . The force per unit area which one part exerts on the other through A p can be expressed as where N + and N − are the number of particles in H + p and H − p , respectively. The summation involves only the set of bonds passing through (or ending at) the cross section A p = Δxh from the positive H + p side (one may refer to [40] for further details). The normal and tangential components of the traction vector t( p , * ) defined above are the normal and shear stress given by where * denotes the unit vector perpendicular to * . As illustrative example, Eq. 83 is applied to compute the 22 stress field in the near-tip zone with assigned mixity factor K = 0 . Figure 8 shows that numerical results are consistent with the analytical solution by Sih et al. [109].
Considering Eqs. 22, and according to Voigt [11], the couple stress field is instead null, under the assumed conditions. At this point it should be mentioned that, in general, the CM model (for equivalent Cauchy materials) solution converges to the classical elasticity solution as the horizon size reduces (keeping the density factor m constant) [113]. Moreover, for a given elastic problem and in the limit of going to zero, it can be demonstrated that the rotation field ( ) is related to the macrorotation of classical elasticity such that where ijk is the alternating symbol, while 12 ( ) and 21 ( ) are the non-zero components of the (infinitesimal) rotation tensor = skew[ ( )].
Regarding the displacement map computed for the elastic problem under consideration, it is noted a fast

Elastic Behavior of Architected Meterials
The CM mechanical formulation is applied to model the effective elastic behavior of periodic mechanical metamaterials homogenized as anisotropic Cauchy continua. In particular, we consider here the two-dimensional tetrachiral blocky system by Bacigalupo and Gambarotta [76,98] as illustrative case. This architected material is made up of square blocks of side a inclined by the chirality angle (0 ≤ < ∕4) , and connected by elastic interfaces of length b and thickness ̆ , as shown in Fig. 10. Assuming for simplicity interface made of homogeneous isotropic elastic material with Young modulus E and Poisson's ratio , and plane stress-conditions, the following relations hold [76] where and are the normal and tangential overall stiffness of the interface, respectively.
In the case of symmetric macro-fields, the continualization scheme outlined in in [98] leads to the effective fourth order elasticity tensor of the square symmetry group [70] (symmetry rotations 3 ), whose components in the reference basis 1 , 2 are given by which depend on the chirality angle as well as on the constitutive parameters characterizing the interfaces. In fact, the overall anisotropic behavior of the here considered mechanical metamaterial (which is defined by the effective macromoduli C ijhk in Eqs. 87-90) arises from geometrictopological features and specific material properties of its elastic constituents [76].
The accuracy of the CM model is assessed by considering frequency analyses and comparing homogeneous CM solutions to heterogeneous finite element solutions of a microstructured geometry. The geometry studied consists of an array of 25× 25 tetrachiral periodic square blocks with a=20mm and chirality angle = tan −1 (1∕2) . In this case, since b = (1 − tan )a [98], the interface length is equal to half side of the square block. The distance between the centroids of two neighboring blocks is denoted by = a∕ cos . Isotropic material with Young modulus E = 10MPa and Poisson's ratio ≈ 1∕2 is considered for the interfaces. Uniform density c = 1000 Kg/m 3 are assumed instead for both rigid blocks and interfaces. The sample was chosen to be large enough in relation to the periodic cells to reduce overall size effects.
In the FE microstructured model, rigid-body constraints are applied to the square blocks of the tetrachiral system, while the interfaces of (vanishing) thickness ̆= a∕500 , are meshed using eight-node quadratic elastic elements, resulting in a global model composed of 380311 nodes and 115962 finite elements [76]. Regarding instead the CM computational framework, by omitting the applied external load vector in Eq. 35, and assuming simple harmonic motion =̆ exp(̆t) , the eigenvalue problem is obtained where the eigenvalues ̆ are the circular natural frequencies of the system. The lowest five eigenvectors ̆ (modal shapes) of the square lamina with homogenized elastic and inertial properties are produced by assembling the global stiffness and mass matrices and of the discretized CM model as detailed in Sect. 2.1. Hence, the modal shapes and natural frequencies of the heterogeneous microstructured square lamina obtained from finite element analysis are compared to those from homogeneous continuum-molecular simulations. The conceived anisotropic model for equivalent Cauchy materials predicts modal shapes and natural frequencies that are in excellent agreement with finite element solutions (see Figs. 11 and 12). The -convergence analysis [114] (i.e. horizon reduction keeping m fixed) shows that the CM natural frequencies seem not to be affected by the grid spacing Δx adopted (see Fig. 11). In fact, reducing the horizon size by eight times (while keeping m = 5 constant) results in minimal changes in the continuum-molecular numerical solution. However, it is noted that increasing the density parameter m (keeping = 5∕4a ≈ fixed), a slight reduction in the difference between the CM heterogeneous finite elements solution is observed. In this case, a m-convergence procedure [114] allows the meshfree numerical solution to converge to the continuum-molecular non-local solution corresponding to = [76]. An energy-based failure criterion was firstly introduced by Foster et al. [115] for isotropic materials in the theoretical framework of state-based peridynamics. Such approach was then adapted to pair potentials based continuum-molecular models with non-central interactions by Diana and Ballarini [40] for modeling fracture in orthotropic elastic materials with theoretically uniform resistance to fracture [116]. Orientation dependent fracture energy within a mechanism-based energetic failure criterion accounting for different degrees of anisotropy has been introduced instead in [72], and the resulting model successfully applied to the study of crack propagation in cortical bone tissues. In this context, the need for considering failure criteria accounting for both axial and shear pairwise deformations firstly relies on crack front which is in general locally associated to a mixed-mode deformation in anisotropic materials [40]. Therefore, a failure criterion based on critical elongation alone is not, in general, recommended for use with non-central elastic pair potentials, since it does not take into account the shear/ angular deformation of the virtual fiber. This aspect may lead to reduced accuracy in kinking angle predictions or incorrect energy dissipation during crack extension [40]. Conversely, in isotropic brittle materials at atmospheric pressure the crack front is locally associated with a mode I deformation, and critical deformation based criteria generally lead to well simulated failure conditions and realistic crack paths even in the case of overall mode II loading [36,41,[117][118][119][120]. Nevertheless, based on our experience, an energy-based failure criterion guarantees a more accurate estimation of the peak load in presence of non-central pairwise interactions [72]. A fundamental advantage of using energy-based failure criteria for materials with directiondependent properties is that they allow for decoupling anisotropic elasticity and material fracture resistance, with evident simplification of the micro parameters identification procedure [72].
The energy-based criterion adopted considers the instantaneous bond failure when the scalar-valued micropotential energy function w( , � , t) attains its critical value w c ( , ).
As previously declared (see Sect. 2), the rotational material micromodulus is not defined when modeling Cauchy-elastic materials, thus w ( , � , t) = 0 in our case. Assuming that the material resistance to fracture is orientation dependent, the critical value of the bond stored energy density is not direction-invariant, and rather depends on the virtual fiber angle . For a given orientation of the material denoted by the angle as defined in Sect. 2.2.1, the critical value w c ( ) is obtained analytically Fig. 10 Periodic square of the tetrachiral block-lattice material [98] made of rigid units connected by elastic interfaces Fig. 11 Tetrachiral block lattice with homogeneous isotropic elastic interfaces ( ≈ 1∕2 ): Natural frequencies obtained with the microstructured finite element model (black triangles) and the anisotropic continuum-molecular model for different discretization parameters by equating the fracture energy G c ( ) = G( ) to the total work required to break all the bonds per unit of fracture surface, when assumed normal to the unit vector (see Fig. 13).
Hence, we have Considering that experimental measures of the fracture energy in materials of orthotropic symmetry (symmetry rotations 3 ) is, in general, given for two orthogonal orientation of the principal material basis ̌ 1 ,̌ 2 , namely = 0 and = ∕2 (with G =0 > G = ∕2 ), w c ( , ) may be written as function the bond orientation angle by a general two parameter law defined as [72] where n = 2N , with N ∈ ℕ + , while w c1 and w c2 are the maximum and minimum values of the critical micropotential energy, namely the critical pairwise stored energy density of virtual fibers aligned with the unit vectors ̌ 1 and ̌ 2 , respectively. By substituting Eq. 93 in Eq. 92 and considering = 0 and = ∕2 , we get Adopting n = 2 , a system of two equations in the two unknown w c1 and w c2 is derived solving the integrals in Eqs. 94

Fig. 12
Tetrachiral block lattice with homogeneous isotropic elastic interfaces ( ≈ 1∕2 ): Eigenmodes obtained with the microstructured finite element model (first row) and the anisotropic continuum-molecular model adopting = 5∕4a and m = 5 (second row) [76] In the case of fracture energy independent of , G =0 = G = ∕2 = G , = 1 , and Eqs. 97-99 (same for any exponent n) reduce to the well known [40,115] with w c ( , ) = w c representing in this case a constant angle-independent function. Critical parameters w c1 and w c2 cannot be less than zero otherwise bond rupture generates energy, thus, a restriction exists on the maximum value of the ratio = G =0 ∕G = ∕2 that can be considered for a given exponent n. For instance, assuming n=2 and n=8 in Eq. 93, we get G =0 < 2G = ∕2 and G =0 < 3.65G = ∕2 , respectively, whereas n=32 leads to < 7.14 . Hence, n can be used to model materials with varying levels of fracture resistance anisotropy, as Eqn. 93 with the corresponding values of the parameters w c1 and w c2 , can handle values up to 2 and 3.65, respectively. As n is increased for a given anisotropy ratio , the transition between the two extreme values w c1 and w c2 of w c ( , ) becomes increasingly less smooth Fig. 14), and then larger values of the density parameter m may be required for an accurate angular discretization of the continuous trigonometric function in Eq. 93 [72]. Conversely, larger values of for a given exponent n result in a progressive increase of the effective anisotropy of the critical micropotential function w c ( , ) (Fig. 15). In any case, w c1 and w c2 are calculated to match, for any n, the given experimental fracture energy values of the material, namely G =0 and G = ∕2 (Figs. 16 and 17). However, since the integrals in Eqs. 94 and 95 are referred to (two) specific orientations of the principal material axes (or conversely, referred to specific crack surface orientation for = 0 ), the effective fracture energy of the material model takes intermediate values for ∕2 > > 0 which depends indirectly on the critical micropotential function adopted, as shown in Figs. 16 and 17. A similar conceptual idea can be found in phase-field models of brittle fracture, in which directiondependent crack propagation in solids can be modeled in a variational framework by considering anisotropic fracture energy functions satisfying given symmetry conditions [121,122]. However, differently from the phenomenological phase-field approach, in our case fracture is described as a mechanism-based process, since cracks nucleate and grow when a number of bond failures coalescence into a surface and propagate [115]. It should be noted that other direction dependent critical micropotential energy laws w c ( , ) can be assumed, depending on the specific behavior of the material to be modeled, and on the experimental data at disposal. When considering fracture resistance properties invariant to ∕2 3 rotations of the material reference system, Eq. 93 may be replaced by where, in this case, the extremum values w c1 and w c2 of the critical micropotential function are referred to bond directions that differ of Δ = ∕4 , whereas the fracture energies from Eq. 92 are denoted by G =0 = G = ∕2 and G = ∕4 (with G =0 < G = ∕4 ). As for Eq. 93, higher values of the exponent n correspond to higher degrees of anisotropy of the overall fracture resistance of the material which is possible to model 7 . Adopting n = 4 , the system of two equations in the two unknown w c1 and w c2 is given by which leads to 13 Schematics for computing the equivalent fracture energy of the material associated to a fracture plane orthogonal to the horizontal and corresponding to a given orientation of the material reference system where Eqs. 103 allow to model anisotropy ratios up to = G = ∕4 ∕G =0 ≈ 1.2 . A typical effective reciprocal fracture energy angular-dependency (resulting from Eq. 101), is shown in Fig. 18. As previously declared, an important aspect of the energy-based failure criterion here illustrated, is that the directional dependent elastic and failure parameters of each interaction are theoretically independent of each other. In fact, when equating the total work per unit of fracture surface and material fracture energy, the integrand function in Eq. 92 is not given in terms of microelastic moduli or micromoduli functions k n ( ) and k t ( ) , resulting in a simplification of the fracture model critical parameters calculation with respect to deformation-based failure criteria.
The status of each virtual fiber is specified by a historydependent pairwise scalar valued function μ( , � , t) [18] where w( , � , t) = w( , t) . Hence, a local damage scalar variable can be defined at each position and time t where d( , t) is the point-wise local damage variable associated to the energy-based failure criterion.
It is worth noting that circular functions in Eqs. 93 and 101 can also be adopted for modeling elastic direction dependency in orthotropic materials, with exponent n modulated to handle different degrees of material anisotropy, while conserving positive-definiteness of the axial and shear elastic micromoduli.

Numerical Experiments
In this subsection four numerical fracture experiments are described, assuming an increasing level of complexity of the constitutive and failure models adopted. Results obtained are compared with theoretical and experimental predictions, when available. In the first two benchmark problems proposed, we show the general capabilities of the CM approach in modeling spontaneous crack nucleation, crack propagation and curving in brittle materials assuming for the purposes of the discussion isotropic elastic and failure properties.
Then, the model is applied to predict crack kinking in elastic anisotropic bodies with orientation independent resistance to fracture, and subjected to different loading conditions. Finally, the proposed formulation is adopted to model fracture in anisotropic media considering both directional dependent elastic and failure properties of the material.
Numerical simulations are performed using an implicit non-linear quasi-static scheme in displacement control with adaptive step refinement [72,125]. Moreover, to avoid fracture in compressive zones a further condition of positiveness of the pairwise deformation s is considered to determine the  current status of each interaction. Hence, according to [40], the degradation function in Eq. 104 applies only to those interactions for which s is not negative. This is in a certain sense similar to the tension-compression split of the stored energy functional usually considered in phase field models of fracture [126,127].

Crack Nucleation, Kinking and Curving
The first example is a model problem introduced in Bourdin et al. [129] to highlight the ability of the variational approach to fracture [130] to recover initiation phenomena and complex crack patterns in solids. The problem has been revised in [123] and proposed by a number of authors as a paradigmatic example to illustrate the potentiality of computational frameworks based on the aforementioned theoretical scheme [127,131]. A two-dimensional square, brittle and elastic matrix with edge-length a is bonded to a rigid circular fiber inclusion of diameter a/3, as shown in Fig. 19. The fiber remains fixed, while a uniform displacement field u 2 is imposed on the upper side of the square; the remaining sides are traction free, as specified in [123]. The elastic matrix is characterized by Young modulus E = 4000MPa and Poisson's ratio = 0.2 . The assumed fracture energy of the material is instead G = 0.01N/mm, whereas a = 30 mm is adopted for computational convenience [127]. The square domain is partitioned into 57600 particles, the grid spacing adopted being Δx = a∕240 . Fracture is not allowed at constrained borders. This forbids the development of fractures exactly at the boundary of the inclusion, and may well interpret the confining effects offered in a real experimental set-up by fractional contact or gluing of the supports [131]. An incremental applied vertical displacement Δu 2 = Δu * = 10 −4 mm is considered at each pseudo-time step of the simulation, with maximum applied displacement at the end of the numerical experiment u * =ũ = 0.02mm. Results obtained show that the matrix remains purely elastic until a crack of finite length appears near the top of the inclusion (Fig. 20. The crack curves almost instantaneously and symmetrically around the circular inclusion, the onset of the cracking process being brutal because the crack appears at a theoretically non-singular point [123]. The crack then continues to propagate horizontally until complete failure of the specimen. The crack trajectory obtained is practically identical to those obtained by Bourdin et al., as shown in Figs. 19 and 20. It is interesting to note that, consistently with results reported in [123,131] the crack appear at a small distance of the order of the length-scale parameter of the fracture model (the horizon in our case and the regularization length in the computational framework of the variational approach to fracture [131]). In any case, contrary to the classical variational approach to fracture, and consistently to the real physics of brittle materials, the conceived fracture model is non symmetric under tension and compression (see Sect. 3.1).
In the second example, a notched plate with hole is examined and compared to the experimental crack path obtained by Ambati et al. [124], and with numerical results obtained by a phase-field (PF) model [128]. This is a well-known numerical experiment widely used in literature to validate computational models for mixed-mode fracture. The layout of the problem is shown in Fig. 21 where all the dimensions are expressed in terms of the length a = 5mm, with out of plane thickness h = 3a . Assuming plane stress conditions, the isotropic material parameters adopted are the following: Young Modulus E = 5980MPa, Poisson's ratio = 0.221 , and fracture energy G = 2.28N/mm. The rectangular domain is discretized using a regular grid spacing Δx = a∕10 , which leads, after removing particles from holes, to a model composed of 29345 particles.
Zero displacement boundary conditions u 1 = u 2 = 0 , are imposed to all particles in the boundary of the lower pin, whereas the vertical displacements of all particles in the boundary of upper pin are kinematically constrained to have the same vertical displacement u 2 = u * [128]. At each pseudo-time step of the simulation an initial incremental applied vertical displacement Δu 2 = Δu * = 10 −3 mm is considered, which may be progressively reduced by the adaptive-pseudo time step refinement algorithm to reach numerical convergence. The maximum applied displacement at the end of the simulation is u * =ũ = 1.2mm. As shown in Fig. 22, a main crack develops from the notch. The crack then curves and reach the large off-centered hole. Later, a secondary straight crack appears from the hole to the vertical right edge. Results obtained agree very well with experimental data reported in Fig. 21, where the envelope of the experimental crack paths found for the four tested samples by Ambati et al. [124] is marked by the gray area. Moreover, it is noted from Fig. 23 that load-displacements curves corresponding to different CM density parameter m = ∕Δx are consistent with results by Kakouris & Triantafyllou [128] obtained using a phase-field model with the same displacement step increment, boundary conditions and material properties.
In the third set of numerical experiments the problem of crack kinking in anisotropic elastic bodies is examined. The layout of the problem is the same as in Fig. 5, and used for validating the anisotropic elastic model in Sect. 2.3.1. By increasing the magnitude of the prescribed far-field displacements given by Eqs. 80 and 81 so that the failure/fracture criterion is satisfied, the "pre-existing" crack of Fig. 5) extends in a direction dictated by the angle K (Fig. 24). Given a mode mixity ratios K = tan −1 (K ∞ II ∕K ∞ I ) , the crack extension direction depends on the elastic anisotropy of the material and may be influenced also by potential anisotropy and spatially random distribution of its surface energy [40]. For the purposes of this discussion, and in order to guarantee consistency among the results compared and corresponding to different materials, illustrative examples assume always anisotropic elasticity, with homogeneous and (directionally) uniform resistance to fracture. This is similar to the approach used in [108,116]. The simulations are performed considering a squaresymmetric material and a fully anisotropic material subjected to several mode mixity ratios K ranging from 0 • to 90 • . The orthotropic material of square-symmetry is a silicon crystal of orientation ⟨100⟩ , whose two-dimensional plane stress elastic tensor reads [132] and whose resistance to fracture is assumed uniform within the surface plane considered [132].
The fully anisotropic material instead is characterized by same ratio among elasticities as that considered in Sect. 2.3.1 (assuming C = 100 GPa). Specifically, the elastic tensor in this case is The magnitude of the K-field displacements applied at the circumferential boundary of the near tip-region is progressively increased until a small extension of the main crack appears. Such crack extension should be on one hand small enough not to alter significantly the K-field overall conditions transmitted by the boundary layer during crack propagation, and on the other should be sufficiently large to allow for its direction measurement. Taking into account these considerations, a square of edge a/2 and centered at the crack-tip is taken as a reference for kinking angle K measurement. The computed kinking angles are also compared with theoretical predictions by maximum-hoop stress intensity factor criterion (HSIF-criterion) [108,133], according to which, the crack will extend along a path whose tangent forms an angle perpendicular to the maximum circumferential (hoop) stress. Actually, this fracture criterion relies on the definition of two stress intensity factors; the so-called "hoop" K (HSIF) and "shear" K r (SSIF) stress intensity factors, respectively (cfr. Sect. 2.3.1 and Fig. 5) [108,134]. The first involves the circumferential stress and the latter the shear stress, HSIF and SSIF being (for anisotropic materials) more convenient quantities than the commonly used Modes I and II stress intensity factors, since they uncouple the Modes I and II on planes at suitable angles relative to the main crack [108,135]. Kinking angles corresponding to HSIF-criterion can be determined analytically as function of the apparent stress intensity factors ( K ∞ I , K ∞ II ), and material parameters, as specified in [40,108]. Results obtained are reported in Figs. 25, 26 and 27. It is noted that CM crack extension directions are consistent with HSIF K angle predictions. Another consideration is that, for the specific materials

Anisotropic Fracture in Cortical Bone
In this section, the CM formulation is adopted to model fracture in cortical bone tissues. Bone is a complex hierarchical composite which has natural mechanisms to resist fracture [137]. Actually, the simultaneous activation of toughening mechanisms by collagen fibres, lamellar structure of collagen fibres and osteon (Haversian canals) at various length scales provides enduring strength and toughness (Fig. 28).As a result, the material exhibits specific anisotropic properties which limit the possibility of formation of fracture planes normal to the anatomical axis of bone which defines the so called longitudinal direction [137]. Numerical analyses are performed to simulate experimental crack paths obtained by Behiri and Bonfield [138,139], which highlight the orientation dependence of the fracture properties of bovine cortical bone tissues. The authors used the compact tension (CT) test to produce stable fracture propagation in fluid-saturated specimens at various orientations of the material, whose principal directions are those parallel and perpendicular to the osteons (Fig. 28). Since cortical bone is, in general, stronger in the longitudinal direction, the load conditions used in CT tests result in final crack paths almost parallel to the osteons [63,72,139].
Orthotropic elastic moduli of bovine bone were obtained by Van-Buskirk et al. adopting a pulse transmission ultrasound method [140]. Assuming the principal orthonormal basis of unit vectors ̌ 1 and ̌ 2 directed along the longitudinal ̌l , and transverse radial ř directions (Fig. 28), respectively, the two-dimensional elastic tensor in compact form is The orientation of the principal material reference system with respect to the reference basis is denoted, as usual, by the anisotropy angle (Fig. 28). Regarding fracture properties, energy-release rates G =0 and G = ∕2 as defined in  Sect. 3 are calculated, consistently with the elastic moduli adopted, by [72,116] (109) where S iiii , i = 1, 2 are the conventional compliance moduli and 1 , 1 or 2 , 2 are the conjugate pairs root of Eqn. 79 with S 1112 =S 2212 = 0 . According to Fig. 13, the critical stress intensity factors K cľ and K cř in Eq. 109 refer to longitudinally ( = 0 ) and transversely ( = ∕2 ) oriented bovine cortical tissues, respectively. The values K cľ = 6.5MPa ⋅ m 1∕2 and K cř = 4.0MPa ⋅ m 1∕2 adopted in this study are those measured experimentally by Behiri and Bonfield using grooved CT specimens of bovine bone [138,139]. Same values can be also found in [63]. In fact, elastic and failure properties of this specific material, homogenized as a Cauchy orthotropic material, have been considered in other previous studies for validation of anisotropic models for fracture based on pairpotentials [63,64]. However, differently from [63,64], the model here considered does not suffer of restrictions affecting the number of independent micromoduli which define the effective in-plane elastic orthotropy of the material [72]. This allows to set all the four elastic macro-moduli given by the experimental study by Van-Buskirk et al. [140] without introducing approximations imposed by the presence of internal constrains. Moreover, the anisotropic fracture model adopted in this work is based on the definition of a energybased failure criterion, with the advantage of decoupling elasticity from failure properties of the material in the analytical identification of the model parameters. Compact tension tests on cortical bone specimens at different orientations , as performed by Behiri and Bonfield [138], have been simulated in displacement control and quasi-static conditions [125] using the conceived anisotropic CM formulation. The microelastic moduli are calculated using Eqs. 65-70, the micromoduli functions being defined by Eqs. 77 and  78. Considering the peculiarities of this material, critical micropotential energy function was determined using Es. 93 along with Eq. 99 [72]. The results, obtained adopting an almost uniform grid spacing Δx = 0.3 mm and m = 5 resulting in a model of 12312 particles, are shown in Figs. 29, 30, and 31 along with experimental final crack paths corresponding to angles of anisotropy = 0 , = ∕4 , = ∕3 and = ∕2 . The crack paths predicted using the CM model are in very good agreement with experimental data [138]. As can be seen from Fig. 29, experiments show quite complex fracture patterns especially in case of = ∕4 and = ∕2 . The continuum-molecular model is able to describe the progressive reduction of the average inclination angle of the fracture surface in cortical bone specimens inclined at ∕4 , as well as the initial deflection that the crack profile undergoes before propagating straight along the strongest direction in case of = ∕2 . It is noted that even with coarser meshes, crack profiles correlate well with experimental data [72]. In case of = 0 , experimental and numerical data show an almost collinear extension of the pre-existing main crack of the CT specimen. The theoretical peak load P can be then calculated by [63,141] where h is the thickness, ȃ is a geometrical parameter depicted in Fig. 28, whereas Y = 6.554 is a non-dimensional constant depending on the geometry of the specimen [141]. The peak load P = 98N is determined from Eq. 110 for K cř = 4.0MPa ⋅ m 1∕2 and thickness h = 1 mm. This theoretical prediction, calculated analytically on the basis of (110) P = K cř hȃ 1∕2 Y explicitly given experimental data [141], is in very good agreement with the peak load value P CM = 96 N obtained by the CM model (Fig. 32).

Diffusion-Based Problems: Heat Transfer and Electrical Conduction
Let us consider the two-dimensional continuum solid body Ω described in Sect. 2, focusing here on its thermal conductive properties. Material points and ′ within a finite distance , are postulated to interact with each other through pairwise exchanges of heat energy, function of pairwise conductivity properties and pairwise thermal parameters. The homogeneous body configuration is now described by the temperature field T( , t) defined over Ω . The rate of introduction of heat energy per unit volume is denoted by T ( , t) . The constant parameter c T denotes instead the specific heat capacity of the material. The thermal diffusion model is based on the definition of a pairwise scalar thermal parameter defined by that can be considered as the equivalent of a pairwise deformation measure in elasticity, and is function of the bond temperature gradient g( , � , t) = g( , � , t) , with as defined in Sect. 2. We define a specific pair potential function for heat transfer such that the thermal response function or pairwise thermal action (having dimensions of a heat current per unit volume squared) can be derived as where T = (T � − T) , while Λ = Λ(‖ ‖) is the attenuation function defined in Sect. 2. In Eqs. 112 and 113, k T = k T ( , � ) is the pairwise conductivity (assumed not to depend on temperature for simplicity), namely microconductivity of the heat transfer model, having the dimensions of a thermal conductivity per unit volume. As for the microelastic functions in elasticity, the microconductivity function in the general case of non-isotropic materials depends on the spatial orientation = Arg ⋅ 1 + ⋅ 2 of the pairwise interaction considered, such that k T ( , � ) = k T ( ).
The pairwise micro-heat flux density (vector) can be then defined according to Fourier's law, particularized to our case, as The thermal macro-potential density function is obtained by integrating the pair potential w T over the horizon region H of radius The integral I can be defined [142][143][144]

Fig. 29
Compact tension test on bovine cortical bone tissues: Experimental crack paths [138] (first row and gray-filled triangles) and normalized vertical displacement map of cracked specimens obtained by CM anisotropic model (second row) corresponding to different anisotropy angles [72] If we let δI = 0 in Eq. 116, while assuming Ṫ as invariant [142,145,146], after some manipulations, the field equation for time-dependent heat transfer at can be obtained 8 where f T (T, T � , , � ) is given by Eq. 113. Equation 117 represents the conservation of energy law for a continuummolecular model based on pair potentials. Steady-state conditions are obtained setting the right term of Eq. 117 equal to zero. For the steady-state electrical conduction problem we can consider the same formal structure of the heat-transfer problem, replacing Eq. 111 by where denotes the electric potential while t represents here the pseudo-time step in equivalent quasistatic conditions. The pairwise bond electric field is then e( , � , t) = e( , � , t) . The electrical micropotential is instead defined as where k e is the microconductivity having the dimension of an electric conductivity per unit volume. The electrical pairwise interaction between pairs or material points has the dimensions of an electric current per unit volume squared, and can be derived in a manner analogous to that of the pairwise thermal interaction. The first derivative of the micropotential function in Eq. 119 gives The balance equation at is where e is the net current per unit volume at . The current density and the heat flux can be derived as non-primary quantities using the approach described in [44,53].

Discretized Equations: Transient and Steady-State Problems
Using a meshfree discretization procedure, the domain Ω is divided, as usual, into N sub-domains Δ p associated to particles p. Particles q belonging to the p-centered horizon region H p are then selected as described in the previous sections. Considering the heat-transfer problem, the primary scalar field variables of the material particles p and q are denoted by z p = T p and z q = T q , respectively. If the particle q ∈ H p , the primary field variables can be collected column-wise in vector form as The pairwise compatibility equation relating the pairwise thermal parameter to the interacting particles temperatures can be written in compact form as where The scalar pairwise constitutive equation of the discretized system is instead where k pq T is the microconductivity characterizing the specific interaction between particles p and q. The algebraic system of governing equations of the transient heat transfer problem and the equivalent stiffness operator can be derived invoking δI = 0 subject to the condition specified above (see [142,146]), and referred to a discretized system. Hence, for the third term of Eq. 116, we obtain where the variation is taken keeping the quantity Ṫ p fixed [142,145,148]. The variation of the second term of Eq. 116 can be treated in discrete form as where Tp represents the heat source term at particle p. Finally, for the macro-potential term in Eq. 116 ⊤ is a specific topology incidence vector for point-wise defined vectors and scalars. As detailed in Sect. 2, the subscripts p and q are associated with different degrees of freedom of the system, thus the assembling symbol in Eq. 129 denotes that the algebraic objects have to be assembled properly through superposition. After global assembling, Eq. 129 can be rewritten in compact form as where T is given by whereas the equivalent stiffness operator corresponding to the whole body is given by and the global heat source vector is Moreover, the global primary field vector is given by (127) In steady-state problems, Eq. 130 reduces to T = T . The steady-state discretized equations can be applied directly to steady-state electrical conduction, provided that the primary field in this case is represented by the electric potential. Eq 122 is then replaced by Moreover, due to the convention used for the direction of the electric field, Eq. 124 should be multiplied by minus one.

Micro-Macro Conductivities
In the classical theory of homogeneous anisotropic materials, the thermal (electrical) conductivity is direction-dependent and the heat-flux (electric current density vector) is not necessarily parallel to the temperature gradient (electric field).
Given the usual reference orthonormal basis 1 , 2 , the classical constitutive equation (Fourier's law) relating the heat flux q and the thermal potential (temperature) gradient T( , t) = reads hence as Considering a generic basis ̂ 1 ,̂ 2 rotated by an angle with respect to the reference basis, Eq. 137 can be rewritten as where q = ⊤ T q and ̂= ⊤ T , with The thermal conductivity tensor ̂ T in Eq. 138 is defined as According to Eqs. 138-140, the off-axis conductivity K T 11 can be written as function of , in terms of the three conductivities K T 11 , K T 12 and K T 22 by Z (134) 141) K T 11 ( ) =K T 11 cos 2 + K T 22 sin 2 + 2K T 12 sin cos As in Sect. 2.2, to preserve material symmetries, we can assume here that the microconductivity characterizing each interaction of the proposed continuum-molecular model exhibits a directional dependency as K T 11 ( ) described by Eqs. 141 where K T 11 and K T 22 are the micro-conductivities of bonds aligned with the unit vectors 1 and 2 , respectively, whereas, K T 12 is the microcondictivity related to the non-symmetric conductive properties in anisotropic materials. Thus, three independent microcothermal moduli are defined for reproducing the overall conductive behavior of a two-dimensional anisotropic solid. The relation between the microthermal moduli, K T 11 , K T 22 , K T 12 (For electrical conduction, the three micro-conductivity parameters or microelectrical moduli are denoted by K e 11 , K e 22 , K e 12 ) and the classical macro-conductivities of anisotropic continua is obtained following a general approach consistent with that detailed in Sect. 4.2.
Consider a general time-independent two-dimensional homogeneous thermal gradient of components 1 , 2 . Under these conditions, g( ) =̂1( ) , thus the pairwise thermal measure g( ) of a fiber inclined at , is related to the homogeneous potential-gradient components i according to ̂= ⊤ T . Hence, we have where, consistently with Eq. 53, g( ) = n ( ) . Substituting Eq. 143 in Eq. 115, and assuming k T ( ) by Eqs. 142, Eq. 115 gives whose general solution is where h denotes, as usual, the out of plane thickness and N T is a scalar-valued function of the horizon , related to the specific influence function Λ considered. In case of dimensionless influence functions, N T can be expressed as (142) k T ( ) =K T 11 cos 2 + K T 11 sin 2 + 2K T 12 sin cos (143) g( ) = 1 cos + 2 sin (144) where Ξ T ∈ ℝ + . In the following, as for the mechanical model, we assume without loss of generality that Λ(‖ ‖) = 1 , thus The thermal macro-potential function W T has the dimensions of a rate of entransy per unit volume, which is related to the specific heat transfer ability of the material [149]. The corresponding quantity in classical continuum physics can be computed by Comparing Eqs. 148 and 145, and collecting the terms multiplying the same temperature gradient components i , three independent equations expressing the macro-conductivities in terms of micro-conductivities are obtained where the constant C T is By solving the algebraic system given by Eqs. 149-151 we obtain In the special case of isotropy, we have K T 12 = 0 and K T 11 = K T 22 = K T . Hence K T 12 = 0 and K T 11 = K T 22 = K T , where which is consistent with results for isotropic materials reported in [44,53] and obtained using different approaches and procedures. Given the principal material reference defined by the orthonormal basis ̌ 1 ,̌ 2 , inclined at with respect to the reference basis, Eq. 142 takes the equivalent form wherě K T 11 and Ǩ T 22 being the eigenvalues of the thermal conductivity tensor T , that according to the Jacobi rotation method, are given by The angle of anisotropy denoting the orientation of the principal material axes is instead Adopting Eqs. 157-162, the effective material behavior is defined by the two microconductivities K T ii , i = 1, 2 , and the anisotropy angle . As for elasticity, the continuummolecular model provides a mechanism-based description of direction-dependent conductive properties of materials, as the overall anisotropy is the result of pairwise properties assigned at a lower scale.
In the eventuality of negative values of the micro-conductivity moduli resulting from Eqs. 158-159, the equivalent stiffness operator may result non-positive definite. This may introduce material instability, especially in presence of non-linearities. At this point it should be reminded that Eq. 157 results from the assumption of trigonometric dependency of k T ( ) on (for a given anisotropy angle) as that of K T 11 ( ) . However, other different circular functions describing the direction-dependent microconductivity can be assumed, provided that the material symmetries are respected. For instance, a more general description of k T ( , ) as function of the bond spatial orientation can (157) k T ( , ) = K T 11 cos 2 ( − ) + K T 22 sin 2 ( − ) be obtained adopting a trigonometric law similar to that of Eq. 93 where n = 2N , with N ∈ ℕ + (Eq. 157 is got assuming n = 2 in Eq. 163). Larger values of the exponent n are associated with higher order anisotropic microconductivity functions k T ( , ) , whose moduli K T 11 and K T 22 have to be calculated, as usual, to reproduce the effective conductive properties of a material. This reflects itself in the possibility to model materials with different levels of effective anisotropy, as Eq. 163 leads to non-negative micro-conductivities for Ǩ T 11 ∕Ǩ T 22 ≤ n + 1 . As an example, assuming n=4 we obtain whereas n = 12 gives Therefore, adopting Eqs. 164-165 and 166-167 together with Eq. 163, the model is able to handle overall anisotropy ratios Ǩ T 11 ∕Ǩ T 22 up to 5 and 13, respectively, while avoiding negative parameters [72]. Different exponents n with their corresponding positive-definite micro-conductivities lead, in principle, to same results since k e1 and k e2 are, in any case, derived analytically to match the effective anisotropy of the material. However, given an anisotropy ratio Ǩ T 11 ∕Ǩ T 22 , incremental values of n result in microconductivity functions k T ( , ) characterized by a progressively faster transition between its extremum values. For this reason, when considering discrete lattices of particles instead of material points in a continuum, the exponent n should not be too large with respect to Ǩ T 11 ∕Ǩ T 22 , in such a way to ensure an accurate angular discretization of the microconductivity function for a given density parameter m = ∕Δx and arbitrary grid orientations (see Sect. 2.2.1). Another consideration is that when n = 2 in Eq. 163, the function k T ( , ) takes the form of that proposed in [18,63] for elasticity models based on pairwise central-forces. Moreover, the directional-dependent function obtained using eight-order spherical harmonic expansions proposed in [63] for two parameters orthotropic elasticity, can be obtained from the proposed two parameters anisotropic function by setting n = 12 in Eq. 163.
As stated before, the equations derived in this subsection apply directly to electrical conduction and can be extended to other diffusion-based problems. Consistently with previous sections, in what follows the density parameter is m = 5 , whereas the volume method is adopted as surface correction algorithm.

Benchmark Problems
In order to demonstrate the effectiveness and accuracy of the CM model for diffusion-based problems in anisotropic materials, two illustrative numerical examples are considered in this section. The first example consists of a Dirichlet problem on a circular anisotropic two-dimensional domain, whereas the second set of numerical experiments deals with damage sensing and crack-length determination via potential drop technique in bone tissues.

Heat Transfer in an Anisotropic Unit Disc
Given the usual reference basis 1 , 2 , let us consider an anisotropic medium with thermal conductivity tensor defined by [150,151] According to Eq. 162, the principal material reference system defined by ̌ 1 ,̌ 2 is inclined at = ∕8 with respect to the reference basis (Fig. 33). The Dirichlet problem is solved in a circular plane domain Ω = (x 1 , x 2 ) ∶ x 1 2 + x 2 2 < 1 of unit radius a. The following temperature distribution can be shown to satisfy the corresponding governing heat conduction equation [150] which leads to inhomogeneous temperature gradient field defined over the spatial circular domain Ω . Adopting a regular grid spacing Δx = a∕100 , the circular domain is discretized into 31428 particles. Micro-macro conductivities correspondence is established by Eq. 163, where n = 34 is adopted, the effective anisotropy ratio being Ǩ T 11 ∕Ǩ T 22 ≈ 33.96 . Eq. 169 is used to impose the boundary conditions to a circumferential boundary layer of thickness = mΔx = a∕20 . Results provided for temperature T distribution are reported in Figs. 34 and 35, where numerical solution is compared with exact analytical solution from Eq. 169. As for stress in elasticity, pairwise constitutive equations in continuum-molecular formulations for heat transfer do not involve a point-wise defined heat flux. However, considering for simplicity a discrete lattice of particles, it can be computed as postprocessing quantity by where H T+ denotes the number of particles q within the p-particle's horizon that have higher temperature with respect to particle p [48,53].
In classical continuum physics, the heat flux vector is instead given by Eq. 136, which considering Eqs. 168-169 gives where The heat flux magnitude ‖q‖ computed by CM model using Eq. 170 within the circular domain under consideration is reported in Fig. 36. Exact solution obtained from Eqs. 171 and 172 is compared with numerical predictions in Fig. 37, which demonstrates the consistency between Eqs. 136 and 170.

Damage Sensing in Conductive Tissues
The internal structure of bone affects not only its mechanical properties but also its effective conductive behavior [72,152]. Reddy and Saha [153] evaluated the electrical conductivity of compact (cortical) bone in three orthogonal directions parallel to the principal directions (e.g. longitudinal, radial, and circumferential) of a long bone, indicating the overall anisotropic character of this material. The electrical properties were found to be substantially dependent on bone moisture content and signal frequency applied to the test specimen, but essentially frequency independent below 10 kHz [153]. According to Reddy and Saha, in fluid-saturated (wet) compact bovine bone at lowfrequency or DC conditions, the effective resistivity in the radial direction is roughly three times that in the longitudinal direction [153], which is consistent with Chakkalakal et al. results [152]. In particular, as reported in [153], electrical resistivity in the radial direction is ř =540 Ω m, whereas in longitudinal direction ̌l =170 Ω m, with ř ∕̌l ≈ 3.2 . Such electrical resistivity anisotropy is expected, compact bone material being characterized by higher porosity in the longitudinal direction due to the presence of Haversian canals [153].
Considering the principal orthonormal basis of unit vectors ̌ 1 and ̌ 2 directed along the longitudinal ľ and transverse radial ř directions (Fig. 38), the two-dimensional conductivity tensor is then given by The angle of anisotropy is instead denoted by (Fig. 38). We revisit here the numerical experiment in Sect. 3.1.2 by coupling the CM mechanical and electrical conduction models in the hypothesis of small deformations, [53] to study the effect of evolving discontinuities in cortical bone tissues effective conductive behavior. Potential drop crack length determination technique [154,155] is used to correlate electrical output data to mechanical damage parameters.
The electromechanical setup under consideration is detailed in Fig. 38, which shows the arrangement of direct electrical current input leads and potential probe output leads. Dimensions of the specimens, mechanical boundary conditions, elastic and fracture parameters of the material as well as discretization parameters of the model are the same considered in Sect. 3.1.2.
The electrical microconductivity function k e ( , ) is determined using Eqs. 163, 164 and 165. The set of discretized equations of the quasi-static coupled electromechanical problem under consideration are obtained by superposition of Eqs. 35 and 130 neglecting inertial and time dependent terms. In particular, the coupled stiffness operator characterizing each pairwise interaction is given by where, according to Eqs. 37  while ̆ m and ̆ e are specific topology incidence matrices for coupled problems [72]. As can be seen from Eq. 174, the pairwise scalar valued function μ defined in Sect. 3 affects As can be observed from Fig. 40, the average slope of the Δ ∕Δ 0 curves depend on the angle of anisotropy , on the crack angle and crack pseudo-velocity. The latter effect can be excluded by reporting Δ ∕Δ 0 as function of the crack length b k (see Fig. 38) measured during the simulation. In this way, by detecting the average slopes of the Δ ∕Δ 0 -b k ∕a relationships (Fig. 41), it would be possible to get a direct determination of the actual crack length b k from the value of Δ registered by potential probe output leads.

Coupled Phenomena and Further Developments
For thermo-mechanical problems, the scalar-valued pairwise force function f n ( , � , t) defined in Eq. 8 can be modified including the temperature effect as where T = T ( , � ) is the linear thermal expansion coefficient along the direction of the bond between and ′ , whereas T is the average temperature change of the virtual fiber itself [58,156]. Field equations are then modified accordingly to include the thermoelastic constitutive relation and the contribution from deformational heating and cooling, as detailed by Oterkus et. al [58].
(176) f n ( , � , t) = Λk n [s( , � , t) − TT ( , � , t)]  A piezoresistive behavior could be modeled instead considering the electrical microconductivity k e ( , ) as function of an elastic deformation measure [50,53]. In this case, the actual length of the fiber may be required in the definition of the pairwise electric field e as well as in Eq. 118 [50,53].
Regarding elasticity, the definition of the third pairwise deformation parameter in Eq. 3, allows the model to be extended to microstructured materials homogenized as polar continua, for which a micro-macro correspondence involving bending moduli is also rerquired. In that case, it would be which leads to elastic size effect proper of two-dimensional discrete heterogeneous materials, the resulting constitutive internal length being, in principle, independent of the size of the horizon. In fact, these non-local effects are not due to the integral nature of balance Eqs. 17 and 18, and then do not vanish when the horizon goes to zero. This fundamental feature characterizing the conceived continuum-molecular framework would also give the possibility to model Cosserat anisotropic materials and to include chiral effects even associated to overall isotropic elastic properties. Besides, since piezoelectricity may be related to chiral asymmetry [157], this aspect would permit to estabilish a mechanism-based description of such complex phenomenon.
Consider here for simplicity the theoretical case of an orthotropic material within the context of couple stresselasticity [80]. Assuming unit vectors of the reference basis aligned with the principal material axes, the constitutive

Conclusions
In this paper we have proposed a unified scheme based on pair-potentials for continuum-molecular modeling of anisotropic elasticity, fracture and diffusion-type problems within the framework of a revised peridynamic theory with

Fig. 37
Steady-state heat transfer in an anisotropic unit disc: Detail of the heat flux magnitude ‖q‖ along two different abscissae Fig. 38 Damage sensing in conductive tissues: Layout of the electromechanical problem considered and arrangement of current input leads and potential probe leads for potential drop crack length determination oriented material points. The governing equations of the model for both elasticity and diffusion-based problems have been derived using a variational formalism. The obtained elasticity equations are not affected by internal restrictions involving the number of independent material constants. An implicit implementation strategy based on a meshfree approach has been also detailed together with analytical expressions of the equivalent stiffness operators. Particular attention has been given to establish a general micro-macro moduli correspondence between the continuum-molecular and classical continuum physics material parameters, in such a way that the micromoduli and their corresponding angular functions reproduce the overall anisotropy of the material. This theoretical scheme couples the intuitive simplicity of a purely pairwise description of (generalized) deformation and constitutive laws (removing Cauchy relations restrictions in elasticity), with the mathematical formalism of a continuum formulation for anisotropic materials. Moreover, the presence of oriented material points in the mechanical formulation allows to establish a natural relation with polar continua, as well as with real lattice-like and full-discrete systems. The accuracy of the method has been assessed by several numerical experiments, whose results have been compared with both analytical solutions and experimental data. Further studies are needed to extend the conceived anisotropic continuum-molecular model with oriented material points including multi-body potentials and other physical behaviors/physical problems such as piezoelectricity, Cosserat elasticity, cohesive fracture, plastic deformations, advection-diffusion and wave propagation.