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.

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) [812], 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 [1820], 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 particles , where 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 of the particles is obtained as follows:where m is the total number of elements that contain particle ; is the volume of elements containing ; and is the weight of node to , which is determined by the geometry of the grid.

2.2. Interaction Force

The definition of the interaction force among different particles is one of the key characteristics of PEM. For particle , only its surrounding particles release an interaction force, and all the surrounding particles belong to the elements containing particle . The interaction force between particle and one of its surrounding particle can be defined as follows:where is the total number of surrounding particles to ; is the distance between and ; and and are the volumes of and , respectively.

As shown in Figure 2, the interaction force is similar to the van der Waals force in molecular dynamics; that is, when is smaller than ( is the balance distance), the interaction force takes the form of repulsion, which prevents the infinite contraction of an object and results in increased repulsion. On the contrary, when is larger than , the interaction force manifests as attraction, which maintains a solid object with a certain shape; when is relatively large, the attraction force becomes near zero. Balance distance is a special distance, in which the repulsion and attraction between two particles are the same in value and the interaction force is zero.

In PEM, the constitution function is related to the material properties and the stress states of materials. Thus, material tests are necessary when determining the exact constitution function . 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 Figure 2, when s is relatively large, the attraction force becomes near zero. This special numeric feature is not consistent with the actual computation. Hence, cutoff distance is introduced. In the actual calculation, when the distance of two particles exceeds , the interaction force is zero and a crack appears between particles and .

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 aswhere expresses the number of interaction forces that cannot transmit attraction in region and 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 particles and , as shown in Figure 3. The interaction force between and 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 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: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.

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 Figure 4, which indicates model complexity. This model describes the elastic stage and damage, as expressed as follows:where and are the line slopes in different stages, as shown in Figure 4. All these parameters can also be obtained through experiments.

Second, a plastic constitutive model should be defined, considering that the external load may cause irreversible plastic deformation in many cases. Figure 5 depicts the relationship between the interaction force and distance . The most notable feature of this model is that balance distance 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 in this plastic model, where , , and are related to the material properties and can be obtained through the experiment.

2.6. Governing Equation

The external load that worked on particle in the discrete system can be described as , and the interaction force that acts on it can be expressed aswhere is the total number of elements that contain node and is the interaction force between particle and its surrounding particle . Therefore, the balance equation of particle can be expressed as

In (8), the displacements of particles are the basic unknown variables. For the elastic constitutive model, this equation can be rewritten aswhere is the distance between particles and . 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.where is the external load vectors and 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 particle can be written aswhere is the acceleration vector of the particle and is the volume of the particle. According to Newton’s law, can be expressed as

The central difference method is used in this method, and (11) can be expressed as follows:where , , and are the displacement, velocity, and acceleration of the particle, respectively; and determines the computation speed. To guarantee the convergence of the calculation, (9) should be modified as

3. Numerical Scheme of PEM

Although PEM’s basic theory is introduced in Section 2 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 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 . 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 Step 2.

Step 7. Plot the displacement nephogram and crack growth path of every step.

The above procedure is illustrated in a flowchart in Figure 6.

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 Figure 7(a), an 8-m long and 1-m wide cantilever beam is under distributed loads (). The materials of this cantilever beam are described using the linear elastic constitutive model, and the relevant material parameters are set as and [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.

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 Figure 10(a), this disc has a straight-through crack in the center; the angle between this crack and the vertical direction is ; the crack length is (r is the radius of the disc); and the pressure acting on each end of the disc is . The material parameters of this disc are the same as those of the cantilever beam in the previous example.

Figure 10(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 Figure 11. 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]: , , and ; , , , and ; and , , and , where superscripts , , and represents the aggregate, cement mortar, and the interface layer between the aggregate and the cement mortar, respectively.

Figure 12 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.

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.

Conflicts of Interest

The grants supported by the National Natural Science Foundation-Outstanding Youth Foundation and National Key Research and Development Plan will not lead to any conflicts of interest regarding the publication of this manuscript.

Authors’ Contributions

All the authors listed have approved the manuscript.


The National Key Research and Development Plan supported this work under Grant no. 2017YFC1501103-02 and the National Natural Science Foundation-Outstanding Youth Foundation also supported this work under Grant no. 51509142. The authors appreciate the KGS for their linguistic assistance during the preparation of this manuscript.