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

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) 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 𝐍=1000𝑧−𝑧_𝑔𝑦_𝑔−𝑦010𝑧_𝑔−𝑧0𝑥−𝑥_𝑔001𝑦−𝑦_𝑔𝑥_𝑔−𝑥0. (2) 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]: 𝐓=𝐃𝐋⁽¹⁾𝐮⁽²⁾−𝐮⁽¹⁾=𝐃𝐋^∗𝐍^∗𝐮^∗_𝑔, (3) where the upper marks 1, 2 express the identifiers of the discrete block elements on the two sides of the interface, respectively; 𝐓=[𝜎_𝑛,𝜏_𝑠1,𝜏_𝑠2]^𝑇 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 𝐋⁽²⁾=−𝐋⁽¹⁾; 𝐍^∗=^𝐍(1)𝟎_𝟎𝐍(2); 𝐮^∗_𝑔=[𝐮_𝑔⁽¹⁾,𝐮_𝑔⁽²⁾]^𝑇; 𝐃 is elastic matrix, whose form is 𝐷𝐃=^𝑛000𝐷_𝑠1000𝐷_𝑠2. (4)

Figure 1 

The interface stress element model.

If the neighboring elements are two different media, we have 𝐷_𝑛=ℎ₁𝜆₁+2𝐺₁+ℎ₂𝜆₂+2𝐺₂⁻¹,𝐷_𝑠1=𝐷_𝑠2=ℎ₁𝐺₁+ℎ₂𝐺₂⁻¹, (5) where 𝜆_𝑖=[𝐸𝜇/(1+𝜇)(1−2𝜇)]_𝑖, 𝐺_𝑖=[𝐸/2(1+𝜇)]_𝑖, 𝐸 is Young’s modulus, 𝜇 is Poisson’s ratio, and  𝑖=1,2 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 𝐷𝐃=^𝑛000𝐷_𝑠1000𝐷_𝑠2. (6)

2.3. Governing Equation

As we know, the expression of virtual work principle can be written by _Ω𝜎_𝑖𝑗𝛿𝑒_𝑖𝑗𝑑Ω=_Ω𝑓_𝑖𝛿𝑢_𝑖𝑑Ω=_𝑠𝜎𝑝_𝑖𝛿𝑢_𝑖𝑑𝑠. (7)

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 𝐓 of the discrete block elements has made external-virtual work on the relative displacement. In this case, (9) can be written as follows: _𝑒Ω𝑒𝑓_𝑖𝛿𝑢_𝑖𝑑Ω+_𝑠𝑒𝜎𝑝_𝑖𝛿𝑢_𝑖𝑑𝑠+_𝑠𝑒0𝑇_𝑖𝛿𝑢_𝑖𝑑𝑠=0, (8) 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 −_{𝑒𝑠𝑒0}𝛿𝐮^𝑇_𝑙𝐓𝑑𝑠=_𝑒Ω𝑒𝛿𝐮^𝑇𝐟𝑑Ω+_𝑠𝑒𝜎𝛿𝐮^𝑇𝐩𝑑𝑠, (9) where 𝐮_𝑙 is the relative displacement on the interface expressed by local coordinate. Accordingly, we have 𝐮_𝑙=−𝐋^∗𝐍^∗𝐮^∗_𝑔. (10) With the equation above and interface stress formula (3), (9) becomes _𝑒𝛿𝐮_𝑔^(∗)𝑇_𝑠𝑒0𝐍^(∗)𝑇𝐋^(∗)𝑇𝐃𝐋^∗𝐍^∗𝐮^∗_𝑔=𝑑𝑠_𝑒𝛿𝐮^𝑇_𝑔Ω𝑒𝐍^𝑇𝐟𝑑Ω+_𝑠𝑒𝜎𝐍^𝑇𝐩𝑑𝑠. (11) By using 𝐮_𝑔=𝐂_𝑒𝐔 and 𝐮^∗_𝑔=𝐂^∗_𝑒𝐔, the equation above can be written as 𝛿𝐔^𝑇_𝑒𝐂_𝑒^(∗)𝑇_𝑠𝑒0𝐍^(∗)𝑇𝐋^(∗)𝑇𝐃𝐋^∗𝐍^∗𝑑𝑠⋅𝐂^𝑇_𝑒𝐔=𝛿𝐔^𝑇_𝑒𝐂^𝑇_𝑒Ω𝑒𝐍^𝑇𝐟𝑑Ω+_𝑠𝑒𝜎𝐍^𝑇𝐩𝑑𝑠, (12) 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 𝐊𝐔=𝐐, (13) where 𝐊 represents global stiffness matrix assembled by each interface’s stiffness 𝐤_𝑗, 𝐤_𝑗=∫_𝑠𝑒0𝐍^(∗)𝑇𝐋^(∗)𝑇𝐃𝐋^∗𝐍^∗𝑑𝑆; 𝐐 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: 𝑥=⁴_𝑖=1𝐿_𝑖𝑥_𝑖,𝑦=⁴_𝑖=1𝐿_𝑖𝑦_𝑖,𝑧=⁴_𝑖=1𝐿_𝑖𝑧_𝑖, (14) 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.

Figure 2 

The discrete element with arbitrary shape and its division.

Figure 3 

Volume coordinates of tetrahedron.

The volume integration of 𝐪 can also be expressed by volume coordinate. We have _Ω𝑒𝐍^𝑇𝑥,𝑦,𝑧𝐟𝑑Ω=^𝑚_{𝑖=1Ω𝑒𝑖}𝐍^𝑇𝐿₁,𝐿₂,𝐿₃,𝐿₄𝐟𝑑Ω. (15) We can calculate (15) by applying the following formula mechanically: _Ω𝐿^𝑎₁𝐿^𝑏₂𝐿^𝑐₃𝐿^𝑑₄𝑑Ω=𝑎!𝑏!𝑐!𝑑!(𝑎+𝑏+𝑐+𝑑+3)!6𝑉. (16)

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: 𝑥=³_𝑖=1𝐴_𝑖𝑥_𝑖,𝑦=³_𝑖=1𝐴_𝑖𝑦_𝑖,𝑧=³_𝑖=1𝐴_𝑖𝑧_𝑖, (17) where 𝐴₁∼𝐴₃ are the defined area coordinates. The surface integration of 𝐪 can also be expressed by area coordinate. We have _𝑠𝑒𝜎𝐍^𝑇𝑥,𝑦,𝑧𝐩𝑑𝑆=^𝑛_{𝑖=1(𝑠𝑒𝜎)𝑖}𝐍^𝑇𝐴₁,𝐴₂,𝐴₃𝐩𝑑𝑆. (18) We can calculate (18) by applying the following formula mechanically: _𝑠𝐴^𝑎₁𝐴^𝑏₂𝐴^𝑐₃𝑑𝑆=𝑎!𝑏!𝑐!(𝑎+𝑏+𝑐+2)!2𝐴, (19) 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 FEM-ISEM 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 𝜹_𝑙=𝜹_𝑙⁽¹⁾+𝜹_𝑙⁽²⁾=−𝐋⁽¹⁾𝐍_ISE𝐮_ISE−𝐋⁽²⁾𝐍_FE𝐮_FE=𝐋⁽¹⁾𝐍_FE𝐮_FE−𝐍_ISE𝐮_ISE, (20) where 𝐍_ISE and 𝐮_ISE are the shape function and the centroid displacement of discrete block element in ISEM and 𝐍_FE and 𝐮_FE are the ones of FEM, respectively. Let 𝐍^∗∗=[𝐍_FE,−𝐍_ISE], 𝐔^∗∗=[𝐮_FE,𝐮_ISE]^𝑇, and (20) becomes 𝜹_𝑙=𝐋⁽¹⁾𝐍^∗∗𝐔^∗∗. (21) The interface stress of 𝑠^𝑜_𝑗 is 𝐓_𝑙=𝐃𝜹_𝑙=𝐃𝐋⁽¹⁾𝐍^∗∗𝐔^∗∗. (22) 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 𝑊₀=_𝑆0𝑗12𝜹^𝑇_𝑙𝐃𝜹_𝐿𝑑𝑠. (23) The stiffness matrix 𝐾⁰_𝑗 on 𝑠^𝑜_𝑗 accordingly is 𝐤⁰_𝑗=_𝑆0𝐽𝐍^∗∗𝑇𝐋^(1)𝑇𝐃𝐋⁽¹⁾𝐍^∗∗𝑑𝑠. (24) The global stiffness matrix and load array of the mixed model are 𝐊=_𝑗𝐂^∗_𝑒^𝑇𝐤_𝑗𝐂^∗_𝑒+_𝑗0𝐂_𝑒^∗∗𝑇𝐤⁰_𝑗𝐂_𝑒^∗∗+_𝑛𝑒𝐂^𝑇_𝑒𝐤𝐂_𝑒,𝐐=_𝑗𝐂^𝑇_𝑒𝐑^𝑒_ISE+_𝑛𝑒𝐂^𝑇_𝑒𝐑^𝑒_FE, (25) respectively, where 𝐂^∗_𝑒, 𝐂_𝑒^∗∗, and 𝐂_𝑒 are the choice matrices of ISE, transitional interface element, and FE, respectively, and 𝐑^𝑒_ISE and 𝐑^𝑒_FE 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 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 𝐸=2.0×10⁴ 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.

Figure 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 𝛼=𝑡𝑔⁻¹(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³ Kg/m³ and that of the liquid be 𝛾=1.0×10³ Kg/m³. 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.

Figure 6 

Wedge and the mesh.

Figure 7 

FEM mesh.

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.

Table 2 

Normal stress 𝜎_𝑦 of section 𝑦 = 4 m (MPa).

Table 3 

Normal stress 𝜎_𝑦 of section 𝑦 = 8 m (MPa).

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.

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).