The Partitioned Mixed Model of Finite Element Method and Interface Stress Element Method with Arbitrary Shape of Discrete Block Element

Based on the model of rigid-spring element suitable for homogeneous elastic problem, which was developed by Japanese professor Kawai, the interface stress element model (ISEM) for solving the problem of discontinuous media mechanics has been established. Compared with the traditional finite element method (FEM), the ISEM is more accurate and applicable. But the total number of freedom degree of ISEM in dealing with three-dimensional problems is higher than that of FEM, which often brings about the reduction on efficiency of calculation.Therefore, it is necessary to establish amixedmodel by gathering the advantages of ISEM and FEM together. By making use of the good compatibility of ISEM and introducing the concept of transitional interface element, this paper combines ISEM and FEM and proposes a mixed model of ISEM-FEM which can solve, to a large extent, the contradictions between accuracy and efficiency of calculation. In addition, using natural coordinate, algorithm of ISEM for block elements of arbitrary shape has been performed. Numerical examples show that the method proposed in this paper is feasible and its accuracy is satisfactory.


Introduction
In recent decades, as the most efficient numerical method of computational mechanics, finite element method (FEM) has been widely used in solving engineering problems [1].With the increase of calculation complexity and application of in-depth, FEM also has gradually exposed some defects.For widespread displacement discontinuity in the practical problems, FEM is difficult to solve because of its compatible displacement model.Although FEM can simulate discontinuous deformations by using contact elements [2], it usually relates to the determination of the calculated parameters.Besides, too many contact elements will make the system of equations illconditioned to a certain degree.As a consequence, the stability of numerical results will not be ensured.
Recently, with the wide use and improvement of continuous mechanics numerical method, the discontinuous mechanics and the analysis methods accordingly develop quickly.Numerical analysis methods based on discontinuous mechanics have been proposed, such as discrete element method [3], element-free method [4,5], discontinuous deformation analysis method (DDA) [6], manifold method [7], extended finite element method [8], natural element method [9,10], and so on [11,12], and they are becoming one of the cutting edge research fields in computational mechanics.
On the base of the rigid body-spring element model that solves the homogeneity elastic problems by professor Kawai [13], we found the interface stress element model (ISEM) of the discontinuous mechanics [14].We have made certain achievements on the static and kinetic problems, field problems, and the application in engineering [15,16].The study proves that the ISEM is more accurate and applicable compared with the traditional FEM.But generally, the total number of freedom degree of ISEM in dealing with threedimensional problems is higher than that of FEM, which often brings about more computational cost.Therefore, it has important academic significance and applied value to Mathematical Problems in Engineering establish a mixed discrete model by absorbing the advantages of ISEM and FEM.
Using natural coordinate, the paper proposes the ISEM with arbitrary shape of discrete block element, then brings in transitional interface element, connects the ISEM and FEM, and establishes the ISEM-FEM mixed model with arbitrary shape of discrete block element, which cannot only enhance the simulation of the analysis system but also solve the contradiction between accuracy and efficiency of the calculation.

The Displacement Pattern of Discrete Block Element.
Assume that the discrete block element only has rigid displacement.We can construct the structural displacementfield with piecewise rigid displacement pattern and have where u = [, V, ]  is the displacement of arbitrary point in the discrete block element.u  = [  , V  ,   ,   ,   ,   ]  expresses the three linear displacements and three angular displacements of the element centroid, which are the primary unknowns to be solved.N is the shape function matrix, whose form is The elements in N should make the displacement components obtained by the previous formula (1) satisfy the condition that the strain components in discrete elements are zero.

Interface Stress Formula.
As shown in Figure 1,  1 is the centroid of discrete block element  1 and  2 is the centroid of discrete block element  2 , respectively. is the interface of the adjacent discrete block elements.Assume that we express stress field by piece interface stress, then the interface stress T of an arbitrary point  on the interface can be obtained by the relative deformation of the adjacent discrete block elements and the material properties, which is expressed as follows [15]: where the upper marks 1, 2 express the identifiers of the discrete block elements on the two sides of the interface, respectively; T = [  ,  1 ,  2 ]  indicates the array assembled by the normal and shear stress of an arbitrary point on the interface; L * = [L (2) , L (1) ]

𝑇
, where L () is a directional cosine matrix of interface local coordinate vector in the global coordinate system.Note that the local coordinate directions on interface of two adjacent elements are opposite, so we have L (2) = −L (1) ; N * = [ N (1) 0 0 N (2) ]; If the neighboring elements are two different media, we have where Young's modulus,  is Poisson's ratio, and  = 1, 2 stands for different media of the two sides of the interface.ℎ 1 and ℎ 2 are perpendicular distances from centroids of discrete elements to interface.
If the neighbour elements are of the same material, we can get

Governing Equation.
As we know, the expression of virtual work principle can be written by The left and right terms of ( 7) are internal-virtual work and external-virtual work, respectively.If ISEM is used to discrete the considered structure, the internal-virtual work of discrete block element is zero because of no strain components in discrete block elements, and the interface stress T of the discrete block elements has made external-virtual work on the relative displacement.In this case, (9) can be written as follows: where Ω  is the domain of discrete element,    is the stress boundary of discrete block element, and   0 is the interface of the adjacent discrete block elements.Equation ( 8) can be expressed by matrix where u  is the relative displacement on the interface expressed by local coordinate.Accordingly, we have With the equation above and interface stress formula (3), ( 9) becomes By using u  = C  U and u *  = C *  U, the equation above can be written as where C  and C *  are the choice matrices and U is global displacement array constituted by the displacement of each discrete element centroid.Since U  is arbitrary, the governing equation for ISEM is where K represents global stiffness matrix assembled by each interface's stiffness k  , k  = ∫   0 N ( * ) L ( * ) DL * N * ; Q is general load array acting on each element centroid, which can be assembled by equivalent load array q of each discrete element, q = ∫ Ω  N  fΩ + ∫    N  p.Equation ( 13) stands for an equilibrium equation at the centroid of each discrete element and can also be established by incomplete generalized variation theory that loosens the displacement compatibility condition on interface or weak solution form of partial differential equation.
From the process of the construction of the ISEM model, we can see that since we adopt piecewise rigid displacement pattern, the displacement can be discontinuous on interface of the discrete element.This can well express some deformation properties, such as slide, fraction, and crack.Besides, by setting plastic cells, viscous cells, contact cells, and so forth on the interface, we can solve different nonlinear problems conveniently.What is more, since the interface stress has algebraic relation with the relative deformation of adjacent elements, the precision of stresses generally is not lower than that of displacements, which enhances the reliability of stress state's criterion, and makes sure that the nonlinear solution will not drift.

The ISEM with Arbitrary Shape of Discrete Block Element
In the interface stress element method, the stiffness matrix only has relations with interface, not with the shape of element (except the load array q), and ISEM has no limitation and requirement for the shape of the discrete block element, as long as we deal with q appropriately.As to the discrete block element with arbitrary shape in space (see Figure 2), we will do the following: divide it into  tetrahedra, the volume, and centroid of which can be solved easily.When calculating the volume, integration of q, we can convert it to the sum of that of each tetrahedron.By introducing the volume coordinates, the , , and  of an arbitrary point in the tetrahedron can be expressed as follows: where  1 ∼  4 are the defined volume coordinates. 1 =  1 /,  2 =  2 /,  3 =  3 /, and  4 =  4 /,  is the volume of the considered tetrahedron 1234, and  1 ,  2 ,  3 , and  4 are the volume of tetrahedron P234, P134, P124, and P123, respectively, as shown in Figure 3.
The volume integration of q can also be expressed by volume coordinate.We have We can calculate (15) by applying the following formula mechanically: For calculating surface integration in q, in the same manner, we will divide each lateral of the discrete element with arbitrary shape into  triangles.The area and centroid of which can also be solved easily.By introducing the area coordinates, the , , and  of an arbitrary point in the triangle can be expressed as follows: where  1 ∼  3 are the defined area coordinates.The surface integration of q can also be expressed by area coordinate.We have We can calculate (18) by applying the following formula mechanically: where  is the area of the considered triangle.As for calculating integration in k  , it is the same as calculating surface integration in q.In this way, we can do integration directly and get results accordingly.
The ISEM with arbitrary shape discrete block element makes the geometric simulation of complicated structure much more convenient.

Subregion Mixed FEM-ISEM Model
By introducing the transitional element, Figure 4 shows the relation between ISE mesh and FE mesh.   is the crosssection of the two meshes.The relative displacement of an arbitrary point on    along the local coordinate is the sum of the displacement caused by the left element of ISE and the deformation induced by the right element of FE.Take an arbitrary point  on    ; for example, suppose that its relative displacement is   , we have   =  (1)   +  (2)   = −L (1)  where N ISE and u ISE are the shape function and the centroid displacement of discrete block element in ISEM and N FE and u FE are the ones of FEM, respectively.Let becomes The interface stress of    is Then on the interface    , the potential energy accumulated by the work that is induced by the interface stress resisting the relative displacement on this interface is The stiffness matrix  0  on    accordingly is The global stiffness matrix and load array of the mixed model are respectively, where C *  , C * *  , and C  are the choice matrices of ISE, transitional interface element, and FE, respectively, and R  ISE and R  FE are the equivalent loads acting on discrete element centroid in ISEM and the equivalent loads acting on nodes in FEM, respectively.
After solving U by KU = Q, we can get the displacement of discrete element centroid and that of FEM nodes with the choose matrix above and obtain the stress of each ISE, transitional interface element, and FEM by applying the stress formula of them, respectively.

Numerical Examples
Example 1.As shown in Figure 5, a concentrated load  acts on the free end of a cantilever beam with Young's modulus  = 2.0 × 10 4 MPa and Poisson's ratio  = 0.20.Discrete this cantilever beam by employing FEM, ISEM, and the mixed model of FEM-ISEM, respectively (the section which point B belongs to is transitional interface element, the left ten rows of it are FE, and the right ones are ISE).Table 1 shows the results of horizontal normal stress of A and B, the vertical displacement of C, and the comparison of them and theoretic solutions.
From the results we can see, in this problem, that the accuracy of ISEM is higher than that of FEM and the accuracy of transitional interface element is lower than ISEM  s but higher than that of FEM.
Example 2. Consider a wedge with the left edge vertical and the right one oriented at an angle  =  −1 (1/2) from the vertical edge (Figure 6) subjected to the action of gravity and the pressure of impounded liquid.The height of liquid is the same as that of the wedge.Let the density of the wedge be  = 2.4 × 10 3 Kg/m 3 and that of the liquid be  = 1.0 × 10 3 Kg/m 3 .Discrete the wedge with FEM, ISEM, and the mixed model of FEM-ISEM, respectively.As for the mixed model of FEM-ISEM, divide the wedge equally along  coordinate of which we use ISE with arbitrary shape of discrete block element within  ∈ [0, 6], transitional interface element on the section where  = 6 m and FE within  ∈ [6,8], as shown in Figure 6.The mesh of ISEM is the same as that of the mixed model of FEM-ISEM of which we use ISE in all computational domain.The mesh of FEM is shown in Figure 7.
In Tables 2 and 3, we indicate the normal stress   of some points on the sections where  = 4 m and  = 8 m, respectively, and make a comparison between the calculation results and the theoretic ones.
From the results we can see that the accuracy of the normal stress on section  = 4 is higher than that of  = 8 m, which attributes to the fact that the accuracy of FEM is lower than that of ISEM and that of transitional interface element.

Figure 2 :
Figure 2: The discrete element with arbitrary shape and its division.

Table 1 :
The calculate results of cantilever beam.

Table 3 :
Normal stress   of section  = 8 m (MPa).By making use of the compatibility of ISEM and introducing the concept of transitional interface element, this paper combines the counting methods of ISEM and FEM and proposed a mixed model of ISEM and FEM with arbitrary shape of discrete element, which can solve in a large degree the contractions between accuracy and efficiency of calculation.The examples prove the applicability of the mixed model and indicate the bright future of its application in engineering.