A 4-Mid-Node Plane Model of Base Force Element Method on Complementary Energy Principle

Abstract

Using the base forces as fundamental variables to describe the stress state and the displacement gradients that are the conjugate variables of the base forces to describe the deformation state for the two-dimensional elasticity problems, a 4-mid-node plane model of base force element method (BFEM) based on complementary energy principle is proposed. In this paper, the complementary energy of an element of the BFEM is constructed by using the base forces. The equilibrium conditions are released by the Lagrange multiplier method, and a modified complementary energy principle described by the base forces is obtained. The formulation of the 4-mid-node plane element of the BFEM is derived by assuming that the stress is uniformly distributed on each edge of the plane elements. A procedure of the BFEM on complementary energy principle is developed using MATLAB language. The numerical results of examples show that this model of the BFEM has high precision and is free from mesh sensitivity. This model shows good performances.

1. Introduction

The finite element method based on the assumed displacement field has become a method of choice for the solution of a wide variety of problems in structural mechanics. The advantages of the finite element method include the efficient and accurate modeling of domains with complex geometric configurations and varying material properties and the capability of accurately analyzing both geometrically and materially nonlinear problems. Shortcomings of the displacement model of FEM have been observed in the analyses of certain classes of problems, such as the large deformation, the treatment of nearly incompressible materials, the bending of thin plates, and the moving boundary problems. The developments of finite elements and corresponding formulations have been a subject of extensive research, which have been summarized in textbooks by Zienkiewicz [1], Bathe [2, 3], and Cook [4].

The attempts to extend and generalize the complementary energy methods that have been proposed over the last almost 50 years for small deformation solid mechanics problems have led to the development of several complementary energy principles and corresponding finite element models for nonlinear elastic solid/structural mechanics problems [5–7].

Recently, based on the complementary energy principle, a hybrid stress-function element method is proposed by Fu et al. [8] and Cen et al. [9–12] using the mesh distortion immune technique for developing plane 8-node elements. It starts from the principle of minimum complementary energy and employs the fundamental analytical solutions of the Airy stress function as the trial functions (analytical trial function method). So as long as the element edges keep straight, the resulting 8-node element HSF-Q8 can also produce the exact solutions for quadric displacement fields, and even the element shape degenerates into triangle and concave quadrangle.

A new concept was introduced by Gao [13], who used the concept of the “base forces” to replace various stress tensors for the description of the stress state at a point. These base forces can be directly obtained from the strain energy. For large deformation problems, when the base forces were adopted, the derivation of basic formulae was simplified by Gao [14, 15]. Based on the concept of the base forces, precise expressions for stiffness and compliance matrices for the FEM were obtained by Gao [13]. The applications of the stiffness matrix to the plane problems of elasticity using the four-edge plane element and the polygonal element were researched by Peng et al. [16]. Using the concept of base forces as state variables, a three-dimensional formulation of base force element method (BFEM) on complementary energy principle was proposed by Peng and Liu [17] for geometrically nonlinear problems. And the new finite element method based on the concept of base forces was called the base force element method (BFEM) by Peng and Liu [17]. A three-dimensional model of base force element method (BFEM) on complementary energy principle was proposed by Liu and Peng [18] for elasticity problems. The application of 2D base force element method (BFEM) to geometrically nonlinear analysis was proposed by Peng et al. [19].

The objective of the present research is to present a new two-dimensional formulation of the base forces element method (BFEM) for elasticity problems. A 4-mid-node plane element model of the BFEM for elasticity problems will be derived. In the present formulations of the BFEM, the “base forces” are treated as unknown variables, and the basic equations are constructed by means of the complementary energy principle. The element equilibrium conditions are fulfilled using the Lagrange multiplier method. Explicit expressions for the control equations are provided, and a procedure of the present method is developed. A number of elasticity problems are solved using the present formulation, and the results are compared with corresponding analytical solutions.

2. Basic Equations

Consider a two-dimensional domain of solid medium, let 𝐱^𝛼 (𝛼=1,2) denote the Lagrangian coordinate system, where 𝐏 and 𝐐 are the position vectors of a material point before and after deformation, respectively. Two triads for original and current configurations can be defined as 𝐏_𝛼=𝜕𝐏𝜕𝑥^𝛼,𝐐_𝛼=𝜕𝐐𝜕𝑥^𝛼. (1) Let 𝐮 denote the displacement of a point; then 𝐮=𝐐−𝐏. (2)

The gradient of displacement 𝐮_𝛼 can be written as 𝐮_𝛼=𝜕𝐮𝜕𝑥^𝛼=𝐐_𝛼−𝐏_𝛼. (3)

In order to describe the stress state at a point 𝐐, a parallelogram with the edges d𝑥¹𝐐_𝟏 and d𝑥²𝐐₂ is shown in Figure 1. Define 𝐓^𝛼=d𝐓^𝛼d𝑥^𝛼+1,d𝑥^𝛼⟶0, (4) where 3=1 for indexes. Quantities 𝐓^𝛼 (𝛼=1,2) are called the base forces at point 𝐐 in the two-dimensional coordinate system 𝑥^𝛼.

Figure 1 

Base forces on plane problem.

The 𝐓^𝛼 and 𝐮_𝛼 can directly be expressed by strain energy. Thus, 𝐮_𝛼 is just the conjugate variable of 𝐓^𝛼. It can be seen that the mechanics problem can be completely established by means of 𝐓^𝛼 and 𝐮_𝛼.

According to the definitions of Cauchy stress tensors, the relation between the base forces and Cauchy stress tensors can be given as 𝝈=𝐴_𝑄⁻¹𝐓^𝛼⊗𝐐_𝛼, (5) where ⊗ is the dyadic symbol, the summation rule is implied, and 𝐴_𝑄 is the current base area 𝐴_𝑄=𝐐₁×𝐐₂. (6)

For the isotropic material, the complementary energy density can be expressed as follows: 𝑊_C=(2𝐸)⁻¹(1+𝜈)𝐽_2𝑇−𝜈𝐽²_1𝑇 (7) in which 𝐸 is Young’s modulus, 𝜈 is Poisson’s ratio, and 𝐽_1𝑇 and 𝐽_2𝑇 are the invariants of 𝐓^𝛼, so 𝐽_1𝑇=𝐴_𝑃⁻¹𝐓^𝛼⋅𝐏_𝛼,𝐽_2𝑇=𝐴²_𝑃⁻¹𝐓^𝛼⋅𝐓^𝛽𝑝_𝛼𝛽, (8) where 𝐴_𝑃 is the original base area, 𝑝_𝛼𝛽 is a metric tensor in the original configurations, and there are the following relations: 𝐴_𝑃=𝐏₁×𝐏₂,𝑝_𝛼𝛽=𝐏_𝛼⋅𝐏_𝛽. (9)

3. Compliance Matrix

Now, consider a 4-mid-node plane element as shown in Figure 2. Let 𝐼, 𝐽, 𝐾, and 𝐿 denote its edges and 𝐓^𝐼, 𝐓^𝐽, 𝐓^𝐾, and 𝐓^𝐿 the force vectors acting on each of the edges.

Figure 2 

Four-mid-node plane element.

Substituting (8) into (7), the complementary energy of an element can be reduced as follows: 𝑊^𝑒_𝐶=1+𝜈𝐓2𝐸𝐴^𝐼⋅𝐓^𝐽𝑝_𝐼𝐽−𝜈𝐓1+𝜈^𝐼⋅𝐏_𝐈², (10) where 𝐴 is the original area of the element, 𝐓^𝐼 is the force vectors acting on the center of the boundary edge 𝐼 of the element, 𝑝_𝐼𝐽 is the dot product of position vectors 𝐏_𝐼 and 𝐏_𝐽 of points 𝐼 and 𝐽, and there are 𝐓^𝐼=𝑇^𝐼1𝐞₁+𝑇^𝐼2𝐞₂,𝐓^𝐽=𝑇^𝐽1𝐞₁+𝑇^𝐽2𝐞₂𝑝_𝐼𝐽=𝐏_𝐼⋅𝐏_𝐉,𝐏_𝐼=𝑃_𝐼1𝐞₁+𝑃_𝐼2𝐞₂,𝐏_𝐽=𝑃_𝐽1𝐞₁+𝑃_𝐽2𝐞₂ (11)(12) in which 𝐞₁ and 𝐞₂ are the unit vectors.

From (10), we can obtain the generalized displacement corresponding to 𝐓^𝐼𝛿_𝐼=𝜕𝑊^𝑒_𝐶𝜕𝐓^𝐼=𝐂_𝐼𝐽⋅𝐓^𝐽, (13) where 𝐂_𝐼𝐽 is a second-order tensor that is called the compliance matrix 𝐂_𝐼𝐽=1+𝜈𝑝𝐸𝐴^𝐼𝐽𝜈𝐔−𝐏1+𝜈^𝐼⊗𝐏_𝐽(𝐼,𝐽=1,2,3,4), (14) where 𝐔 is the unit tensor, and it can be expressed as 𝐔=𝐞₁⊗𝐞₁+𝐞₂⊗𝐞₂. (15)

Substituting (12) and (15) into (14), the compliance matrix of an element can be reduced as follows: 𝐂_𝐼𝐽=1+𝜈1𝐸𝐴𝑃1+𝜈^𝐼1𝑃_𝐽1+𝑃_𝐼2𝑃_𝐽2𝐞₁⊗𝐞₁−𝜈𝑃1+𝜈^𝐼1𝑃_𝐽2𝐞₁⊗𝐞₂−𝜈𝑃1+𝜈^𝐼2𝑃_𝐽1𝐞₂⊗𝐞₁+𝑃_𝐼1𝑃_𝐽1+1𝑃1+𝜈^𝐼2𝑃_𝐽2𝐞₂⊗𝐞₂, (16) or 𝐂_𝐼𝐽=𝐶_𝐼1𝐽1𝐞₁⊗𝐞₁+𝐶_𝐼1𝐽2𝐞₁⊗𝐞₂+𝐶_𝐼2𝐽1𝐞₂⊗𝐞₁+𝐶_𝐼2𝐽2𝐞₂⊗𝐞₂, (17) or 𝐂_𝐼𝐽=1+𝜈×1𝐸𝐴𝑃1+𝜈^𝐼1𝑃_𝐽1+𝑃_𝐼2𝑃_𝐽2−𝜈𝑃1+𝜈^𝐼1𝑃_𝐽2−𝜈𝑃1+𝜈^𝐼2𝑃_𝐽1𝑃_𝐼1𝑃_𝐽1+1𝑃1+𝜈^𝐼2𝑃_𝐽2(𝐼,𝐽=1,2,3,4). (18)

For a strain problem, it is necessary to replace 𝐸 with 𝐸/(1−𝜈²) and 𝜈 with 𝜈/(1−𝜈) in (16)–(18).

4. Governing Equations

The total complementary energy of the elastic system can be written as Π_𝐶=_𝑛𝑊^𝑒_𝐶−𝐮_𝐼⋅𝐓^𝐼. (19)

According to the complementary energy principle, Π_𝐶^𝑛𝑒 takes stationary value under equilibrium conditions and stress boundary conditions for the real displacement and stress state.

The equilibrium conditions can be released by the Lagrange multiplier method. For the stress boundary condition and the corresponding condition of the face forces between elements, computational techniques can be applied to achieve a rotation. Using the Lagrange multiplier method, a new function for an element can be introduced as follows: Π^𝑒∗_𝐶𝐓,𝜆,𝜆₃=Π^𝑒_𝐶(𝐓)+𝜆⁴_𝐼=1𝐓^𝐼+𝜆₃𝐓^𝐼×𝐏_𝐼 (20) in which arbitrary vectors 𝜆, 𝜆₃ are the Lagrange multipliers, and 𝜆 can be expressed as 𝜆=𝜆₁𝐞₁+𝜆₂𝐞₂. (21)

For the elastic system, there is Π^∗_𝐶=_𝑛Π^𝑒_𝐶^∗𝐓,𝜆,𝜆₃. (22) And by means of the modified complementary energy principle 𝛿Π^∗_𝐶=_𝑛𝛿Π^𝑒_𝐶^∗𝐓,𝜆,𝜆₃=0. (23)

Further, (23) can be expressed as 𝜕Π^∗_𝐶𝐓,𝜆,𝜆₃𝜕𝐓=0,𝜕Π^∗_𝐶𝐓,𝜆,𝜆₃𝜕𝜆=0,𝜕Π^∗_𝐶𝐓,𝜆,𝜆₃𝜕𝜆₃=0. (24) Equations (24) are the compatibility equations and displacement boundary conditions for the elastic system. These are the governing equations of the BFEM.

An explicit expression of stress for the 4-mid-node element, according to (5), can be written as 𝜎=𝐴⁻¹𝐓^𝐼⊗𝐏_𝐼 (25) or 𝜎=𝐴⁴⁻¹_𝐼=1𝑇^𝐼1𝑃_𝐼1𝐞₁⊗𝐞₁+𝑇^𝐼1𝑃_𝐼2𝐞₁⊗𝐞₂+𝑇^𝐼2𝑃_𝐼1𝐞₂⊗𝐞₁+𝑇^𝐼2𝑃_𝐼2𝐞₂⊗𝐞₂ (26) or 1𝜎=𝐴⁴_𝐼∑₌₁𝑇^𝐼1𝑃⁴_𝐼1∑_𝐼=1𝑇^𝐼1𝑃_4𝐼2𝐼∑₌₁𝑇^𝐼2𝑃⁴_𝐼1∑_𝐼=1𝑇^𝐼2𝑃_𝐼2. (27)

An explicit expression for displacement of node, based on the governing equation of element, can be developed. For a 4-mid-node plane element, the governing equation can be written as 𝛿_𝐼=𝜕Π^𝑒_𝐶^∗𝐓,𝜆,𝜆₃𝜕𝐓^𝐼. (28)

Substituting (20) into (28) and considering (19), (10), and (14), the explicit expression for displacement of nodes can be written as 𝛿_𝐼=𝐂_𝐼𝐽⋅𝐓^𝐽+𝜆+𝜆₃𝜀⋅𝐏_𝐼, (29) where 𝜀 is the alternating tensor and 𝜀 can be written as 𝜀=𝜀^𝛼𝛽𝐞_𝛼⊗𝐞_𝛽 (30) or 𝜀=𝐞₁⊗𝐞₂−𝐞₂⊗𝐞₁. (31)

Substituting (11), (17), (21), and (31) into (29), the explicit expression for displacement of nodes can be written as 𝛿_𝐼=𝐶_𝐼1𝐽1𝑇^𝐽1+𝐶_𝐼1𝐽2𝑇^𝐽2+𝜆₁+𝜆₃𝑃_𝐼2𝐞₁+𝐶_𝐼2𝐽1𝑇^𝐽1+𝐶_𝐼2𝐽2𝑇^𝐽2+𝜆₂−𝜆₃𝑃_𝐼1𝐞₂ (32) or 𝛿_𝐼=𝐶_𝐼1𝐽1𝑇^𝐽1+𝐶_𝐼1𝐽2𝑇^𝐽2+𝜆₁+𝜆₃𝑃_𝐼2𝐶_𝐼2𝐽1𝑇^𝐽1+𝐶_𝐼2𝐽2𝑇^𝐽2+𝜆₂−𝜆₃𝑃_𝐼1. (33)

5. Numerical Examples

To verify the validity and accuracy of the 4-mid-node plane model of base force element method proposed in this paper, a program on the BFEM is made, and a number of the numerical examples on the BFEM for the typical plane elasticity problems are presented in this section. The results obtained from the 2D BFEM are compared with corresponding theoretical solutions. The responses obtained using the traditional displacement finite element method are also given for some problems in order to assess the advantages of the 2D BFEM. All calculations of the traditional displacement finite element method are performed using the four-node isoparametric element (Q4 model).

Example 1 (pure shear of rectangular plate). Consider the problem of a rectangular plate subjected to uniformly distributing shearing forces as shown in Figure 3. Parameters about the calculated specimen are taken as 𝑎=2.0 m, 𝑏=1.0 m, 𝐸=10⁵ MPa, 𝜈=0.3, and 𝜏₀=1.0 MPa for a plane stress problem. The calculation is done using two different element meshes with the center nodes of each of the edges as shown in Figure 4 successively.

Figure 3 

Rectangular plate subjected to shearing forces.

(a) Uniform 4-side elements

(b) Arbitrary 4-side elements

Figure 4 

Two kinds of meshes for the pure shear problem.

The numerical solutions are that the stress component of various elements is 𝜏_𝑥𝑦=1 for two kinds of different meshes. The results are consistent with the theoretical results [1].

Example 2 (thick cylinder with fixed outer surface under internal pressure). Consider a thick cylinder with inner radius 𝑎 and outer radius 𝑏 under inner pressure 𝑝_𝑖 as shown in Figure 5.

Figure 5 

Thick cylinder with fixed outer edge under internal pressure.

The calculated parameters are taken to be 𝑎=1.0m, 𝑏=5.0m, 𝐸=10⁶MPa, 𝜈=0.3, and 𝑝_𝑖=1.0MPa for a plane strain problem. The calculation is done using three different element meshes with the center nodes of each of the edges as shown in Figure 6 successively.

(a) 20 × 10 elements

(b) 20 × 7 elements

(c) 20 × 5 elements

Figure 6 

Three kinds of meshes for cylinder.

The values of stress components 𝜎_𝑟, 𝜎_𝜃, and 𝑢_𝑟 for various radiuses 𝑟 are listed in Tables 1, 2, and 3, respectively. Comparisons of the results from the theoretical solution [1] and traditional 4-node quadrilateral isoparametric element (Q4 model) are also given in Tables 1, 2, and 3, respectively. The numerical results of the present model are consistent with those of the theoretical solutions and have given very good performance compared with Q4 model.

Table 1 

The values of 𝜎_𝑟.

Table 2 

The values of 𝜎_𝜃.

Table 3 

The values of 𝑢_𝑟.

Example 3 (a cantilever under uniformly distributed load). It consists of a straight cantilever beam under uniformly distributed load, as shown in Figure 7. The beam has length 𝐿=10m, cross-sectional dimensions 𝑏=1m and ℎ=1m, elastic modulus 𝐸=10⁹ MPa, and Poisson’s ratio 𝜈=0, respectively. While the applied total load level is 𝑝=10 MPa. The calculated specimen was divided into the 4-mid-node base force elements with the center nodes of each of the edges, as shown in Figure 8.

Figure 7 

Cantilever beam under uniformly distributed load.

(a) 7 × 7 elements

(b) 14 × 7 elements

(c) 20 × 7 elements

Figure 8 

Meshes of cantilever beam under uniformly distributed load.

The vertical displacement solutions of BFEM at the tip of the beam with different meshes are listed in Table 4. The vertical displacement of BFEM at the tip of the beam is compared with that provided by traditional 4-node quadrilateral isoparametric element (Q4 model) and the analytical solution.

Table 4 

Vertical displacement of the cantilever beam under uniformly distributed load.

Table 4 shows that the results of the present formulation agree well with those of the analytical solution, and the 4-mid-node element of BFEM has given very good performance compared with Q4 model for the large aspect ratio of elements.

Example 4 (a cantilever with large aspect ratio under a concentrated end load). A clamped cantilever with a concentrated end load is carried out as shown in Figure 9. The beam has length 𝐿=5 m, cross-sectional dimensions 𝑏=1 m and ℎ=0.1 m, elastic modulus 𝐸=10⁹ N/m², and Poisson’s ratio 𝜈=0. While the applied total load level is 𝑃=1 N. The calculated specimen was divided into the 4-mid-node base force elements with the center nodes of each of the edges, as shown in Figure 10.

Figure 9 

Cantilever beam under a concentrated end load.

(a) 60 × 9 elements

(b) 30 × 9 elements

(c) 15 × 9 elements

Figure 10 

Meshes of cantilever beam under a concentrated end load.

The vertical displacement of BFEM at the tip of the beam is compared with that provided by traditional 4-node quadrilateral isoparametric element (Q4 model) and the analytical solution in Table 5. The results show that the 4-mid-node element of BFEM has given very good performance compared with Q4 model for the large aspect ratio of elements.

Table 5 

Vertical displacement of the cantilever beam with a concentrated end load.

Example 5 (a cantilever with mesh distortion under a concentrated end load). A clamped cantilever with a concentrated end load is carried out as shown in Figure 9. The beam has length 𝐿=10 m, cross-sectional dimensions 𝑏=1 m and ℎ=1 m, elastic modulus 𝐸=10⁶ Pa, and Poisson’s ratio 𝜈=0, while the applied total load level is 𝑃=1 N. The calculated specimen was divided into the 4-mid-node base force elements with the center nodes of each of the edges, as shown in Figure 11.

(a) Mesh distortion parameters 𝑑=0

(b) Mesh distortion parameters 𝑑=0.5

(c) Mesh distortion parameters 𝑑=1.0

(d) Mesh distortion parameters 𝑑=1.5

Figure 11 

Distortion meshes of cantilever beam.

The normalized vertical displacement of BFEM at the tip of the beam is compared with that provided by traditional 4-node quadrilateral isoparametric element (Q4 model) and the analytical solution in Table 6. The results show that the 4-mid-node element of BFEM has given very good performance compared with Q4 model for the distortion mesh problems.

Table 6 

Normalized vertical displacement with the effect of mesh sensitivity.

6. Conclusions

In this paper, a 4-mid-node plane model of base force element method (BFEM) on complementary energy principle is proposed for elasticity problems. Based on the complementary energy principle, the equilibrium conditions are released by the Lagrange multiplier method, and a modified complementary energy principle described by the base forces is obtained. The 4-mid-node plane element model is derived by assuming that the stress is uniformly distributed on each edge of a plane element. A BFEM procedure is developed using MATLAB language. The chief features of the method are that the model does not introduce an interpolating function and can be used in any coordinate system and is not necessary to introduce Gauss’ integral for calculating the stiffness coefficient at a point.

The present model is applied for the solution of some elasticity problems. The comparisons are made with the corresponding analytical solutions, and a good agreement of results is observed. A series of numerical experiments are performed in order to assess the relative performances of the present model; based on these studies, it is found that the 4-mid-node plane model of BFEM with complementary energy principle provides reliable predictions for elasticity problems. A good performance of the 4-mid-node plane element of BFEM is observed for some problems.

The 4-mid-node plane element model of BFEM is also compared with the standard displacement method. It has been demonstrated by the results of numerical calculation that the new BFEM formulations proposed in this paper are effective and accurate. The 4-mid-node plane element model can be used efficiently for the elasticity analysis with increasing element aspect ratios and distortion meshes. The results presented here revealed a better overall performance of the BFEM in calculations. These results provide further confirmation that the BFEM can be used successfully and efficiently in the structural analysis, thus serving as an alternative method, and certain applications, being the method of choice for the analysis.

Conflict of Interests

The authors declare that there is no conflict of interests.

Acknowledgment

This work is supported by the National Science Foundation of China (nos. 91016026, 10972015, and 11172015).