Flaw-tolerance of nonlocal discrete systems and interpretation according to network theory

Discrete systems are modeled as a network of nodes (particles, molecules, or atoms) linked by nonlinear springs to simulate the action of van der Waals forces. Such systems are nonlocal if links connecting non-adjacent nodes are introduced. For their topological characterization, a nonlocality index (NLI) inspired by network theory is proposed. The mechanical response of 1D and 2D nonlocal discrete systems is predicted according to finite element (FE) simulations based on a nonlinear spring element for large displacements implemented in the FE programme FEAP. Uniaxial force-displacement responses of intact and defective systems (with links or nodes removed) are numerically simulated. Strain localization phenomena, size-scale effects and the ability to tolerate defects are investigated by varying the degree of nonlocality.


INTRODUCTION
onlocal continuum theories based on gradient models, integral formulations or fractional calculus have been widely explored in mechanics to describe long-range interactions (see e.g.[1][2][3][4][5][6][7][8][9][10], among others).At the same time, discrete systems composed of particles or molecules have been proposed in the physics community to analyze the behaviour of materials at the nano or micro scales.Lattice beam models [11], albeit suffering from mesh dependency due to the local nature of the bonds/links, have been extensively used to simulate the meso-scale behavior of concrete.Efforts accounting for nonlocal effects, such as the three-dimensional Born model [12], have been proposed to study the distribution of broken bonds within a homogeneous discrete mechanical system.With the progress in computer technology, wide tridimensional discrete systems can nowadays be modelled by molecular dynamics (MD), accounting for nonlinear interatomic potential laws and nonlocal interactions among the discrete molecules.Attempts to couple MD with FEM have been explored in [13], with the intent to perform multi-scale analyses where discrete systems are connected to continuum finite elements.In the present study, the ability of nonlocal molecular systems to tolerate flaws is investigated.To this aim, a nonlinear spring element for large displacements whose constitutive relation is ruled by van der Waals forces is implemented in the N finite element analysis programme FEAP.The tensile and compressive responses of 1D and 2D discrete systems of molecules with or without flaws are numerically simulated, by varying the range of nonlocal interactions defined by a novel nonlocality index NLI inspired by graph and network theories.For each system, the distribution of the internal forces at different load levels is analyzed in relation to the topological properties of the underlying network, in order to understand the force redistribution mechanisms in the presence of defects.

FINITE ELEMENT FORMULATION
onsidering a discrete system composed of atoms or molecules subjected by nonlinear interactions governed by generalized Lennard-Jones potentials, graph theory can be efficiently invoked for their characterization.Atoms can be topologically represented by nodes and the nonlinear interactions between them can be modelled by a set of nonlinear springs.In this framework, each link can be implemented as a finite element defined by two nodes whose coordinates in the initial global undeformed reference system are: Each i-th node (i=1,2) has two translational degrees of freedoms in case of 2D problems, u i and v i , collected in the element displacement vector u: The nodal coordinates vector in the updated configuration is given by: A local reference system is now introduced, since the constitutive relation is usually defined in this frame.Hence, the normal vector n pointing from node 1' to node 2' and the vector t perpendicular to the link are introduced, see Fig. 1.The origin of the local reference system is in the node 1'.Such a frame is rotated by an angle  with respect to the axis of the global reference system.The following rotation matrix R can be introduced: Figure 1: Reference coordinate systems for the large displacements analysis.
The gap vector g collecting the relative opening and sliding displacements between the nodes 1' and 2' is given by the product between a matrix operator L and the position vector x containing the coordinates of the element nodes in the updated configuration: where the matrix L takes the form:

C
The gap vector between the nodes in the local reference system is obtained by the product between the rotation matrix R and the gap vector g: The constitutive law used here stems from a generalization of the Lennard-Jones (L-J) potentials: where  is the depth of the potential well, r 0 is the distance at which the potential between two neighbouring atoms is equal to zero, , are the parameters defining the shape of the L-J potential.Keeping the exponents , as free parameters was suggested in [14][15][16] to bypass the rapid decay of the law after the equilibrium point and being able compute non-bonded interactions between nodes also at longer distances, or for coarse-grained representations.
The corresponding displacement-force relation can be obtained by differentiating the L-J potential with respect to the normal opening displacement g n : In the proposed formula, whose graphical trend is depicted in Fig. 2  The force in Eq. ( 9) directed along the normal vector n is collected in the local force vector, F loc  (F LJ ,0) T .Starting from the Principle of Virtual Work where each element contribution is related to the scalar product between the force vector and the gap vector as for a classical interface element [17,18], its virtual variation and fully consistent linearization have to be performed in order to derive the residual vector p and the stiffness matrix K for the application of the Newton-Raphson iterative procedure.After a lengthy calculation, which is here omitted for the sake of brevity, we obtain: where C is the tangent operator stemming from the consistent linearization of the local force vector Floc with respect to the local gap vector gloc: whereas the term R / u related to the differentiation of the rotation matrix with respect to the displacements vector is a third order tensor.This mathematical formulation has been implemented in the Finite Element Analysis Program FEAP [19].A flow chart of the operations is shown in Tab. 1.

NUMERICAL INVESTIGATION: MONODIMENSIONAL MODEL
he first step of the study is the understanding of the nonlinear behaviour of a L-J nonlocal system of atoms for different degrees of nonlocality.Let us consider a 1D system of nodes representing atoms equally spaced by a length l 0 along a straight line.A local system would be represented by the case where only local interactions take place between each node and its nearest neighbours (Fig. 4(a)).This is just a mechanical system of nonlinear springs in series, where the removal of a link (failure of one spring) leads to the failure of the entire chain (Fig. 4(b)).Long range interactions can be modelled by introducing additional links between a generic node and the nodes located at distances 2 l 0 , 3l 0 ,… and so forth.Various systems with different degree of nonlocality can be generated by changing the cut-off distance for nonlocal interactions.For comparison purposes, it is useful to introduce a dimensionless number NLI (nonlocality index) synthesizing the topological properties of the system: where l i is the length of the i-th link between two connected nodes and l 0 is the length of the local link between two adjacent nodes.Therefore, NLI=1 corresponds to a local system, whereas NLI>1 refers to systems with long range interactions.
As a first example, let us consider the uniaxial tension test for a system composed of 5 nodes and NLI=2 (Fig. 5  For a given set of stiffnesses, it is possible to compute analytically the relation between the displacement  of the whole system and the displacement  i of each single spring, satisfying the equilibrium and the congruence: The relations between the spring forces and the total imposed displacement  are: These analytical results provide some qualitative information on the solution of the fully nonlinear system.In general, comparing  1 and  2 in Eq.( 13), we expect  1 >  2 , which is a main difference with respect to a purely local system where all the local springs experience the same elongation.Moreover, if a spring is removed, then the system has still an equilibrium configuration, whereas a purely local system would fail.The solution of the fully nonlinear problem is discussed in the next subsections.

Monodimensional model without defects: the effect of nonlocality and system size
Let us consider the discrete system in Fig. 6 with NLI=1, 2 or 3.In the local system characterized by springs in series, the forces and the elongations are the same for all the springs.Each spring presents an increasing branch until the achievement of the maximum force and then a softening, according to the L-J constitutive equation (9).In case of NLI=2, the situation is much more complex.Due to the symmetry, it is sufficient to comment on the behaviour of the springs 1, 2, 5 and 6.The forces in such springs vs. the total imposed displacement  are shown in Fig. 7(a).The positions of the nodes of each spring vs.  are illustrated in Fig. 7(b).From these diagrams we note that springs 1 and 5 have elongations monotonically increasing during the simulation.The force level of the local spring 1 is much higher than that of the nonlocal spring 5.This is due to the fact the spring 5 is connecting nodes that are at an initial distance twice larger than that of the spring 1. Springs 2 and 6, after an initial positive relative displacement between their connected nodes, then experience a progressive reduction of their elongations with a consequent elastic unloading.A similar response takes place for the system with NLI=3.Springs 1, 5 and 6 have an elongation monotonically increasing with time, whereas springs 2 and 6 experience an elastic unloading after an initial increased separation.
A comparison between the force-displacement curves of the various systems depending on NLI is shown in Fig. 8(a).The nonlocality increase the peak force observed in case of NLI=1.Moreover, a gain in the capacity of tolerating defects takes place, as we will quantify in the next subsection.The difference between the responses of the systems with NLI=2 and NLI=3 is quite small, in line with the observation that nonlocal interactions are rapidly decaying by increasing the nodal distance, see Fig.For a given degree of nonlocality, e.g., NLI=2, the effect of the size of the system can be explored by increasing the number of nodes at a distance l0 for each pair (see Fig. 8(b)).As a trend, the peak force does not change, whereas the post-peak response becomes more and more brittle by increasing h.

Monodimensional model with defects
Let us consider a system with 5 nodes ( h  4 l 0 ) and with a broken link (see Fig. 9(a)).This can be the case in nanoscopic systems when there is an atomic vacancy or a weakened interatomic potential due to dislocations.The problem is solved numerically under displacement control and a very complex behaviour of the various springs is observed.To comment on that, it is useful to analyze the force-nodal distance curve of each spring, Fig. 9(b).All the springs in series (1,3,4) are initially in the initial branch of the van der Waals forces and all of them experience an elastic unloading (see the direction of the arrows in Fig. 9(b)).Springs 5, 6 and 7 used to model long-range interactions are all in the softening branch of the van der Waals force diagram.The spring 7 experiences an increasing in the corresponding load level due to the shortening of springs 3 and 4. On the other hand, 5 and 6 increase their elongation monotonically in time.The global response, shown in Fig. 10(a) presents softening due to the heavy damage induced by the removal of a local link.However, the defect can be tolerated due to the nonlocality and the load-carrying capacity of the system is an increasing function of NLI, see Fig. 10(a).By increasing the system size and removing only the second spring in all the cases, the force is an increasing function of h as shown in Fig. 10(b), as the same was observed in case of a defect-free system.

Computational complexity
s is the monodimensional case, a bidimensional system can be modelled as a collection of nodes and links, thus forming a grid in the plane.An hexagonal disposition of nodes is adopted to simulate materials like graphene.Also in this case the NLI can be defined as the maximum cut-off radius for long range interactions and the length of the link for pure local connections.The numerical construction of the finite element meshes for numerical simulations requires the definition of the list of nodes and the element connectivity matrix.For the latter, a pre-processor written in MATLAB ® has been developed where, centering in the i-th node of the grid, all the other nodes within the cut-off radius defined by the product NLI l 0 are identified and spring elements are inserted.Clearly, by increasing NLI and the size of the system, the computational complexity drastically increases.As an illustrative example, let us consider square grids with l 0 =3.5Å and with lateral size 14x14 Å, 31x31 Å, and 70x70 Å, respectively.For each grid, the NLI can range from 1 up to a maximum value NLI max , given by the longest distance between two nodes belonging to the system.For the three system sizes considered, we have NLI max  6,12 and 29.The maximum number of links in case of NLI max is equal to 600 for the smallest system, 4,851 for the intermediate system, and 126,756 for the largest one.In dimensionless form, Fig. 11, the ratio between the number of links and the maximum value has a nonlinear dependency with the ratio between NLI and NLI max .The deviation from linearity observed for all the cases when NLI is approaching NLI max is due to the finite size of the system.A As far as the numerical simulations are concerned, the increase of NLI leads to a transition from a sparse connectivity matrix to a very dense one, see Fig. 12 where systems with NLI=1, 2 and 3 are compared.

Flaw-tolerance
In this section, the behaviour of the bi-dimensional system subjected to uniaxial tension in the vertical direction is investigated.We compare the case where all the links are present with other three cases presenting defects: (1) removal of a single node in the centre of the specimen; (2) removal of 2 not adjacent nodes; (3) removal of 2 adjacent nodes.Three different values of NLI are considered for each system: NLI=1 (local system), NLI=2, NLI=3.A sketch of the three configurations with defects are shown in Fig. 13 for the case with NLI=1.The global force vs. displacement curves for all the examined configurations are shown in Fig. 14.For a given value of NLI, we note that Case 3 leads to the higher load and stiffness reduction.Case 1 is the less critical one, whereas Case 2 provides an intermediate response.This trend is in line with the fact that the removal of two adjacent nodes is particularly critical and induces a significant number of removed links.The effect of NLI is particularly important.The system with NLI=1 is initially compressed, due to the fact that the nodes are initially at a distance slightly smaller than the equilibrium one given by the L-J potential.By increasing NLI, further links in a tension state are introduced and reduce this initial compressive stress state.Although the nodes removal has a higher impact on the systems characterized by higher nonlocality indices since a larger number of links is removed, the force is an increasing function of NLI, thus showing that nonlocality has a very important role for flaw-tolerance.

CONCLUSIONS
aking into consideration nonbonded interactions between nodes according to the L-J generalized potential, the influence of nonlocality on the mechanical response of discrete systems has been investigated.In the monodimensional case, we have observed a complex behaviour of the links depending on the nonlocality index NLI.In the case of a link absent or removed, a phenomenon analogous to strain localization observed in continuum systems has been evidenced.Results on bi-dimensional systems have also shown the important role exerted by the nonlocality on the load carrying capacity of the system and its ability to tolerate defects or vacancies.The proposed approach and its further extension to 3D can be applied to predict strength and stiffness of discrete systems as, e.g., 2D graphene layers of carbon-carbon bonds, or nanotubes, represented by a 3D network of carbon atoms.Moreover, a constitutive law based on L-J potential is definitely useful to study problems where interactions between two different materials are governed by nonbonded forces.This is for example the case of an interphase boundary interaction between collagen and graphene.Future developments regarding the implementation of bonded interactions, including not only stretching energy, but also bending and torsion energy contributions to describe physical links between atoms or molecules, are envisaged.
(a), the first term in brackets represents the repulsive force contribution between atoms (Pauli repulsion), while the second one describe the long-range attractive force (van der Waals force).Different force-displacement curves are shown in Fig.2(b) by varying the parameters of the model.In the numerical simulations discussed in the next sections, the selected parameters are   13 and   2.

Figure 2 :
Figure 2: (a) Lennard-Jones potential law: attractive and repulsive branches.(b) Different force-displacement forces by varying and  : the red curve refers to L-J potential in Fig. 2(a).The blue curve is the modified law used in this study.

Figure 4 :
Figure 4: (a) System with only local links.(b) Collapse of the system due to failure of just one link.
(a)).Four local links are present and three additional links due to long-range interactions are introduced.For a preliminary qualitative analysis of the system, a linearized version would consist in springs with different stiffnesses ki, depending on T LOOP until convergence of the Newton-Raphson scheme CALL to the element subroutine INPUT: X initial nodal coordinates, u nodal displacements COMPUTE: updated coordinates x=X+u rotation matrix R and gap vector g (Eq.5) gap vector in the local frame g loc (Eq.7) local force vector F loc and the tangent operator C (Eq. 11) OUTPUT: residual vector p and tangent stiffness matrix K (Eq.10) END LOOP the distance between the connected nodes.Due to the symmetry of the problem, only half model can be analyzed (Fig.5(b)).

Figure 8 :
Global force vs. global displacement curves.(a) The effect of NLI for a system with a size h  4 l 0 .(b) The effect of the size h of the system with NLI=2.

Figure 9 :Figure 10 :
Figure 9: (a) Geometry of the 5 nodes system without the link no. 2. (b) Force-displacement diagram of the springs.

Figure 11 :
Figure 11: Number of links vs. NLI in dimensionless form, for 3 systems with different sizes.

Figure 12 :
Figure 12: (a) Representation of the bidimensional hexagonal geometric network by varying NLI.(b) Connectivity matrices for the same NLI.

Figure 13 :
Figure 13: Representation of the bidimensional hexagonal network with different types of flaws (NLI=1).

Figure 14 :
Figure 14: Comparison of the force-displacement curves varying NLI, depending on the type of defect.

Table 1 :
Sequence of element operations.