PSI - Issue 72

Sreten Mastilovic / Procedia Structural Integrity 72 (2025) 538–546

540

Hence, each particle’s motion is uniquely determined by its initial conditions. In the present study, particles interact through a simple central-force law, which fully defines the strain energy density function based on their positions r . The system of differential equations of motion (2) is usually approximated by a corresponding finite difference scheme and solved using any of established numerical algorithms, Mastilovic (2022), Krajcinovic (1996) and Allen et al. (1987). For example, the Verlet equation for computing the next position of particle i can be written in a modified form:             j i i ij i j i i i i m f r r t t r t η r t r t t r t , ) ( ) ( ( ) ) ( 2    (3) adapted for quasistatic loading by including a damping coefficient η . This time-reversible scheme is commonly employed in PD with a damping coefficient η = 0.95 to dissipate kinetic energy and reach a quasi-equilibrium state reasonably fast. Due to the characteristics of interparticle forces, the magnitude of the central force exerted by particle j on its neighboring particle i , f ij (and vice versa), as well as the total force acting on particle i from all its first-nearest neighbors, F i , can be determined using the following expressions:       j ij ij ij i ij ij ij r f f r r F f ,  (4) To improve computational efficiency, it is common practice in MD, Allen et al. (1987), to create and maintain a neighbor list for each particle. Between periodic updates of the neighbor list, the program limits interaction checks to particles already on the list, rather than evaluating every particle pair in the system—significantly reducing computation time. The time necessary to evaluate interactions among all N particles using Verlet algorithm (3) is proportional to N2, Allen et al. (1987). In PD models, the interparticle potential must be specified from the outset. It is important to emphasize that, together with the spatial and temporal scales of the simulation, this potential fundamentally distinguishes PD from MD. Nonetheless, both approaches share a key aspect: selecting an appropriate potential is crucial for two main reasons. First, the accuracy and relevance of the potential directly influence the quality of the simulation results. Second, the potential’s complexity affects the program’s computational efficiency and execution time. The hybrid potential, widely used in particle models, Mastilovic (2022 and 2008), to simulate the behavior of brittle and quasi-brittle materials at the mesoscale, is expressed as follows:           1 1 1 , 2 2       a a r r r f r k r r k r r r       (5) 2 This potential combines the Hookean (5)₁ and Born-Mayer (5)₂ forms, where the superscripts a and r denote the attractive and repulsive branches of the interaction, respectively. The bond stiffness kij has an average value set according to the desired elastic modulus Mastilovic (2022 and 2008), while the adjustment parameter B controls the steepness of the repulsive wall and can be chosen, in principle, using the ballistic equation of state, Mastilovic (2022). The force–distance relationship (5) is designed to capture key features of the deformation behavior typical of quasi-brittle materials, including brittle tensile response, increased shock wave velocity, and decreased compressibility under rising pressure Mastilovic (2022 and 2008). It is important to emphasize that all parameters of potential (5) are identifiable a priori using a macroscopic, top-down approach rather than from microscopic-scale material behavior (which is also an option). 2 1 0 0 ij ij ij 0 ij ij ij ij ij ij ij ij ij (5) 1   ij r      1  rij    ij      1  rij    1 2 , 1 2 B B 2 0 ij ij 1 2 0 ij ij                      ij ij r ij f r ij r r r e k r r e k r    B B B 