Research Article  Open Access
The Partitioned Mixed Model of Finite Element Method and Interface Stress Element Method with Arbitrary Shape of Discrete Block Element
Abstract
Based on the model of rigidspring 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 threedimensional 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 ISEMFEM 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.
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 indepth, 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], elementfree 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 bodyspring 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 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 ISEMFEM 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 where is the displacement of arbitrary point in the discrete block element. expresses the three linear displacements and three angular displacements of the element centroid, which are the primary unknowns to be solved. is the shape function matrix, whose form is The elements in 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 Figure 1, is the centroid of discrete block element and is the centroid of discrete block element , respectively. is the interface of the adjacent discrete block elements. Assume that we express stress field by piece interface stress, then the interface stress 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; indicates the array assembled by the normal and shear stress of an arbitrary point on the interface; , where 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 ; ; ; is elastic matrix, whose form is
If the neighboring elements are two different media, we have where , , is Young’s modulus, is Poisson’s ratio, and stands for different media of the two sides of the interface. and are perpendicular distances from centroids of discrete elements to interface.
If the neighbour elements are of the same material, we can get
2.3. Governing Equation
As we know, the expression of virtual work principle can be written by
The left and right terms of (7) are internalvirtual work and externalvirtual work, respectively. If ISEM is used to discrete the considered structure, the internalvirtual work of discrete block element is zero because of no strain components in discrete block elements, and the interface stress of the discrete block elements has made externalvirtual 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 is the interface of the adjacent discrete block elements. Equation (8) can be expressed by matrix where 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 and , the equation above can be written as where and are the choice matrices and is global displacement array constituted by the displacement of each discrete element centroid. Since is arbitrary, the governing equation for ISEM is where represents global stiffness matrix assembled by each interface’s stiffness , ; is general load array acting on each element centroid, which can be assembled by equivalent load array of each discrete element, .
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 array ), and ISEM has no limitation and requirement for the shape of the discrete block element, as long as we deal with 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 , 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 are the defined volume coordinates. , , , and , is the volume of the considered tetrahedron 1234, and , , , and are the volume of tetrahedron P234, P134, P124, and P123, respectively, as shown in Figure 3.
The volume integration of can also be expressed by volume coordinate. We have We can calculate (15) by applying the following formula mechanically:
For calculating surface integration in , 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 are the defined area coordinates. The surface integration of 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 , it is the same as calculating surface integration in . 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 FEMISEM Model
By introducing the transitional interface 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 where and are the shape function and the centroid displacement of discrete block element in ISEM and and are the ones of FEM, respectively. Let , , and (20) 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 on accordingly is The global stiffness matrix and load array of the mixed model are respectively, where , , and are the choice matrices of ISE, transitional interface element, and FE, respectively, and and are the equivalent loads acting on discrete element centroid in ISEM and the equivalent loads acting on nodes in FEM, respectively.
After solving by , 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 Figure 5, a concentrated load acts on the free end of a cantilever beam with Young’s modulus MPa and Poisson’s ratio . Discrete this cantilever beam by employing FEM, ISEM, and the mixed model of FEMISEM, 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 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 Kg/m^{3} and that of the liquid be Kg/m^{3}. Discrete the wedge with FEM, ISEM, and the mixed model of FEMISEM, respectively. As for the mixed model of FEMISEM, divide the wedge equally along coordinate of which we use ISE with arbitrary shape of discrete block element within , transitional interface element on the section where m and FE within , as shown in Figure 6. The mesh of ISEM is the same as that of the mixed model of FEMISEM 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 m and 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 is higher than that of 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.
Acknowledgment
This work was supported by the National Natural Sciences Foundation of China (51179064, 11372099, and 11132003).
References
 O. C. Zienkwicz and R. L. Taylor, Finite Element Method, ButterworthHeinemann, Oxford, UK, 5th edition, 2000.
 R. E. Goodman, R. L. Taylor, and T. L. Brekke, “A model for the mechanics of jointed rock,” Journal of Soil Mechanics and Foundation Division, vol. 94, pp. 637–659, 1968. View at: Google Scholar
 P. A. Cundall, “A computer model for simulating progressive large scale movements in blocky rock systems,” in Proceedings of the Rock Fracture Proceeding of the International Symposium on Rock Mechanics, Nancy, France, 1971. View at: Google Scholar
 T. Belytschko, Y. Y. Lu, and L. Gu, “Elementfree Galerkin methods,” International Journal for Numerical Methods in Engineering, vol. 37, no. 2, pp. 229–256, 1994. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Kryslc, “Meshless methods: an overview and recent developments,” Computer Methods in Applied Mechanics and Engineering, vol. 139, pp. 3–47, 1996. View at: Publisher Site  Google Scholar
 G.H. S. GenHua Shi and R. E. Goodman, “Generalization of twodimensional discontinuous deformation analysis for forward modelling,” International Journal for Numerical & Analytical Methods in Geomechanics, vol. 13, no. 4, pp. 359–380, 1989. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 G. H. Shi, “Manifold method of material analysis,” in Transactions of the 9th Army Conference on Applied Mathematics and Computing, pp. 51–76, Minneapolis, Minn, USA, 1992. View at: Google Scholar
 T. Belytschko and T. Black, “Elastic crack growth in finite elements with minimal remeshing,” International Journal For Numerical Methods in Engineering, vol. 45, pp. 601–620, 1999. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. Braun and M. Sambridge, “A numerical method for solving partial differential equations on highly irregular evolving grids,” Nature, vol. 376, no. 6542, pp. 655–660, 1995. View at: Publisher Site  Google Scholar
 E. Cueto, N. Sukumar, B. Calvo, M. A. Martínez, J. Cegoñino, and M. Doblaré, “Overview and recent advances in natural neighbour galerkin methods,” Archives of Computational Methods in Engineering, vol. 10, no. 4, 2003. View at: Google Scholar  Zentralblatt MATH
 Z. Chuhan, O. A. Pekau, J. Feng, and W. Guanglun, “Application of distinct element method in dynamic analysis of high rock slopes and blocky structures,” Soil Dynamics and Earthquake Engineering, vol. 16, no. 6, pp. 385–394, 1997. View at: Publisher Site  Google Scholar
 J.S. Lin, “A meshbased partition of unity method for discontinuity modeling,” Computer Methods in Applied Mechanics and Engineering, vol. 192, no. 1112, pp. 1515–1532, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 T. Kawai, “A new discrete model for analysis of solid mechanics problem,” Seisan Kenkyn, vol. 29, pp. 204–207, 1977. View at: Google Scholar
 J. S. Zhuo, Q. Zhang, and N. Zhao, “Interface stress element method for deformable body with discontinuous medium such as rock mass,” in Proceedings of the 8th International Society for Rock Mechanics (ISRM '95), T. Fujll, Ed., pp. 939–941, Tokyo, Japan, 1995. View at: Google Scholar
 J. S. Zhuo and Q. Zhang, Interface Stress Element Method of Discontinuous Mechanics, Science Press, Beijing, China, 2000.
 Q. Zhang and Z. Q. Wang, “Interface stress element method and its application in failure analysis of concrete gravity dam,” Science China Technological Sciences, vol. 55, no. 12, pp. 3285–3291, 2012. View at: Publisher Site  Google Scholar
 Z. L. Xu, Applied Elasticity, Higher Education Press, Wiley Eastern Limited, 1992.
Copyright
Copyright © 2013 Zhang Qing et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.