# Particle Equilibrium Method for Crack Propagation Simulation

Authors: Bo-Bo Xiong; Bin Tian; Xiao-chun Lu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

Crack propagation simulation is constantly of great significance. This study presents the particle equilibrium method (PEM) to achieve this goal. PEM is based on the idealization of the problem domain as an assemblage of distinct particles, which release interaction forces to their surrounding particles. Once the distance between two particles exceeds the extreme distance, a permanent crack will emerge between them. The function of interaction force describes the mechanical response of elastic and plastic materials. The calculated structure attains its final state to the external load when all the particles reach the equilibrium condition of force. A calculation program is developed based on the proposed method’s theory and is applied to three numerical examples with reliable calculation results.


---

## Body

## 1. Introduction

Crack propagation is one of the most common phenomena in structures, and this remarkable feature is consistently associated with the damage or destruction of structures. Thus, crack propagation simulation is of great significance, and several methods have been applied to address this problem. The finite element method (FEM) is a well-known numerical analysis method that has been extensively used in predicting crack paths [1]. However, FEM requires several remeshing processes of the finite element model to represent arbitrary and complex crack paths, which is slightly difficult; moreover, automatic remeshing can result in highly distorted elements, which worsens the performance of the FEM [1, 2]. To improve the FEM in terms of fracture simulation, extended finite element method (XFEM) [3], fractal finite element method [4], node-based smoothed extended finite element method (NS-XFEM) [5], and edge-based smoothed finite element method (ESFEM) [6] were proposed. However, these FEM versions are mesh-based approximation methods and have certain inherent limitations [7].Thus, several meshless methods have been developed. Meshless methods discretize the problem domain only with a set of unstructured points or nodes without using a predefined element mesh. In contrast to the “element” concept of the FEM, nodal connectivity in meshless methods is enforced with the “influence-domain” concept [1]. This feature makes the meshless method suitable for modeling fracture propagation. Meshless methods that commonly use the weak form of Galerkin can be classified as approximation meshless methods and interpolation meshless methods [1]. Approximation meshless methods were proposed first, including smooth particle hydrodynamics (SPH) [8–12], diffuse element method (DEM) [13], and element-free Galerkin method (EFG) [14, 15]. Although these meshless methods have been widely used, they have difficulty in imposing essential and natural boundary conditions. Consequently, interpolation meshless methods, such as nature element method (NEM) [16], point interpolation method (PIM) [7], radial PIM (RPIM) [17], and nature neighbor RPIM [18–20], have been developed. These methods have improved the solution for the problem.The present study proposes an approach for modeling structural responses to external loads, especially for simulating the crack growth phenomenon in brittle materials. This method begins by discretizing the problem domain as a set of particles, which interact solely with their surrounding particles through an interaction force. The concept of “element” is used in this method to record the relative position relationship between particles and their surrounding particles. This interaction force appears as a repulsion when two particles are relatively near each other and as attraction when the distance between two particles is far. In this study, the interaction force is expressed by a function, and different functions are established to describe the mechanical properties of elastic and plastic materials. If the distance between two particles exceeds the extreme distance, then a crack is generated between them and their interaction force would be interrupted. Thus, the crack propagation of the material can be naturally simulated and the damage rate of the material from the cracks can be obtained simultaneously. When all particles reach the equilibrium condition of forces, the system attains a stable state, and this is the final response to the external load. On the basis of the principle, a calculation program is developed and applied on three examples.In the following sections, the basic principle of the proposed method is introduced. Second, the development of the numerical calculation program of the proposed method is presented. Finally, the performance of the proposed method is evaluated using three numerical examples.


## 2. Basic Principle of the Particle Equilibrium Method (PEM)

Suppose that a target object occupies a reference configuration in region Ω and this object consists of many particles that are not mutually in contact. Any particle may undergo the effect of various loads, such as the interaction force among the surrounding particles, as well as the body (gravity or seepage force) and boundary loads, if this particle is in boundaryΓ. The particles within the target object may develop displacement under the applied loads. If all particles reach balance, then this displacement is the final response of the object to the external load.


### 2.1. Discrete Particle Systems

The target object should be divided into many nontouching particlesxi(xi∈Ω,i=1,2,3…n), where n is the total number of particles that constitute the target object, to determine the response of the object under an external load. The discretization procedure is the same as the conventional FEM. Thus, PEM can use the data of the finite element mesh, which is easily acceptable to numerous users. The divided particles split the object volume; however, the volume concentrates on the particle. As shown in Figure 1, the volume Vi of the particles xi(xi∈Ω) is obtained as follows:(1)Vi=∑j=1mλijVij,wherem is the total number of elements that contain particle xi; Vij is the volume of elements containing xi; and λij is the weight of node xi to Vij, which is determined by the geometry of the grid.Figure 1
Interaction force in the calculation system.


### 2.2. Interaction Force

The definition of the interaction force among different particles is one of the key characteristics of PEM. For particlexi(xi∈Ω), only its surrounding particles release an interaction force, and all the surrounding particles belong to the elements containing particle xi. The interaction force fij between particle xi and one of its surrounding particle xj(j=1,2,3…k) can be defined as follows:(2)fij=fs,Vi,Vj,where k is the total number of surrounding particles to xi; s is the distance between xi and xj; and Vi and Vj are the volumes of xi and xj, respectively.As shown in Figure2, the interaction force fij is similar to the van der Waals force in molecular dynamics; that is, when s is smaller than se (se is the balance distance), the interaction force fij takes the form of repulsion, which prevents the infinite contraction of an object and results in increased repulsion. On the contrary, when s is larger than se, the interaction force fij manifests as attraction, which maintains a solid object with a certain shape; when s is relatively large, the attraction force becomes near zero. Balance distance se is a special distance, in which the repulsion and attraction between two particles are the same in value and the interaction force fij is zero.Figure 2
Interaction force function.In PEM, the constitution functionf() is related to the material properties and the stress states of materials. Thus, material tests are necessary when determining the exact constitution function f(). This special form of constitution function describes the relationship between force and displacement in a new manner, which avoids derivatives of the element displacement. The displacement field associated with this method can be discontinuous. Therefore, the cracking behavior of materials and the slip movements among phases are easy to describe using PEM.


### 2.3. Fracture Criterion

As shown in Figure2, whens is relatively large, the attraction force becomes near zero. This special numeric feature is not consistent with the actual computation. Hence, cutoff distance sc is introduced. In the actual calculation, when the distance of two particles exceeds sc, the interaction force fij is zero and a crack appears between particles xi and xj.Therefore, the aforementioned fracture criterion indicates that this PEM is convenient for and easy to use in describing the fracture characteristics of materials, such as crack initiation, propagation, and arrest without a definite complex crack stress or mesh rezoning. Furthermore, PEM can efficiently stimulate multiple cracking in three-dimensional spaces, which enables it to have a broad application prospect.The damage rate of materials can be defined as the number of interaction forces that cannot transmit attraction to the total number of interaction forces in the computational domain. This damage rateϕΩ is expressed as(3)ϕΩ=1-NbNt,where Nb expresses the number of interaction forces that cannot transmit attraction in region Ω and Nt is the total number of interaction forces in the calculation domain. The cracking and damage rates of materials can be considered natural when analyzing their damage failure process using the proposed method.


### 2.4. Simulation of Existing Cracks

The cracks in the target object that exist originally or are generated during the process of responding to the external load will considerably change the mechanical properties of the materials. The performance of materials with cracks should be defined. A crack is supposed to exist between particlesxi and xj, as shown in Figure 3. The interaction force between xi and xj and other interaction forces involved in the crack should be redefined according to the actual distance between them. In PEM, once the interaction force between two particles is interrupted because of the exceeding distance or the original crack, this pair of interaction forces ceases to transmit attraction; hence, crack impairment is an irreversible process. However, if the distance between these two particles is less than se during the following deformation, then repulsion exists between them and produces friction when shearing deformation occurs (Figure 3). This friction constantly occurs along the normal direction of the line that connects two particles, and this friction is assumed to act on the two particles precisely. The exact value of this friction is obtained as follows:(4)T=fij·μ,where μ is the friction factor. The proposed method aims to study the performance of brittle materials and thus regards the dislocation in the cracks as relatively small, thereby avoiding the definition of a complicated contacting relationship.Figure 3
Schematic of cracks.Crack propagation is the only factor that results in the failure of the target object in PEM. When part of the calculation object separates from itself completely because of the cracks, the calculation object exceeds its limit-bearing capacity and the computational process is terminated.


### 2.5. Constitutive Model

For the convenience of application, two of the most well-known constitutive models are introduced. First, the widely used linear elastic constitutive model is easy to use and represents most characteristics of one material, which works in this elastic stage in most cases. The relationship between the interaction force and distance in this model is plotted in Figure4, which indicates model complexity. This model describes the elastic stage and damage, as expressed as follows:(5)fij=V1.V2se2k1s-se,s≤smaxk2s-sc,smax<s≤sc0,s>sc,where k1 and k2 are the line slopes in different stages, as shown in Figure 4. All these parameters can also be obtained through experiments.Figure 4
Schematic of elastic constitutive model.Second, a plastic constitutive model should be defined, considering that the external load may cause irreversible plastic deformation in many cases. Figure5 depicts the relationship between the interaction force fij and distance s. The most notable feature of this model is that balance distance se is variable because of plastic deformation. Consequently, the material produces inelastic, irrecoverable deformation under excessive pressure or tension. Therefore, this plastic constitutive model can properly describe a macroscopically mechanical behavior. Equation (6) feasibly expresses the interaction force fij in this plastic model, where b, n, and m are related to the material properties and can be obtained through the experiment.(6)fij=V1.V2se2s-seasesn-bsesms<scFigure 5
Schematic of conventional plastic constitutive model.


### 2.6. Governing Equation

The external load that worked on particlexi in the discrete system (xi∈Ω) can be described as pi, and the interaction force that acts on it can be expressed as(7)Fi=∑j=1mfij,where m is the total number of elements that contain node xi and fij is the interaction force between particle xi and its surrounding particle xj(j=1,2,3…m). Therefore, the balance equation of particle xi can be expressed as(8)pi=∑j=1mfij.In (8), the displacements of particles are the basic unknown variables. For the elastic constitutive model, this equation can be rewritten as(9)pix=∑j=1mV1.V2se2kd-sesix-sjxdpiy=∑j=1mV1.V2se2kd-sesiy-sjyd,where d is the distance between particles xi and xj. When all particles reach the equilibrium condition of forces (as expressed in (10)), the system attains a stable state, which is the final response to the external load.(10)∑j=1nFi=P,onΩ,where P is the external load vectors and n is the total number of the discrete particles.An explicit solution technique is used in the proposed method to track the crack propagation process. The dynamic equilibrium equation of the particlexi(xi∈Ω) can be written as(11)ami=pi-Fi,where a is the acceleration vector of the particle and mi is the volume of the particle. According to Newton’s law, a can be expressed as(12)a=∂2s∂t2=st¨.The central difference method is used in this method, and (11) can be expressed as follows:(13)st¨=pi-Fimi,(14)s˙t+1/2=s˙t-1/2+Δt·st¨,(15)st+1=st+Δt·s˙t+1/2,where s, s˙, and s¨ are the displacement, velocity, and acceleration of the particle, respectively; and Δt determines the computation speed. To guarantee the convergence of the calculation, (9) should be modified as(16)s˙t+1/2=αs˙t-1/2+Δt·st¨0<α<1.


## 2.1. Discrete Particle Systems

The target object should be divided into many nontouching particlesxi(xi∈Ω,i=1,2,3…n), where n is the total number of particles that constitute the target object, to determine the response of the object under an external load. The discretization procedure is the same as the conventional FEM. Thus, PEM can use the data of the finite element mesh, which is easily acceptable to numerous users. The divided particles split the object volume; however, the volume concentrates on the particle. As shown in Figure 1, the volume Vi of the particles xi(xi∈Ω) is obtained as follows:(1)Vi=∑j=1mλijVij,wherem is the total number of elements that contain particle xi; Vij is the volume of elements containing xi; and λij is the weight of node xi to Vij, which is determined by the geometry of the grid.Figure 1
Interaction force in the calculation system.


## 2.2. Interaction Force

The definition of the interaction force among different particles is one of the key characteristics of PEM. For particlexi(xi∈Ω), only its surrounding particles release an interaction force, and all the surrounding particles belong to the elements containing particle xi. The interaction force fij between particle xi and one of its surrounding particle xj(j=1,2,3…k) can be defined as follows:(2)fij=fs,Vi,Vj,where k is the total number of surrounding particles to xi; s is the distance between xi and xj; and Vi and Vj are the volumes of xi and xj, respectively.As shown in Figure2, the interaction force fij is similar to the van der Waals force in molecular dynamics; that is, when s is smaller than se (se is the balance distance), the interaction force fij takes the form of repulsion, which prevents the infinite contraction of an object and results in increased repulsion. On the contrary, when s is larger than se, the interaction force fij manifests as attraction, which maintains a solid object with a certain shape; when s is relatively large, the attraction force becomes near zero. Balance distance se is a special distance, in which the repulsion and attraction between two particles are the same in value and the interaction force fij is zero.Figure 2
Interaction force function.In PEM, the constitution functionf() is related to the material properties and the stress states of materials. Thus, material tests are necessary when determining the exact constitution function f(). This special form of constitution function describes the relationship between force and displacement in a new manner, which avoids derivatives of the element displacement. The displacement field associated with this method can be discontinuous. Therefore, the cracking behavior of materials and the slip movements among phases are easy to describe using PEM.


## 2.3. Fracture Criterion

As shown in Figure2, whens is relatively large, the attraction force becomes near zero. This special numeric feature is not consistent with the actual computation. Hence, cutoff distance sc is introduced. In the actual calculation, when the distance of two particles exceeds sc, the interaction force fij is zero and a crack appears between particles xi and xj.Therefore, the aforementioned fracture criterion indicates that this PEM is convenient for and easy to use in describing the fracture characteristics of materials, such as crack initiation, propagation, and arrest without a definite complex crack stress or mesh rezoning. Furthermore, PEM can efficiently stimulate multiple cracking in three-dimensional spaces, which enables it to have a broad application prospect.The damage rate of materials can be defined as the number of interaction forces that cannot transmit attraction to the total number of interaction forces in the computational domain. This damage rateϕΩ is expressed as(3)ϕΩ=1-NbNt,where Nb expresses the number of interaction forces that cannot transmit attraction in region Ω and Nt is the total number of interaction forces in the calculation domain. The cracking and damage rates of materials can be considered natural when analyzing their damage failure process using the proposed method.


## 2.4. Simulation of Existing Cracks

The cracks in the target object that exist originally or are generated during the process of responding to the external load will considerably change the mechanical properties of the materials. The performance of materials with cracks should be defined. A crack is supposed to exist between particlesxi and xj, as shown in Figure 3. The interaction force between xi and xj and other interaction forces involved in the crack should be redefined according to the actual distance between them. In PEM, once the interaction force between two particles is interrupted because of the exceeding distance or the original crack, this pair of interaction forces ceases to transmit attraction; hence, crack impairment is an irreversible process. However, if the distance between these two particles is less than se during the following deformation, then repulsion exists between them and produces friction when shearing deformation occurs (Figure 3). This friction constantly occurs along the normal direction of the line that connects two particles, and this friction is assumed to act on the two particles precisely. The exact value of this friction is obtained as follows:(4)T=fij·μ,where μ is the friction factor. The proposed method aims to study the performance of brittle materials and thus regards the dislocation in the cracks as relatively small, thereby avoiding the definition of a complicated contacting relationship.Figure 3
Schematic of cracks.Crack propagation is the only factor that results in the failure of the target object in PEM. When part of the calculation object separates from itself completely because of the cracks, the calculation object exceeds its limit-bearing capacity and the computational process is terminated.


## 2.5. Constitutive Model

For the convenience of application, two of the most well-known constitutive models are introduced. First, the widely used linear elastic constitutive model is easy to use and represents most characteristics of one material, which works in this elastic stage in most cases. The relationship between the interaction force and distance in this model is plotted in Figure4, which indicates model complexity. This model describes the elastic stage and damage, as expressed as follows:(5)fij=V1.V2se2k1s-se,s≤smaxk2s-sc,smax<s≤sc0,s>sc,where k1 and k2 are the line slopes in different stages, as shown in Figure 4. All these parameters can also be obtained through experiments.Figure 4
Schematic of elastic constitutive model.Second, a plastic constitutive model should be defined, considering that the external load may cause irreversible plastic deformation in many cases. Figure5 depicts the relationship between the interaction force fij and distance s. The most notable feature of this model is that balance distance se is variable because of plastic deformation. Consequently, the material produces inelastic, irrecoverable deformation under excessive pressure or tension. Therefore, this plastic constitutive model can properly describe a macroscopically mechanical behavior. Equation (6) feasibly expresses the interaction force fij in this plastic model, where b, n, and m are related to the material properties and can be obtained through the experiment.(6)fij=V1.V2se2s-seasesn-bsesms<scFigure 5
Schematic of conventional plastic constitutive model.


## 2.6. Governing Equation

The external load that worked on particlexi in the discrete system (xi∈Ω) can be described as pi, and the interaction force that acts on it can be expressed as(7)Fi=∑j=1mfij,where m is the total number of elements that contain node xi and fij is the interaction force between particle xi and its surrounding particle xj(j=1,2,3…m). Therefore, the balance equation of particle xi can be expressed as(8)pi=∑j=1mfij.In (8), the displacements of particles are the basic unknown variables. For the elastic constitutive model, this equation can be rewritten as(9)pix=∑j=1mV1.V2se2kd-sesix-sjxdpiy=∑j=1mV1.V2se2kd-sesiy-sjyd,where d is the distance between particles xi and xj. When all particles reach the equilibrium condition of forces (as expressed in (10)), the system attains a stable state, which is the final response to the external load.(10)∑j=1nFi=P,onΩ,where P is the external load vectors and n is the total number of the discrete particles.An explicit solution technique is used in the proposed method to track the crack propagation process. The dynamic equilibrium equation of the particlexi(xi∈Ω) can be written as(11)ami=pi-Fi,where a is the acceleration vector of the particle and mi is the volume of the particle. According to Newton’s law, a can be expressed as(12)a=∂2s∂t2=st¨.The central difference method is used in this method, and (11) can be expressed as follows:(13)st¨=pi-Fimi,(14)s˙t+1/2=s˙t-1/2+Δt·st¨,(15)st+1=st+Δt·s˙t+1/2,where s, s˙, and s¨ are the displacement, velocity, and acceleration of the particle, respectively; and Δt determines the computation speed. To guarantee the convergence of the calculation, (9) should be modified as(16)s˙t+1/2=αs˙t-1/2+Δt·st¨0<α<1.


## 3. Numerical Scheme of PEM

Although PEM’s basic theory is introduced in Section2 as a numerical calculation method, the exact calculation procedure is an indispensable part through which this method can solve experimental problems. PEM 1 is the first computer code based on the proposed method. This computer code is written in C#, and some of the main technical points have been previously illustrated in this study. The method is currently in the research and development stage, and a new version with additional practical functions will be created in the future. For convenience, the procedure of the proposed method is systematically reorganized as follows.Step 1.
Input the calculation data, including geometric information (grid and node data), load, constraint, material information, and controlling parameters.Step 2.
Calculate the interaction force of the particles based on the current displacement (using (9)). Then, calculate the acceleration of the particles under the current load (using (13)).Step 3.
Calculate the new location of the particles (using (15) and (16)), which is the kernel of this method.Step 4.
Calculate the current damage rate to determine whether the calculating stepΔt should be changed (using (4)). If excessive cracks emerge under this incremental load, repeat Step 2 and recalculate the location of the target object with a relatively small calculating step Δt. Only through this manner will the computational process be stable and the crack prorogation process be captured accurately.Step 5.
Save the intermediate result, such as the coordinates of the particles under the current calculating step and the damage rate, which can present the dynamic response process.Step 6.
Calculate the velocity of the particles. If the velocity is close to zero for several calculating steps, then consider the calculation process convergent. Otherwise, go back to Step2.Step 7.
Plot the displacement nephogram and crack growth path of every step.The above procedure is illustrated in a flowchart in Figure6.Figure 6
Calculation flowchart of the proposed method.


## 4. Numerical Tests

To verify the feasibility and accuracy of the proposed method, three numerical examples are presented in this section.


### 4.1. Cantilever Beam under Uniformly Distributed Loads

A cantilever beam under distributed loads is analyzed in this example. As shown in Figure7(a), an 8-m long and 1-m wide cantilever beam is under distributed loads (q=1000N/m). The materials of this cantilever beam are described using the linear elastic constitutive model, and the relevant material parameters are set as k1=3×108P,k2=-3×109Pa,smax=1.045se and sc=1.05se [22]. As shown in Figure 7(b), the cantilever beam is divided into 891 particles. FEM (ABAQUS software is used here) and the proposed method are employed to analyze this problem, and the results are displayed in Figure 8. The comparison of the nephograms of the calculated displacements indicates that the distributed displacements of the cantilever beam, which are calculated by the two methods, are the same and the relevant values of the displacement are similar. Given the analytical solution to the vertical displacement of the cantilever beam under uniformly distributed loads (Figure 7(b)), Figure 9 shows the comparison of the vertical displacements of the cantilever beam calculated by ABAQUS, PEM, and the analytical formula; this comparison indicates that the results are consistent with one another. Thus, the proposed method can analyze the response of the structure under external loads.Figure 7
Schematic of the cantilever beam under uniformly distributed loads.
(a) (b)Figure 8
Displacement nephogram: (a) U1 (calculated by the proposed method), (b) U2 (calculated by the proposed method), (c) U1 (calculated by ABAQUS), and (d) U2 (calculated by ABAQUS).
(a) (b) (c) (d)Figure 9
Vertical displacement of cantilever beam.


### 4.2. Simulation of Brazilian Splitting Test

This example simulates the cracking course of the disc under pressure to verify the capability of the proposed method in simulating crack propagation. As shown in Figure10(a), this disc has a straight-through crack in the center; the angle between this crack and the vertical direction is α=25°; the crack length is a=0.8,r=0.4m (r is the radius of the disc); and the pressure acting on each end of the disc is F=2×106N. The material parameters of this disc are the same as those of the cantilever beam in the previous example.Figure 10
Disc-cracking test: (a) schematic of disc-cracking test, (b) crack paths calculated by the proposed method, and (c) crack paths from the experimental test [21].
(a) (b) (c)Figure10(b) presents the distribution of fractures on the disc after the calculation using the proposed method, whereas Figure 10(c) shows the actual cracks on the disc under the experimental conditions [21]. The comparison of Figures 10(b) and 10(c) implies that the crack patterns calculated by the proposed method are consistent with the test results. This example shows the potential of PEM in simulating crack propagation.


### 4.3. Three-Point Bending Beam of Concrete

In this example, the proposed method is used to study the cracking of the concrete structure. The calculation model of a three-point bend beam made of concrete is presented in Figure11. To reflect the composition characteristics of the concrete, aggregate, cement mortar, and the interface layer between the aggregate and the cement mortar of the concrete material are endowed with different properties (as shown in Figure 11). All those three-phase compositions are described using the linear elastic constitutive model, and the relevant material parameters are set as follows [23, 24]: k1A=6×1010P, k2A=-6×1011Pa,smaxA=1.045se, and scA=1.05se; k1C=4×1010P, k2c=-4×1011Pa, smaxc=1.045se, and scc=1.05se; and k1I=2.5×1010P, k2I=-2.5×1011Pa,smaxI=1.045se, and scI=1.05se, where superscripts A, C, and I represents the aggregate, cement mortar, and the interface layer between the aggregate and the cement mortar, respectively.Figure 11
Schematic of the three-point bending beam of concrete.Figure12 presents the crack path on the three-point bend beam after the calculation using the proposed method. Most of the cracks occur along the interface layer between the aggregate and the cement mortar, which corresponds to ordinary rules [25, 26]. This example demonstrates that the proposed method has the potential of becoming a mature technology for simulating complex deformation and cracking behavior.Figure 12
Crack path of the three-point bending beam of concrete.


## 4.1. Cantilever Beam under Uniformly Distributed Loads

A cantilever beam under distributed loads is analyzed in this example. As shown in Figure7(a), an 8-m long and 1-m wide cantilever beam is under distributed loads (q=1000N/m). The materials of this cantilever beam are described using the linear elastic constitutive model, and the relevant material parameters are set as k1=3×108P,k2=-3×109Pa,smax=1.045se and sc=1.05se [22]. As shown in Figure 7(b), the cantilever beam is divided into 891 particles. FEM (ABAQUS software is used here) and the proposed method are employed to analyze this problem, and the results are displayed in Figure 8. The comparison of the nephograms of the calculated displacements indicates that the distributed displacements of the cantilever beam, which are calculated by the two methods, are the same and the relevant values of the displacement are similar. Given the analytical solution to the vertical displacement of the cantilever beam under uniformly distributed loads (Figure 7(b)), Figure 9 shows the comparison of the vertical displacements of the cantilever beam calculated by ABAQUS, PEM, and the analytical formula; this comparison indicates that the results are consistent with one another. Thus, the proposed method can analyze the response of the structure under external loads.Figure 7
Schematic of the cantilever beam under uniformly distributed loads.
(a) (b)Figure 8
Displacement nephogram: (a) U1 (calculated by the proposed method), (b) U2 (calculated by the proposed method), (c) U1 (calculated by ABAQUS), and (d) U2 (calculated by ABAQUS).
(a) (b) (c) (d)Figure 9
Vertical displacement of cantilever beam.


## 4.2. Simulation of Brazilian Splitting Test

This example simulates the cracking course of the disc under pressure to verify the capability of the proposed method in simulating crack propagation. As shown in Figure10(a), this disc has a straight-through crack in the center; the angle between this crack and the vertical direction is α=25°; the crack length is a=0.8,r=0.4m (r is the radius of the disc); and the pressure acting on each end of the disc is F=2×106N. The material parameters of this disc are the same as those of the cantilever beam in the previous example.Figure 10
Disc-cracking test: (a) schematic of disc-cracking test, (b) crack paths calculated by the proposed method, and (c) crack paths from the experimental test [21].
(a) (b) (c)Figure10(b) presents the distribution of fractures on the disc after the calculation using the proposed method, whereas Figure 10(c) shows the actual cracks on the disc under the experimental conditions [21]. The comparison of Figures 10(b) and 10(c) implies that the crack patterns calculated by the proposed method are consistent with the test results. This example shows the potential of PEM in simulating crack propagation.


## 4.3. Three-Point Bending Beam of Concrete

In this example, the proposed method is used to study the cracking of the concrete structure. The calculation model of a three-point bend beam made of concrete is presented in Figure11. To reflect the composition characteristics of the concrete, aggregate, cement mortar, and the interface layer between the aggregate and the cement mortar of the concrete material are endowed with different properties (as shown in Figure 11). All those three-phase compositions are described using the linear elastic constitutive model, and the relevant material parameters are set as follows [23, 24]: k1A=6×1010P, k2A=-6×1011Pa,smaxA=1.045se, and scA=1.05se; k1C=4×1010P, k2c=-4×1011Pa, smaxc=1.045se, and scc=1.05se; and k1I=2.5×1010P, k2I=-2.5×1011Pa,smaxI=1.045se, and scI=1.05se, where superscripts A, C, and I represents the aggregate, cement mortar, and the interface layer between the aggregate and the cement mortar, respectively.Figure 11
Schematic of the three-point bending beam of concrete.Figure12 presents the crack path on the three-point bend beam after the calculation using the proposed method. Most of the cracks occur along the interface layer between the aggregate and the cement mortar, which corresponds to ordinary rules [25, 26]. This example demonstrates that the proposed method has the potential of becoming a mature technology for simulating complex deformation and cracking behavior.Figure 12
Crack path of the three-point bending beam of concrete.


## 5. Conclusion

This study presents a calculation method called PEM, which can simulate the crack propagation of structures under an external load. In this method, the target object is regarded as a collection of discrete particles, and every particle interacts with its surrounding particles though an interaction force. Structures with various material characteristics can be analyzed using this method by defining the different functions of the interaction force. The calculation system achieves its final response to the external load when all particles reach the equilibrium condition of forces. A calculation program is developed based on the basic theory and is applied to the simulations of three calculation cases. The comparison of the results of the proposed method with those of FEM and the experimental observation indicates that PEM can simulate complex deformation and cracking behavior.



---

*Source: 5264319-2018-06-11.xml*