An interactive method of interface boundary elements and partitioned finite elements for local continuous/discontinuous deformation problems

The interactive method of interface boundary elements (IBEs) and partitioned finite elements (PFEs) is proposed for solving such problems as local continuous/discontinuous deformation (i.e., landslide, concrete cracking, rock mass joints, and internal contraction joints) in concrete dams. The system is divided into continuous displacement bodies and continuous stress joints. The continuous displacement bodies are solved using PFE with the nodal displacements treated as variables, and the rigid displacements in each body and the constraining internal forces on the boundary interface are solved using IBE based on the continuous stress condition and the static force equilibrium condition in each body. Each IBE consists of all interface boundary nodes in a body, and the flexibility matrix is formed using PFE or theoretical analysis. Using this method, a nonlinear iteration procedure is carried out only on the possible discontinuous interface, an approach that greatly improves the computational efficiency. Three numerical examples are used to verify the correctness and validity of the proposed method. © 2014 The Authors. International Journal for Numerical Methods in Engineering published by John Wiley & Sons Ltd.


INTRODUCTION
Problems such as concrete cracking, unstable sliding of rock or soil slopes, and opening or closing of internal contraction joints in a concrete dam are board concern for engineering safety. Such problems can be attributed to situations involving local continuous/discontinuous deformation, and the following methods are primarily used to address them. The first method is finite element method with thin layer elements, interface elements, or spring elements reflecting the local discontinuous change in the deformation [1][2][3][4]. The elasto-plastic model is commonly used to reflect the material's nonlinear characteristics in the local discontinuous deformation region. Selected special methods are presented for contact problems [5][6][7][8][9][10]. In those methods, the stiffness matrix should be updated for the entire region. Furthermore, the status conversion of open versus closed for the discontinuous deformation interface is difficult to simulate using the elasto-plastic constitutive law.
The second method is based on the discrete element method [11], and spring elements are used to describe the relative motion between different continuums. Because most of the regions exist in an elastic state, the deformation due to continuous displacement bodies should be considered. Based on 535 this observation, the deformable discrete element method has been proposed in recent years [12][13][14]. This method expresses deformation of the continuous region in the form of mode composition and yields the deformation of the continuum by solving for the coefficient of mode composition. The disadvantages of the method are as follows: (i) a dynamic explicit solution method must be used to solve a static problem; (ii) if the continuum is composed of complex shapes that cover large areas, the mode composition cannot accurately reflect the changes in the displacement in any location, especially in the area near the stress concentration; (iii) the stiffness of the spring element is difficult to determine; and (iv) the method cannot accurately simulate transition problems from continuity to discontinuity, that is, cracks or the process of sliding surface formation.
To overcome this disadvantage, the discrete element method using a cohesive beam model is used to simulate continuous material in [15]. This approach propose a complex methodology designed to tackle the difficulties in simulating such material properties as Young's modulus, Poisson's ratio, and density, to fit the requirements of the static and dynamic mechanical behavior of the structure.
The third method is known as discontinuous deformation analysis (DDA). This method represents blocks in continuum portions as constant strain elements by considering the rigid motion of the continuum at the same time. The displacements in the continuum are chosen as the unknown variables, and the spring is used to describe the relationship between the relative displacements and forces acting on the boundaries of the continuum blocks. Static and dynamic problems can be solved using this approach. Moreover, the implicit solution method can be used for dynamic problems. The disadvantages of the method are as follows: (i) the constant strain mode is used within the continuum, which means that the size of each continuum can only take on the size of a general finite element; the actual continuum must be divided into several smaller continuums, and a reasonable stiffness must be chosen for the spring elements to ensure the original continuity; (ii) the mechanical parameters are difficult to be determined for springs, and the method cannot accurately simulate the transition from continuity to local discontinuity.
Obviously, it is inappropriate to singly adopt FEM, DEM(distinct element method), or DDA for problems involving continuity and discontinuity [16]. A combined method of using FEM and DEM has been widely used in recent years. Two types of combination, that is, internal coupling and external coupling, are commonly applied. In the external combination is that DEM and FEM methods are used for different zones and they are connected by defining a special transition layer that links the discrete zone and the finite element zone [17][18][19]. The internal combination mixes the FEM and DEM together [20,21]. These methods are dedicated to simulating the structural failure process and are complicated by searching the contact area between different discrete bodies.
For most structural problems with local discontinuities, the states prior to the instability or failure receive greater attention. A simple method known as the interactive method of interface boundary element (IBE) and partitioned finite element (PFE) is presented in this paper. The displacement continuous bodies are solved by PFE with the nodal displacements treated as variables, and the rigid displacements in each body and the constraining internal forces on boundary interface are solved by IBE based on the continuous stress condition and the static force equilibrium conditions in each body. Each IBE consisted of all interface boundary nodes in a body, and the flexibility matrix is formed via PFE or theoretical analysis. The method leverages the advantage of FEM, DDA, and DEM, and a nonlinear iteration procedure is carried out only on possible discontinuous interfaces, which greatly improves the computational efficiency.

DESCRIPTION OF LOCAL CONTINUOUS/DISCONTINUOUS DEFORMATION PROBLEMS
As shown in Figure 1, a typical wall constructed by blocks on a foundation is used to describe the local continuous/discontinuous deformation problem. Six blocks are connected by three horizontal joints and two vertical joints located on the foundation. To obtain the displacement response and stress response on the wall resulting from self-weight and the external load, it is advisable to divide the system into two components. The first component is composed of the blocks and the foundation. Assuming that the rigid displacements and the internal forces acting on the joints are known, the displacement and stress for each block and the foundation 536 T. LI, X. LIU AND L. ZHAO  can be solved using the PFE method. For the second component, the joints are primarily affected, and the rigid displacements are taken into account for the blocks and the foundation.

Partitioned finite element equations
Assuming that the entire structure is divided into n blocks by several discontinuous joints, the total displacement ¹ı t º i for each node in any block can be discomposed into the rigid displacement ¹ı g º i , including translation components and rotation components, which move with the block's centroid and displacement ¹ıº i relative to the centroid.
The number of rigid displacements for unknown variables depends on the constraint levels of the block. A total of m groups of discontinuous joints and p j .j D 1; 2; : : : ; m/ pairs of connected points on the jth discontinuous joint are assumed to simulate the discontinuous deformation and to transfer the interactive forces between different blocks. Collectively, there are mp mp D m P j D1 p j ! pairs of connected points. The local coordinate system at each pair of connected points (k D 1; 2; : : : ; mp/ can be expressed as k Á k & k with the normal direction of & k as shown in Figure 2, and the transformation matrix to the global coordinate system of xy´is rr k .
r x r xÁ r x& r y r yÁ r y& r´ r´Á r´& where r ij .i D x; y;´/.j D ; Á; &/ is the cosine of the angle between the two axes of i and j . If the thickness of the fill medium in the discontinuous joint is small, the value and the direction of the stresses crossing the discontinuous medium generally can be considered as unchanged. That is, the continuous stress interfaces can be assumed for the joints. Assuming that the interface forces ¹F cl º D ® f cl 1 f cl 2 f cl k f cl mp¯T in the connected points under the local coordinate system, and the interface forces f cl k D ® f cl k f cl kÁ f cl k&¯T are imposed on the kth pair of connected points, the interface forces under the global coordinate system are written as where f cg k and f cl k have the following relationship Once the interface forces on the discontinuous deformation interfaces are known, the displacement ¹ıº i i D 1; 2; : : : ; n relative to the centroid point in each block can be solved using the general finite element method. The PFE equations will take the form of where ¹F º i is the external load vector applied on the ith block and ¹Rº i is the load vector transferred from the interface forces where np i is the number of points in the ith block and ¹Rº i can be expressed as follows: of which OET i D 2 6 6 6 6 4 t 1;1 t 1;2 t 1;l t 1;mp t 2;1 t 2;2 t 2;l t 2;mp t k;1 t k;2 rc k;l t k;mp : : : t np;1 t np;2 rc np;l t np;mp 3 7 7 7 7 5 i where t k;l is the transformation factor of the interface forces of the lth pair of connected points to the load vector of the kth node (k D 1; 2; : : : ; np i / in the ith block. If the kth node is different from any point of the lth pair of connected points, then t k;l D 0. If the kth node is the same as the first point of the lth pair of contact points, then t k;l D rr T k . If the kth node is the same as the second point of the lth pair of contact points, then t k;l D rr T k .

Nodal rigid displacements in each block
If the displacements in the centroid of the ith block are taken as the rigid displacement variables, the rigid displacements of every point in the block can be described according to the rigid motion rules.
For the 2D problem, the centroid displacements of the ith block, which are composed of translation components and rotation components, are assumed as follows: The rigid displacements of the kth node with a coordinate value of .x; y/ relative to the centroid . N x; N y/ in the ith block will be For the 3D problem, the centroid displacements of the ith block are The rigid displacements of the kth node with a coordinate value of .x; y; ´/ relative to the centroid . N x; N y; Ń/ in the ith block will be

Static force equilibrium equation for rigid blocks
For every block, the dynamic force equilibrium equation of the discrete element is where m i and I i are the mass and rotational inertia relative to the centroid of the ith block. For the static problem, the left side of equation (6) is zero, and the equation can be rewritten as

Boundary nodal displacements due to block deformations and rigid movements
The number of variables ntotv of the relative displacements in the local discontinuous interfaces is equal to the product of mp and the dimension dim. The mp is the number of pairs of connection points. In this work, we take dim D 2, that is, The relative displacement equations at the interface boundaries will be where ® ı L¯i s the displacement vector of the second point in the pair of connection points corresponding to the negative direction in the local coordinate system, and ® ı L¯i s the displacement vector of the first point in the pair of connected points corresponding to the positive direction in the local coordinate system.
Assuming that the number of discontinuous pairs of connection points is nb i in the ith block, its relationship with the global discontinuous pairs can be built using the PFE method. For each block i, the total displacements of each node in the block can be calculated by equations (3), (5), and (7).
Substituting formula (4) into equation (9) and multiplying equation (9) by OET T i i on both sides, we obtain in which under the corresponding local coordinate system in the local discontinuous interfaces due to the ith block, ® N ı¯b i D OET T i OEK 1 ¹F º i is the contributed displacement vector caused by external force imposed on the ith block under the corresponding local coordinate system in the local discontinuous interfaces, and OET T i OE! i ¹ º i is the rigid movement vector on the ith block under the corresponding local coordinate system in the local discontinuous interfaces.
Imitating the finite element method, each item in equation (10) is considered in terms of each item in element i (the node number of the element is nb i /. After accumulating and assembling all of the corresponding items in equation (8) for all elements, the following equation can be obtained.
OEr j

Interface boundary equations
Assuming that the relative displacements change with the interface internal forces, equation (11) can be rewritten as where˛is the flexibility matrix used to reflect the relative displacements and internal forces on the boundary interface. That is, Applying equation (7) to all blocks and considering equation (4), we obtain Integrating equations (12) and (13) results in "

Classical Hertz contact problem
The first example is the solution of the classical Hertz contact problem. Figure 3(a) shows the contact problem of an infinite cylinder on a half space acted upon by the vertical uniform force p on the neutral axis of the cylinder. The analytical solutions are [22] The half contact width , are elastic modulus and Poisson's ratio for cylinder and E 2 ; 2 are elastic modulus and Poisson's ratio for half space, respectively.
In the example, the following parameters are used: E 1 D E 2 D 1:0 10 5 MPa, 1 D 2 D 0:3, r D 100 mm, and p D 2:0 10 4 N/mm. The analytic solution is as follows.
The half contact width b D 6:808 mm, The maximum pressure q max D 1870:277 N/mm 2 The interactive method of IBE and PFE is used to obtain the numerical solution. The 600 mm 300 mm space is simulated for the foundation, and the normal constraints are imposed on the side and bottom boundaries. The mesh is shown in Figure 3(b). The local refinement technique shown in Figure 3(c) is used to accurately simulate the deformation and contact pressure distribution near the contact surface. A total of 77 pairs of nodes are selected on the possible contact interface. The cylinder and foundation are simulated using PFE, and the possible contact interface is simulated by IBE.
The steps for solving this example are described as follows: Step 1: Choose the nodal displacements to represent the rigid movements of PFE. The foundation has no rigid movement, and no additional constraint is required except for the given boundary constraints. The rigid displacements including . N u N v Â / should be constrained with a nonzero value. For PFE, rotation is not variable at any node and can be represented by the vertical displacements of two different nodes. Therefore, the constraint rigid displacement can be substituted by the constraint nodal displacements. In this example, the optional points of A and B in the cylinder are shown in Figure 3(b). The horizontal and vertical displacements of node A and the vertical displacement of node B are constrained and will be solved with the rigid movements of the cylinder.
Step 2: Create 77 pairs of the connection nodes in the possible discontinuous interfaces by finding the local direction cosine to form the transformation matrix rr k for shifting from the local coordinate to the global coordinate system and controlling the area s k , the initial gap of gap0 k , and the initial contact state of state0 k D state k , in which 0 represents opening and 1 represents contacting.
Step 3: For each PFE, form the collection of the boundary nodes nb i that belong to any pair of the connection points in the possible discontinuous interfaces and indicate the position of the first or the second point in the connection pair. Assemble equation (3) to form the displacement vector ® N ı¯b i caused by external forces, the flexibility matrix OEC i in equation (10), and the transformation matrix OET T i OE! i , which transforms the rigid movement of the block to the boundary nodes.
Step 4: Determine the number of unknown contact forces. This value is the sum of the unknown contact forces for pairs with state k D 1. Form equation (14) with˛D 0, and solve it to obtain the contact forces and the rigid movement. The displacements for any node can be solved with equation (9).
Step 5: Update the new contact state. First, the gap for each pair of connection points is determined.
where w k is the total displacement along the normal direction at the second or first point of the kth pair.
The new contact state can be judged as follows: If gap k <D 0, state k D 1, and if gap k > 0, state k D 0.
Step 6: If state k D state0 k for all pairs of connection points, the contact iteration has converged, and the stresses can be obtained with PFE. Otherwise, take state0 k D state k , and repeat Step 4 to Step 6.
Step 7: Repeat Step 4 to Step 6 for addition of all load increments.
The main results are the half contact width b D 6:74 mm and the maximum pressure q max D 1904:80 N/mm 2 .
It is shown that the results are in good agreement with the analytic solutions. Figure 3(d) shows the contact area if the cylinder is subjected to the given load. Figure 3(e) and 3(f) show the stress distributions of the horizontal and vertical stresses near the contacted area. Figure 3(g) shows the pressure distribution at the contacted interface.

The 3D Hertzian contact problem
This example simulates a 3D elastic sphere on a half space. A concentrated force P is applied to the top surface of the sphere, and the distribution of the contact force is analyzed. Compared with the theoretical solution, the calculated solution is proven to be reasonable, which verifies the correctness and validity of the proposed method in the 3D case.
The analytical solutions are [23] Contact half width a D 3P 0 R  The 300 mm 300 mm 120 mm domain is simulated, and the normal constraints are imposed on the side and bottom boundaries. The finite element model is shown in Figure 4(a), and the local refinement mesh of the contact region is shown in Figure 4(b). A total of 3001 pairs of nodes are selected on the possible contact interface. The vertical displacements of node C,D and E, horizontal displacements along x direction of node D and E, and horizontal displacement along y direction of node D are assumed to be solved via the rigid movements of the sphere. The steps for solving this example are the same as in the first example. Following the analytical solutions, the results of this problem are as follows: the half contact width a D 6:00 mm, and the maximum pressure q max D 958:12 N=mm 2 .
Obviously, the results are in good agreement with the analytical solutions. The pressure distribution in the contacted interface is shown in Figure 5, and the maximum pressure q max under a different concentrated force P is presented in Figure 6. With the increase of P , q max exerts a tendency to increase at the same time, reflecting the nonlinear behavior. Similarly, the contact width a under a different concentrated force P is shown in Figure 7; the half contact width a gradually increases with the increase of P . The values of q max and a agree well with those of theoretical solution, which indicates that the method proposed in this paper is accurate. In addition, the distributions of ZZ in the contact area under the concentrated forces 0.1P 0 , 0.4P 0 , and P 0 are illustrated in Figure 8

Cracking process of the notched three-point bending beam
A three-point bending test with a notch depth equal to the half height of the beam is taken as an example. This test was made by Petersson [24] in 1981. The interactive method of IBE and PFE presented in this paper is used to simulate the fracture process, and the fictitious crack model (FCM) presented by Hillerborg [25] is also used to simulate the fracture process.
For the bilinear curve model of Petersson [24], the contact normal stress with relative normal displacement or relative normal displacement with a normal contact force is divided into two components.
For the Cornelissen curve model, the recommended values of c 1 D 1:0 and c 2 D 5:64 are used, and the limit crack width of w 0 D 209 m is computed from the fracture energy G f . The contact normal stress with relative normal displacement is The incremental relative normal displacement with incremental normal contact force is where B 0 .w/ is the derivative of B (w/ with respect to w. The finite element mesh is shown in Figure 9(a), in which the segment of BC represents the given notch, and 51 pairs of contact points are arranged in segment AB.
With the interactive method of IBE and PFE, the specimen is divided into two components, as shown in Figure 9(b) and 9(c). These components are connected via the possible discontinuous interface AB. To obtain the relationship between the load P and the vertical displacement at point A, the specified displacement method is used at point A.
The steps for solving this example are shown as follows.
Step 1: Choose the nodal displacements to represent the rigid movements of PFE. Because of the constraints in both the horizontal and vertical directions at point D and the specified displacement at point A, there is no rigid displacement in the left component. Additionally, only the vertical displacement is constrained at point E, and the rigid movement in the horizontal direction should be solved for the right component. The horizontal displacement at point F, which could be arbitrarily chosen for the right component, is assumed to be obtained through the rigid movement.
Step 2: Create 51 pairs of connection nodes in the possible cracking interface to find the local direction cosine to form the transformation matrix rr k for shifting from the local coordinates to the global coordinates while controlling area s k and the initial contact state of state0 k D state k . When the connecting stress is smaller than the tensile strength, displacements between connected points are completely consistent with state k D 1. If the connecting stress is larger than the tensile strength, the displacements between connected points are determined by the FCM with state k D 3.
Step 3: Same as Step 3 in Example 1. Step 4: Determine the number of unknown contact forces, which is the sum of the unknown contact forces for the pairs. If state k D 1, all components of the contact forces are unknown. When state k D 3, only the contact force along the normal direction is unknown creating. Forming equation (14), in which˛of flexibility matrix to reflect the relative displacement and internal force on the boundary interface is varied with the contact state. The contact forces and the rigid movement are obtained by solving the equation. The displacements for any node can be solved with equation ( Step 5: Update the new contact state according to the normal contact force. If f cl k& <D f t s k , state k D 1, and if f cl k& > f t s k , state k D 3 Step 6: If state k D state0 k for all pairs of connection points, the contact iteration has converged, and the stresses can be obtained by PFE. Otherwise, take state0 k D state k , and repeat Step 4 to Step 7.
Step 7: Repeat Step 4 to Step 7 for additional load increments.
A total of 50 steps with a vertical displacement increment of 0.02 mm at point A are used to simulate the crack process of the specimen via the interactive method of IBE and PFE.
The softening curve of the concrete is shown in Figure 10. Figure 11 gives the comparison of the predicted relationship between the load and the displacement at the action point with that of the test. The result with the softening of the Cornelissen curve or bilinear curve is between the upper and lower envelope lines and is much better than that of the linear curve. Figure 12 shows the relationship between the load and the crack opening displacement at the crack endpoint, and Figure 13 shows the relationship between the load and the rigid horizontal displacement of the right block of the beam. The maximum loads in Figure 12 are the same as those in Figure 13. Figure 14 shows the deformation pattern at the end of the load with the softening of the Cornelissen curve. Figures 15 and 16 show the horizontal and vertical stress distributions around the crack tip, respectively.

Mesh sensitivity analysis.
The mesh size is a factor of great importance in the simulation of crack propagation. To investigate the mesh sensitivity of the results and to verify to what extent the proposed method is dependent on the mesh size, an example of a linear softening model is given to obtain the load-displacement curve with a mesh size ranging from 0.2 to 20 mm, which is approximately 1/50 to 1/5 of the seam height, as shown in Figure 17. As observed from the figure, the results obtained with different mesh sizes are nearly identical if the size is less than 1/10 of the seam height. If the mesh size is greater than 1/5 of the seam height, the load-displacement curve is no longer smooth because the size is set too large. However, in this case, the results still agree well with those obtained from the case in which the mesh size is small, and large errors only exist at the end of the curve. In summary, it is obvious that the method proposed in this paper has only small dependence on the mesh size.

Safety factor for determined sliding surface
A typical slope stability problem is shown in Figure 18, and the main coordinates are marked in the figure. The example in [26] is used to assess the sliding modes and the safety factor with three types       of soils. The results from the rigid limit equilibrium method are between 1.23 in. and 1.52 in. [15]. In this work, the rigid limit equilibrium method is used to assess the method presented in this paper. We assume that the sliding surface is known as an arc with a center coordinate of (19.61, 46.32) and a radius of 21.32, the elastic modulus of the sliding block is 1.0 10 4 kPa, the density is 20 kN/m 3 , and the values of c and tan ' on the sliding surface are assumed to be 3.0 kPa and 19.6 ı , respectively.
The problem is divided into two blocks for the sliding body and foundation as shown in Figures 19  and 20. The full finite mesh is shown in Figure 21, and the finite element mesh for the sliding body is enlarged in Figure 22. A series of coefficients is used to reduce c and tan ' to obtain the safety factor when the sliding surface is in the state of limit equilibrium.
The steps for solving this example are described as follows.     Step 1: Choose the nodal displacements to represent the rigid movements of PFE. The horizontal and vertical displacements of node A and the vertical displacement of node B, which are labeled in Figure 15, are assumed to be solved via the rigid movements of the sliding body. The normal constraint conditions are used for the latter and bottom boundaries.
Step 2: Create 75 pairs of connection nodes on the sliding surface, and calculate its local direction cosine to form transformation matrix rr k for shifting from local coordinates to global coordinates by controlling the area s k , the initial gap of gap0 k , and the initial contact state of state0 k D state k , in which 0 represents opening, 1 represents contacting, and 2 represents sliding.
Step 3: Same as Step 3 in Example 1.
Step 4: Determine the number of unknown contact forces, which is the sum of unknown contact forces for pairs with state k D 1. Construct equation (14) with˛D 0, and solve it to obtain the contact forces and the rigid movement. The displacements for any node can be solved with equation (9).
Step 5: Determine the number of unknown contact forces, which is the sum of the unknown contact forces for a pair with state k > 0. If state k D 1, all components of the contact forces are unknown. If state k D 2, the Mohr-Coulomb yield criterion is obeyed between the shear contact force and the normal contact force. Construct equation (14), in whichį s zero for state k D 1 and a large value in the direction of sliding for state k D 2. The contact forces and rigid movement can be solved using equation (14). The displacements for any node can be solved with equation (9).
Step 6: Update the new contact state according to the contact forces. If Step 7: If state k D state0 k for all pairs of connection points, the contact iteration has converged, and the stresses can be obtained by PFE. Otherwise, take state0 k D state k , and repeat Step 4 to Step 6.
Step 8: Repeat Step 4 to Step 6 for additional load increments.
The finite element mesh for the sliding body is shown in Figure 23, and Figure 24 gives the relationship between the strength reduction coefficient and the horizontal, vertical, and rotational rigid movements for the sliding body. If the strength reduction coefficient is set to 1.37, all rigid movements are suddenly increased, and all contact pairs of nodes are in the sliding state. Therefore, the safety factor can be taken as K D 1:30. Figure 25 gives the deformation mode and the displacement vector for the slope of the critical sliding state.

CONCLUSION
This work presents the interactive method of IBE and PFE for solving problems that include local continuous/discontinuous deformation. The continuous displacement bodies are solved by PFE with the nodal displacements treated as variables, whereas the rigid displacements in each body and the constraining internal forces on the boundary interface are solved by IBE based on the continuous stress condition and the static force equilibrium condition in each body. Each IBE consists of all interface boundary nodes in a body, and the flexibility matrix is constructed with PFE or theoretical analysis. With this method, a nonlinear iteration procedure is carried out only on the possible discontinuous interface, which largely improves the computational efficiency.
The method can be used to predict the structural failure process before crushing, that is, the stability of a landslide, crack propagation of a concrete structure, and opening and closing of rock mass joints or internal contraction joints in concrete dams. This method also can be used to simulate the frictional contact problem with or without initial gaps. The main differences in the algorithm are the number of unknown contact forces and the use of the flexibility matrix to reflect the relative displacement and internal force on the boundary interface. No spring is required for the possible discontinuous surface if it is in a continuous state.
Three numerical examples are used to verify the correctness and validity of the proposed method. The first example is the solution of the classical Hertz contact problem. The equations should be updated in the full domain if using general FEM, and the dynamic equation should be solved if using general DEM. The interactive method of IBE and PFE will simplify the iteration process because the nonlinear updating is only limited to the possible contact interface. The resulting half contact width and the maximum pressure are quite close to the analytical solution. The second example is fracture process simulation for a three-point bending test with a notch depth equal to the half height of the beam. This model can be easily solved using the interactive method of IBE and PFE. The FCM is used on the cracking zone, and the linear, bilinear, and Cornelissen curve models are used to reflect the relationship between the normal contact force and the opening width. It is shown that the result from the softening of the Cornelissen curve or bilinear curve between the upper and lower envelope lines is much better than that of the linear curve. The third example is used to simulate the sliding modes and to find the safety factor for a typical slope stability problem. Using the method proposed in this paper, a series of coefficients is used to reduce c and tan ' to obtain the safety factor when the sliding surface is in the state of limit equilibrium. The result is similar to that of the reference, and the deformation mode or the displacement vector for the slope in the critical sliding state is reasonable.
Based on the numerical examples, it is concluded that the proposed method is suitable for the modeling of local continuous/discontinuous deformation problems.