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

Authors: Zhang Qing; Zhuo Jiashou; Xia Xiaozhou
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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 a mixed model 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.


---

## Body

## 1. 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 three-dimensional problems is higher than that of FEM, which often brings about more computational cost. Therefore, it has important academic significance and applied value to 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.


## 2. Synopsis of ISEM

### 2.1. 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(1)u=Nug,
where u=[u,v,w]T is the displacement of arbitrary point in the discrete block element. ug=[ug,vg,wg,θx,θy,θz]T 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
(2)N=[1000z-zgyg-y010zg-z0x-xg001y-ygxg-x0].
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.


### 2.2. Interface Stress Formula

As shown in Figure1, c1 is the centroid of discrete block element e1 and c2 is the centroid of discrete block element e2, respectively. AB 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 M 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]:
(3)T=DL(1)(u(2)-u(1))=DL*N*ug*,
where the upper marks 1, 2 express the identifiers of the discrete block elements on the two sides of the interface, respectively; T=[σn,τs1,τs2]T indicates the array assembled by the normal and shear stress of an arbitrary point on the interface; L*=[L(2),L(1)]T, where L(i) 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)00N(2)]; ug*=[ug(1),ug(2)]T; D is elastic matrix, whose form is
(4)D=[Dn000Ds1000Ds2].Figure 1
The interface stress element model.If the neighboring elements are two different media, we have(5)Dn=(h1λ1+2G1+h2λ2+2G2)-1,Ds1=Ds2=(h1G1+h2G2)-1,
where λi=[Eμ/(1+μ)(1-2μ)]i, Gi=[E/2(1+μ)]i, E is Young’s modulus, μ is Poisson’s ratio, and  i=1,2 stands for different media of the two sides of the interface. h1 and h2 are perpendicular distances from centroids of discrete elements to interface.If the neighbour elements are of the same material, we can get(6)D=[Dn000Ds1000Ds2].


### 2.3. Governing Equation

As we know, the expression of virtual work principle can be written by(7)∫ΩσijδeijdΩ=∫ΩfiδuidΩ=∫sσpi¯δuids.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:
(8)∑e[∫ΩefiδuidΩ+∫sσepi¯δuids+∫s0eTiδuids]=0,
where Ωe is the domain of discrete element, sσe is the stress boundary of discrete block element, and s0e is the interface of the adjacent discrete block elements. Equation (8) can be expressed by matrix
(9)-∑e∫s0eδulTTds=∑e[∫ΩeδuTfdΩ+∫sσeδuTp¯ds],
where ul is the relative displacement on the interface expressed by local coordinate. Accordingly, we have
(10)ul=-L*N*ug*.
With the equation above and interface stress formula (3), (9) becomes
(11)∑eδug(*)T∫s0eN(*)TL(*)TDL*N*ug*ds=∑eδugT[∫ΩeNTfdΩ+∫sσeNTp¯ds].
By using ug=CeU and ug*=Ce*U, the equation above can be written as
(12)δUT[∑eCe(*)T∫s0eN(*)TL(*)TDL*N*ds·CeT]U=δUT∑eCeT[∫ΩeNTfdΩ+∫sσeNTp¯ds],
where Ce and Ce* are the choice matrices and U is global displacement array constituted by the displacement of each discrete element centroid. Since δUT is arbitrary, the governing equation for ISEM is
(13)KU=Q,
where K represents global stiffness matrix assembled by each interface’s stiffness kj, kj=∫s0eN(*)TL(*)TDL*N*dS; Q is general load array acting on each element centroid, which can be assembled by equivalent load array q of each discrete element, q=∫ΩeNTfdΩ+∫sσeNTp¯dS.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.


## 2.1. 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(1)u=Nug,
where u=[u,v,w]T is the displacement of arbitrary point in the discrete block element. ug=[ug,vg,wg,θx,θy,θz]T 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
(2)N=[1000z-zgyg-y010zg-z0x-xg001y-ygxg-x0].
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.


## 2.2. Interface Stress Formula

As shown in Figure1, c1 is the centroid of discrete block element e1 and c2 is the centroid of discrete block element e2, respectively. AB 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 M 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]:
(3)T=DL(1)(u(2)-u(1))=DL*N*ug*,
where the upper marks 1, 2 express the identifiers of the discrete block elements on the two sides of the interface, respectively; T=[σn,τs1,τs2]T indicates the array assembled by the normal and shear stress of an arbitrary point on the interface; L*=[L(2),L(1)]T, where L(i) 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)00N(2)]; ug*=[ug(1),ug(2)]T; D is elastic matrix, whose form is
(4)D=[Dn000Ds1000Ds2].Figure 1
The interface stress element model.If the neighboring elements are two different media, we have(5)Dn=(h1λ1+2G1+h2λ2+2G2)-1,Ds1=Ds2=(h1G1+h2G2)-1,
where λi=[Eμ/(1+μ)(1-2μ)]i, Gi=[E/2(1+μ)]i, E is Young’s modulus, μ is Poisson’s ratio, and  i=1,2 stands for different media of the two sides of the interface. h1 and h2 are perpendicular distances from centroids of discrete elements to interface.If the neighbour elements are of the same material, we can get(6)D=[Dn000Ds1000Ds2].


## 2.3. Governing Equation

As we know, the expression of virtual work principle can be written by(7)∫ΩσijδeijdΩ=∫ΩfiδuidΩ=∫sσpi¯δuids.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:
(8)∑e[∫ΩefiδuidΩ+∫sσepi¯δuids+∫s0eTiδuids]=0,
where Ωe is the domain of discrete element, sσe is the stress boundary of discrete block element, and s0e is the interface of the adjacent discrete block elements. Equation (8) can be expressed by matrix
(9)-∑e∫s0eδulTTds=∑e[∫ΩeδuTfdΩ+∫sσeδuTp¯ds],
where ul is the relative displacement on the interface expressed by local coordinate. Accordingly, we have
(10)ul=-L*N*ug*.
With the equation above and interface stress formula (3), (9) becomes
(11)∑eδug(*)T∫s0eN(*)TL(*)TDL*N*ug*ds=∑eδugT[∫ΩeNTfdΩ+∫sσeNTp¯ds].
By using ug=CeU and ug*=Ce*U, the equation above can be written as
(12)δUT[∑eCe(*)T∫s0eN(*)TL(*)TDL*N*ds·CeT]U=δUT∑eCeT[∫ΩeNTfdΩ+∫sσeNTp¯ds],
where Ce and Ce* are the choice matrices and U is global displacement array constituted by the displacement of each discrete element centroid. Since δUT is arbitrary, the governing equation for ISEM is
(13)KU=Q,
where K represents global stiffness matrix assembled by each interface’s stiffness kj, kj=∫s0eN(*)TL(*)TDL*N*dS; Q is general load array acting on each element centroid, which can be assembled by equivalent load array q of each discrete element, q=∫ΩeNTfdΩ+∫sσeNTp¯dS.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.


## 3. 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 arrayq), 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 Figure2), we will do the following: divide it into m 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 x, y, and z of an arbitrary point in the tetrahedron can be expressed as follows:
(14)x=∑i=14Lixi,y=∑i=14Liyi,z=∑i=14Lizi,
where L1~L4 are the defined volume coordinates. L1=V1/V, L2=V2/V, L3=V3/V, and L4=V4/V, V is the volume of the considered tetrahedron 1234, and V1, V2, V3, and V4 are the volume of tetrahedron P234, P134, P124, and P123, respectively, as shown in Figure 3.Figure 2
The discrete element with arbitrary shape and its division.Figure 3
Volume coordinates of tetrahedron.The volume integration ofq can also be expressed by volume coordinate. We have
(15)∫ΩeNT(x,y,z)fdΩ=∑i=1m∫ΩieNT(L1,L2,L3,L4)fdΩ.
We can calculate (15) by applying the following formula mechanically:
(16)∫ΩL1aL2bL3cL4ddΩ=a!b!c!d!(a+b+c+d+3)!6V.For calculating surface integration inq, in the same manner, we will divide each lateral of the discrete element with arbitrary shape into n triangles. The area and centroid of which can also be solved easily. By introducing the area coordinates, the x, y, and z of an arbitrary point in the triangle can be expressed as follows:
(17)x=∑i=13Aixi,y=∑i=13Aiyi,z=∑i=13Aizi,
where A1~A3 are the defined area coordinates. The surface integration of q can also be expressed by area coordinate. We have
(18)∫sσeNT(x,y,z)p¯dS=∑i=1n∫(sσe)iNT(A1,A2,A3)p¯dS.
We can calculate (18) by applying the following formula mechanically:
(19)∫sA1aA2bA3cdS=a!b!c!(a+b+c+2)!2A,
where A is the area of the considered triangle.As for calculating integration inkj, 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.


## 4. Subregion Mixed FEM-ISEM Model

By introducing the transitional interface element, Figure4 shows the relation between ISE mesh and FE mesh. sjo is the crosssection of the two meshes. The relative displacement of an arbitrary point on sjo 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 M on sjo; for example, suppose that its relative displacement is δl, we have
(20)δl=δl(1)+δl(2)=-L(1)NISEuISE-L(2)NFEuFE=L(1)(NFEuFE-NISEuISE),
where NISE and uISE are the shape function and the centroid displacement of discrete block element in ISEM and NFE and uFE are the ones of FEM, respectively. Let N**=[NFE,-NISE], U**=[uFE,uISE]T, and (20) becomes
(21)δl=L(1)N**U**.
The interface stress of sjo is
(22)Tl=Dδl=DL(1)N**U**.
Then on the interface sjo, the potential energy accumulated by the work that is induced by the interface stress resisting the relative displacement on this interface is
(23)W0=∫Sj012δlTDδLds.
The stiffness matrix Kj0 on sjo accordingly is
(24)kj0=∫SJ0(N**)T(L(1))TDL(1)N**ds.
The global stiffness matrix and load array of the mixed model are
(25)K=∑j(Ce*)TkjCe*+∑j0(Ce**)Tkj0Ce**+∑neCeTkCe,Q=∑jCeTRISEe+∑neCeTRFEe,
respectively, where Ce*, Ce**, and Ce are the choice matrices of ISE, transitional interface element, and FE, respectively, and RISEe and RFEe are the equivalent loads acting on discrete element centroid in ISEM and the equivalent loads acting on nodes in FEM, respectively.Figure 4
Transitional interface element.After solvingU 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.


## 5. Numerical Examples

Example 1.
As shown in Figure5, a concentrated load P acts on the free end of a cantilever beam with Young’s modulus E=2.0×104 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.Table 1
The calculate results of cantilever beam.
Methods
Stress of point A (10−3 MPa)
Stress of point B (10−3 MPa)
Displacement of C (10−4 m)
Results of FEM
28.72
14.33
3.55
Results of ISEM
29.48
15.02
3.82
Results of FEM-ISEM
29.31
15.56
3.75
Theoretic solutions [17]
30.00
15.00
4.00Figure 5
Cantilever beam and the mesh.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α=tg-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×103 Kg/m3 and that of the liquid be γ=1.0×103 Kg/m3. 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 y coordinate of which we use ISE with arbitrary shape of discrete block element within y∈[0,6], transitional interface element on the section where y=6 m and FE within y∈[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.Figure 6
Wedge and the mesh.Figure 7
FEM mesh.In Tables2 and 3, we indicate the normal stress σy of some points on the sections where y=4 m and y=8 m, respectively, and make a comparison between the calculation results and the theoretic ones.Table 2
Normal stressσy of section y = 4 m (MPa).
Methods
x = 0.25 m
x = 0.75 m
x = 1.25 m
x = 1.75 m
Results of FEM
0.0369
-
0.0194
-
0.0757
-
0.1310
Results of ISEM
0.0362
-
0.0204
-
0.0766
-
0.1329
Result of FEM-ISEM
0.0359
-
0.0196
-
0.0767
-
0.1330
Theoretic solutions [17]
0.0361
-
0.0202
-
0.0764
-
0.1327Table 3
Normal stressσy of section y = 8 m (MPa).
Methods
x = 0.50 m
x = 1.50 m
x = 2.50 m
x = 3.50 m
Results of FEM
0.0786
-
0.0435
-
0.1607
-
0.2512
Results of ISEM
0.0711
-
0.0392
-
0.1542
-
0.2676
Result of FEM-ISEM
0.0701
-
0.0385
-
0.1582
-
0.2703
Theoretic solutions [17]
0.0722
-
0.0404
-
0.1528
-
0.2654From the results we can see that the accuracy of the normal stress on sectiony=4 is higher than that of y=8 m, which attributes to the fact that the accuracy of FEM is lower than that of ISEM and that of transitional interface element.


## 6. Concluding Remarks

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.



---

*Source: 950696-2013-11-07.xml*