- Email: [email protected]
Development and implementation of new triangular finite element based on MGE theory for bi-material analysis
Results in Physics 13 (2019) 102231
Contents lists available at ScienceDirect
Results in Physics journal homepage: www.elsevier.com/locate/rinp
Development and implementation of new triangular finite element based on MGE theory for bi-material analysis Junbao Wanga,b, Weiwei Lia, Zhanping Songa,b, a b
T
⁎
School of Civil Engineering, Xi’an University of Architecture and Technology, Xi’an 710055, China Institute of Tunnel and Underground Structure Engineering, Xi’an University of Architecture and Technology, Xi’an 710055, China
ARTICLE INFO
ABSTRACT
Keywords: MGE theory Triangular finite element Shear boundary layer Internal characteristic length
Compared with the classical continuum theory, the modified gradient elastic theory (MGE theory) better explains the size effect, strain and damage localization of materials. In this paper, we constructed a C1 planar triangular finite element based on the MGE theory. We also derived its finite element format and compiled the corresponding finite element calculation program by the MATLAB software. The program was used to analyze the bi-material shear boundary layer problem with unidirectional infinite length. The results show that the change of the internal characteristic length of the material in the infinite direction (the same direction as the shear stress) has little effect on calculation results of the shear strain, while the change of the internal characteristic length in the finite direction has a certain impact on calculation results of the shear strain, and its impact is even greater on the material with small elastic modulus in the bi-material; when the size of the finite element mesh does not exceed the internal characteristic length of the material, the mesh dependence of calculation results can be eliminated. Compared with the MGE rectangular element, the triangular element constructed in this paper has higher calculation accuracy and better adaptability.
Introduction Experiments indicated that the mechanical properties of metallic materials have distinct size effect at the micrometer and nanoscale [1–3], and the strain and damage localization of geomaterials also have strong size effect [4,5]. However, the classical continuum theory has limitations in explaining the size effect of the mechanical properties of materials, and the numerical simulation for the size effect based on the classical continuum theory excessively depends on the mesh [6]. This is because in the classical continuum theory, the interaction among the microstructures of the materials is not taken into consideration, and it can only describe the macroscopic mechanical behavior of the material, failing to reflect its microscopic nature. Every material has its microstructure, and the macroscopic mechanical behavior of the material is a macroscopic reflection of its microstructural changes. Therefore, if an appropriate theory can be used to explain the macroscopic mechanical behavior of the material on a microscopic level, it is of great significance for the prediction of material damage [7–9]. The emergence of the strain gradient theory provides a new idea for solving the size effect [10–14]. As a non-local theory, the strain gradient theory is introduced with a strain gradient term in its basic equation to elevate the interaction between material particles. Cosserat ⁎
et al. [10] assumed that the material microparticles were rigid particles, and that the particles both generated displacement and rotation during the deformation process. Correspondingly, the coupled stress and microrotation components were introduced into the equilibrium equation of the classical elastic theory, but only the rotation gradient was taken into consideration in the theory. Mindlin and Tiersten [15] introduced the constants related to the microstructure of the material in the constitutive equation of the classical elastic isotropic material to reflect the gradient effect, and solve the stress concentration factor of the infinite central plate with a circular hole affected by the far-field uniform force. This theory is complicated in the formula and does not define the physical meaning of the newly introduced material constant. However, it points out the concept of the internal length parameter, which positively affects the development of the strain gradient theory. Most scholars later have focused on the definition of the internal characteristic length and its physical meaning. Among them, the typical one is the strain gradient elastic theory proposed by Aifantis et al. [16–18], which assumes that the strain energy density is determined by the strain tensor and the strain gradient tensor. Combined with the classical elastic constitutive theory, they derive the relevant formula. Fleck and Hutchinson [19] further developed the theory of coupled stress, deduced the strain gradient theory with a focus on the plasticity of
Corresponding author at: School of Civil Engineering, Xi’an University of Architecture and Technology, Yanta Road 13, Xi’an, Shaanxi 710055, China. E-mail address: [email protected] (Z. Song).
https://doi.org/10.1016/j.rinp.2019.102231 Received 11 December 2018; Received in revised form 11 March 2019; Accepted 22 March 2019 Available online 29 March 2019 2211-3797/ © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/).
Results in Physics 13 (2019) 102231
J. Wang, et al.
materials, and raised the analytical solution of the shear strain of the interface boundary layer. Shu and Fleck [20] studied the interface boundary layer with the finite element method, but the calculation results had mesh dependence. Huang and Zhang et al. [21] studied the interface boundary layer and proposed the CMSG theory. In recent years, scholars have proposed many gradient theories [21–25], but the numerical methods of these theories have the following shortcomings: (1) Some strain gradient theories have no numerical implementation methods; (2) Although some strain gradient theories have numerical implementation methods, their discrete elements are singular; (3) So far, there is no mature numerical simulation software for strain gradient theories. In the finite element numerical analysis of the gradient theory, Zervos [6] provided a finite element method based on the Mindlin microstructural linear elasticity theory, but this method does not eliminate finite element mesh dependence. Askes and Gutiérrez [26] proposed a numerical method based on the Eringen theory. However, it has a low displacement interpolation polynomial order and a slow convergence speed. Soh et al. [27] constructed a planar triangular element with 18 degrees of freedom based on the mesoscopic finite element elastic strain gradient theory. Based on the assumed strain finite element method of the strain gradient theory, Chen [28] analyzed the size effect of the micro-scale beam bending with the 8-node rectangular isoparametric element. Zhao et al. [29] established a strain gradient axisymmetric model and simulated the size effect of the material. Song et al. [30] proposed a modified strain gradient elastic theory (MGE theory), and constructed a C1 planar rectangular element with 4 nodes and 32 degrees of freedom based on the theory. Due to various gradient elasticity theories and diversified finite element structures, the corresponding finite element numerical analysis method needs further improvement and development. In this paper, we constructed a C1 planar triangular finite element (4 nodes and 20 degrees of freedom) based on the modified strain gradient elasticity theory [30], and derived its finite element format. We also compiled a corresponding finite element program by the MATLAB software. Finally, we verified the rationality of the program with numerical examples, analyzed the shear strain distribution of the bi-material shear boundary layer, and simulated the bi-material shear boundary layer effect.
ij
lke =
(k = x , y , z )
(3)
ijk
The strain energy density is developed using the Taylor series and retained to the quadratic term of ij and ijk . In the initial state of the material, ij = kl = mnq = 0 , so the strain energy density can be expressed as
= 0.5Dijkl
ij kl
+ Dijmn lqe (0.5lke
ijk mnq
+
ij mnq )
(4)
where Dijkl and Dijmn are elastic tensors. Substituting Eq. (4) into Eq. (2), we can obtain ij
= Dijmn (
ijk
+ lqe
mn
= Dijmn (
mn
mnq)
+ lqe
(5)
e mnq) lk
(6)
Eq. (5) is the constitutive equation of the MGE theory. When the internal characteristic length is lqe = 0 , the Eq. (5) degenerates into the classical elastic constitutive equation. Principle of virtual work According to Eqs. (1) and (2), the variation U of the total potential energy in the deformed body can be expressed as
U=
V
=
S
(
ij
ij
[nj ( (
V
+
ijk ) dV
ijk
ij, j ij, j
ijk , kj )
Ej (nk
ijk , kj )
ui dV +
ijk )
S
+ El (nl ) nk nj ijk nk nj E (
ijk ]
ui dS
ui ) dS
where V is the area occupied by the deformed body, S is the boundary of the deformed body, ui is the displacement component of the deformed body, nj is the outer normal unit vector on the boundary S of the deformed body, ij, j is the first-order Cauchy stress gradient, ijk, kj is the (·) is the tensor xk (·) nj nk ) x is a tensor k
second-order stress gradient, E (·) = nk
MGE theory
f V i
ui dV +
S
t¯i ui dS +
S
r¯i E ( ui ) dS
(9)
U= W
Substituting Eqs. (7) and (8) into equation (9), the equilibrium equation of the deformed body is
Constitutive equation
( ij,
ijk )
is
ijk ,
ij, j
(1)
nj (
ijk
= =
(10)
+ fi = 0
ij, j
ijk, kl )
Ej (nk
ijk )
+ El (nl ) nk nj
ijk
= t¯i
(11)
The high-order surface force boundary condition is ijk nk nj
= r¯i
(12)
The displacement boundary condition is
ij ijk
ijk, kj
The normal surface boundary condition is
In the isothermal and infinitesimal deformation process, it can be expressed as ij
(8)
where fi is the physical force of the deformed body, t¯i the normal surface force of the deformed body, and r¯i the high-order surface force of the deformed body. Since the temperature and speed changes are not considered, we can obtain from the virtual work principle of the deformed body [31,32] that
The MGE theory assumes that the strain energy density depends on the strain tensor and the strain gradient tensor [30]. It introduces the internal characteristic length vector related with the material microstructure into the constitutive equation to evaluate the influence of the material microstructure on the mechanical behavior of the material.
=
normal differ-
ential operator, and Ej (·) = ( jk tangential differential operator. The variation W of external force virtual work can be expressed as
W=
The MGE theory [30] assumes that the strain energy density determined by the strain tensor ij and the strain gradient tensor namely
(7)
ui = u¯ i
(2)
(13)
In Eqs. (11)–(13), u¯ i , r¯i and t¯i can be determined according to actual problems. Eqs. (11)–(13) constitute all boundary conditions of the deformed body.
where ij represents the Cauchy stress, and ijk the high-order stress. Define an internal characteristic length vectorlke , and 2
Results in Physics 13 (2019) 102231
J. Wang, et al.
Variational principle
L1 = L2 =
Define the potential energy functional of the strain gradient elastic (ui , ij, ijk ) , then material as
(u i ,
ij,
ijk ) = U
W
L3 =
(14)
ij,
ijk )
= U
W= S
V
dV
f V i
ui dV
(t¯i ui + r¯i ui, j ) dS = 0
(15)
According to Eq. (15), the total potential energy of the strain gradient elastic material is
(u i ,
ij,
ijk )
=
V
dV
f ui dV V i
S
(t¯i ui + r¯i ui, j ) dS
(17)
where A is the area of Δ123, and A1, A2 and A3 are the areas of ΔP23, ΔP13 and ΔP12, respectively. The ratios of L1, L2 and L3 are the area coordinates of the P point. In addition, the conversion relation between the area coordinates and the global coordinates is consistent with that in the definition of the classical constant strain triangle element [33]. The element's basic unknown parameter ae contains the displacement of each node and the displacement gradients at the three corners, and di (i = 1, 2, 3, 4) is the basic unknown parameter of node i.
According to Eq. (9), we can obtain
(u i ,
A1 A A2 A A3 A
di = {u, v, ui, x , vi, x , ui, y , vi, y } d 4 = {u , v } ae = {d1, d2, d3, d 4}
(16)
Obviously, Eqs. (10)–(12) are consistent with the standing-value (ui , ij, ijk ) = 0 ) of the potential energy functional condition ( (ui , ij, ijk ) required by the variational principle [33].
(18)
The conversion relation between the element’s basic unknown parameter ae and the overall basic unknown parameter a is as follows (19)
ae = Ga
where G is the transformation matrix. In order to meet the continuity of the C1 element shape function, the Hermite interpolation function [34] is used to determine the shape function of the triangular element. The triangular element shape function N constructed in this paper is expressed as follows
Triangular element programming based on the MGE theory In order to apply the MGE theory to finite element numerical calculations, Song et al. [30] constructed the corresponding C1 planar rectangular finite element (4 nodes and 32 degrees of freedom). Given that the rectangular element has poor adaptability to irregular boundaries, and in order to further enrich the finite element library based on the MGE theory, in this paper we developed a C1 planar triangular finite element based on the MGE theory (4 nodes and 20 degrees of freedom), and achieved its programming.
N = [ N1 I2 N1 I2
(20)
N1 I2 ]
where I2 is a second-order unit matrix, and
N1 = L12 (L1 + 3L 2 + 3L3)
Discretization and finite element format
N2 = L12 (c3 L2
c2 L3) + (c2
N3 = L12 (b2 L3
b3 L 2 ) + (b3
N4 = L22 (L2 + 3L3 + 3L1)
Fig. 1 shows a schematic diagram of the planar triangular element constructed herein. The 4 nodes of the element include the three corner points (points 1, 2, 3 in Fig. 1) and the centroid of the triangle (the point 4 in Fig. 1). Each corner point of the triangular element has 6 degrees of freedom (u , ui, x , ui, y , v , vi, x , vi, y ), and the centroid point has 2 degrees of freedom (u and v ) with a total degrees of freedom of 20. In Fig. 1, x and y are the global coordinates of the triangular element, L1, L2 and L3 are the area coordinates of the triangular element. Wherein, the area coordinates are defined as follows: in the triangular element of Fig. 1, the position of the any point P can be determined according to L1, L2 and L3 , and
N5 = L 22 (c1 L3
c3 L1) + (c3
N6 = L 22 (b3 L1
b1 L3) + (b1
N7 = L32 (L3 + 3L1 + 3L2 ) N8 = L32 (c2 L1
7L1 L 2 L3 c3) L1 L 2 L3 b2 ) L1 L 2 L3 7L1 L 2 L 3 c1) L1 L2 L3 b3) L1 L2 L3 7L1 L 2 L3
c1 L2) + (c1
c2) L1 L 2 L3
N9 = L32 (b1 L2 b2 L1) + (b2 N10 = 27L1 L2 L3
b1) L1 L 2 L3 (21)
Based on the above expression, the displacement of any point in the element can be expressed as follows (22)
u = Nae
With introduction of the differential operators L , L , B and B , the relation between the strain , the strain gradient and the element’s basic unknown parameter ae can be expressed as
3
= L u = B ae = L u = B ae
(23)
Only consider the two-dimensional plane, then
B =LN B =L N
4 1
(24)
where L and L can be expressed as
P
L =
2
x
0
Fig. 1. Schematic diagram of the triangular element. 3
y
y 2
L =
0
x2
0
(25a)
x
0
2
2
x y
x y
2
2
x y
x2
0
0
2
y2
2
2
y2
x y
(25b)
Results in Physics 13 (2019) 102231
J. Wang, et al.
F
F 1
3
2
4
5
7
9
8
10
H L (a)
6 (b)
Fig. 2. The beam structure : (a) Schematic diagram; (b) Discretization model.
Based on the definition above, the characteristic length lke is expressed as a matrix I, and the equation (15) can be expressed as (ae ) =
aT
GT ne
+ aT
GT ne
Ve
Ve
parameters; (Ⅳ) gnf: it stored structural load information; (Ⅴ) gbc: it stored structural boundary constraint information. 3) The finite element model was solved by the SolveModel subprogram. Firstly, the element stiffness matrix and the element load vector were solved to integrate the structural overall stiffness matrix K and the overall load vector P. On this basis, the displacement, strain and stress distribution of the structure under a given load were further obtained. 4) The solved results were post-processed by the Megpost subprogram to obtain the displacement, strain data and the corresponding cloud image, and analyze the finite element numerical of the structure.
(B T DB + B T IDI TB + 2B T DI TB ) dV G a NTfdV +
Se
(NTt¯ +
NTr¯ T ) dS = 0
(26)
where I is the matrix form of D is the elastic tensor matrix, and ne is the number of elements. Define the element stiffness matrix as K e and the element load vector Pe , we can obtain from Eq. (26) that
lke ,
Ke =
Pe =
Ve
Ve
(B T DB + B T IDI TB + 2B T DI TB ) dV
N TfdV +
Se
G TK e G a = ne
(N Tt¯ +
N Tr¯ T ) dS
G TPe ne
(27)
Case verification
(28)
Fig. 2(a) is a schematic view of a beam fixed at both ends. The beam span-height-ratio is 4, the span length (L) is 4000 mm, the height (H) is 1000 mm, and the thickness is 200 mm. We assumed the elastic modulus of the beam material is 2.1 × 105 MPa, the Poisson's ratio is 0.3, the internal characteristic length is 0, and its midspan is subjected to a vertical concentrated force F of 20 kN. MIDAS software is a finite element analysis software for structural design, which mainly includes Midas building, Midas Civil, Midas GTS and Midas XD. In order to verify the accuracy of the program developed in this paper, we analyzed the deflection of the beam shown in Fig. 2(a) using both the program developed in this paper and the Midas GTS software. Midas GTS software [36] is equipped to handle the entire range of geotechnical design applications including deep foundations, excavations, complex tunnel systems, seepage analysis, consolidation analysis, embankment design, dynamic and slope stability analysis. The discretization model of the beam structure and the node number are shown in Fig. 2(b) (the element internal nodes are not marked). Table 1 lists the deflection value of the each node calculated by the two methods, Fig. 3 shows the comparison of the values of each node. It can be seen from Table 1 and Fig. 3 that the results from both the softwares show that the deflection of the node at the middle of the span is the largest, while the deflection of the nodes at the end of the beam is
(29)
Thus, the expressions of the structural overall stiffness matrix K and the overall load vector P are
K= P=
ne
G TK e G
ne
G TPe
Ka = P
(30) (31)
Preparation of finite element program In section 3.1, the finite element format of the triangular element based on MGE theory was derived. Based on the finite element format, the corresponding calculation program was compiled by the MATLAB software. The process of programming includes the following steps: 1) The input file (the format is dat) including the node number, node coordinates, element number, element composition, material properties, node load, and constraints was obtained by pre-processing. The pre-processing phase is performed by the GiD software which is a universal, adaptive and user-friendly pre and post processor for numerical simulations in science and engineering [35]. The preprocessing specifically includes the following four steps: (Ⅰ) building the finite element model using the GiD software; (Ⅱ) exporting the finite element analysis of model information in text form by the GiD software; (III) modifying the derived model information to acquire the input file matching by the MATLAB software; (Ⅳ) reading the input file to the MATLAB software using the FemModel subprogram. 2) We used the FemModel subprogram to define the model and read the model information file (dat file), and obtained the following variables: (Ⅰ) gNode: it stored the node number and coordinate information; (Ⅱ) gElement: it stored the element number and its composition information; (III) gMaterial: it stored material
Table 1 The result of node deflection. Node number
Proposed program Y1 (×10−4 mm)
MIDAS Y2 (×10−4 mm)
Relative error (%)
1 2 3 4 5 6 7 8 9 10
0.00 0.00 1.63 1.64 3.47 3.04 1.63 1.64 0.00 0.00
0.00 0.00 1.65 1.61 3.57 2.96 1.65 1.61 0.00 0.00
0.00 0.00 1.76 2.17 2.75 2.65 1.76 2.17 0.00 0.00
Note: The relative error is calculated as |(Y2 − Y1)/Y2|. 4
Results in Physics 13 (2019) 102231
J. Wang, et al.
y
∞
Interface
2×10-2 m
σxy x
Material 1
Material 2
Analyzing domin
∞ Fig. 4. Bi-material shear boundary layer analysis model. Table 2 Material parameters.
Fig. 3. Comparison of calculated deflection of each node.
zero due to the displacement constraints. The closer the node is to the middle of the span, the greater the deflection. The highest error of 2.75% is obtained in node 5. If the size of the finite element is small enough, the deflection curve of the beam must be smooth with the center of the span symmetrically, and the deflection at the middle of the span is the largest. Therefore, the maximum error can occur only in the midspan of the beam, namely node 5. The errors are zero in nodes 1, 2, 9 and 10 because there are displacement constraints in these nodes (the deflection value of nodes 1,2,9 and 10 are all zero). Overall, the calculation results from the program developed in this paper are close to that of the Midas GTS. Therefore, the finite element program based on the MGE theory proposed in this paper is feasible.
Material number
Elasticity modulus E (MPa)
Poisson’s ratio μ
Characteristic length l
1 2
2 × 104 2 × 103
0.3 0.3
l1x , l1y l2x , l2y
c c c c
Fig. 5. Schematic diagram of the divided triangle mesh element.
Analysis of shear boundary layer calculation
c = 1.25 × 10−4 m. Take the internal characteristic lengths of the two materials in the y direction to be 0, that is, l1y = l2y = 0, and the internal characteristic lengths in the × direction are l1x = l2x = lx = 2.0 × 10-4 m, l1x = l2x = lx = 3.0 × 10-4 m, l1x = l2x = lx = 4.0 × 10-4 m and l1x = l2x = lx = 5.0 × 10-4 m. The influence of the internal characteristic length in the × direction on the calculation results is analyzed. Fig. 6 shows the shear strain distributions of the analyzed area at different internal characteristic lengths. In the figure, the point x = 0.01 m is the interface between the two materials. It can be
The bi-material shear boundary layer is a classic case study in the strain gradient theory. In order to further explain the rationality of the finite element program proposed in this paper, it was used to analyze the bi-material shear boundary layer issue. The effect of internal characteristic length on the calculation results The MGE theory relates the macroscopic deformation characteristics to the microstructure by introducing the internal characteristic length vector of the material, so it can reflect the influence of the material microstructure on the macroscopic mechanical behavior. In order to obtain the influence of the internal characteristic length of the material on the macroscopic mechanical properties, four different internal characteristic lengths were used to analyze the bi-material shear boundary layer. A schematic diagram of the selected analysis model is shown in Fig. 4. The model is a strip-like structure in which two strain gradient materials are bonded (regardless of gravity), with material 1 on the left and material 2 on the right. In the model, the y direction is infinitely long, the left side is a fixed end, and the right side is subjected to a uniform shear stress σxy = 50 MPa. The area between the two dashed lines in Fig. 4 was selected as the analysis area (0.02 × 0.004 m2). Material parameters [30] are shown in Table 2. Discretization (meshing) of the analyzed area is performed using the triangular element proposed in Section 3.1. A schematic diagram of the individual grids divided is shown in Fig. 5. The triangular mesh element in Fig. 5 is a right-angled isosceles triangle, and the size of the triangular element refers to the side length c of the right-angled side. Here, when dividing the mesh, the size of the triangular element is
lx = 2.0×10-4 m lx = 3.0×10-4 m
lx = 4.0×10-4 m
lx = 5.0×10-4 m
x(×10-2m) Fig. 6. Distributions of Shear strain with different lx. 5
Results in Physics 13 (2019) 102231
J. Wang, et al.
Fig. 7. The shear strain curves of the points on the x-axis: (a) Shear strain with different lx; (b) Shear strain curves with different ly.
seen from Fig. 6 that as the internal characteristic length increases (from 2.0 × 10−4 m to 5.0 × 10−4 m), the boundary layer thickness of the material on the left side of the interface (Material 1) changes a little, while the boundary layer thickness of the material on the right side of the interface (Material 2) is increasing. According to the shear strain calculation results of the points on the x-axis (y = 0), the shear strain curves with different lx can be plotted, as shown in Fig. 7(a). It can be seen from Fig. 7(a) that for material 1, the variation of the internal characteristic length has a certain influence in the area near the fixed end of the left side (area I). At the same position in the area, the shear strain decreases slightly as the length of the internal characteristic increases. For material 2, the internal characteristic length has a much greater influence on the strain value in the area close of the material interface (region II), and the material shear strain at the same position in the region decreases with the increase of the internal characteristic length. From the above analysis, we can see that when the MGE theory is used to calculate and analyze the bi-material shear boundary layer, the calculation result of the material with small elastic modulus is more sensitive to the value for the internal characteristic length. Assume the internal characteristic lengths of the two materials in the x direction are both 0, that is, l1x = l2x = 0, and the internal characteristic lengths in the y direction are as follows: l1y = l2y = ly = 2.0 × 10−4 m, l1y = l2y = ly = 3.0 × 10−4 m, l1y = l2y = ly = 4.0 × 10−4 m and l1y = l2y = ly = 5.0 × 10−4 m. Fig. 7(b) shows the shear strain curves with different ly. It can be seen that the shear strain curves for different internal characteristic lengths almost coincide. Therefore, the change in the characteristic length of the y-direction has no effect on the calculation result, which is mainly due to the infinite length of the material in the y direction. In summary, when the MGE theory is used to analyze the bi-material shear boundary layer problem with infinite length in one direction, the change of the internal characteristic length of the material in the infinite direction (in the same direction as the shear stress) has little effect on the shear strain calculation results, while the variation of the internal characteristic length in the finite length direction has a certain influence on the calculation results of the shear strain, and the impact on the material with small elastic modulus in the bi-material is greater.
is analyzed with the program proposed in the paper. The analytical model, material parameters, and analysis areas are shown in Fig. 4 and Table 2, respectively (same as Section 4.1). In order to analyze the mesh dependence of the calculation results, the analyzed regions were divided into five different sizes of meshes. The mesh shape is shown in Fig. 5. The mesh sizes (right-edge length is c) are 2.00 × 103 m (2 × 10 × 2, 40 elements), 1.00 × 10−3 m (4 × 20 × 2, 160 elements), 5.00 × 10−4 m (8 × 40 × 2, 640 elements), 2.50 × 10−4 m (16 × 80 × 2, 2560 elements) and 1.25 × 10-4m (32 × 160 × 2, 10,240 elements). The finite element meshes are shown in Fig. 8. According to the analysis in the Section 4.1, the variation of the characteristic length in the y direction has no effect on the calculation results of the bi-material shear boundary layer. Therefore, the internal characteristic lengths of the two materials in the y direction are both set to 0, that is, l1y = l2y = 0, and the internal characteristic lengths in the x direction are set to l1x = 5 × 10−5 m and l2x = 5 × 10−4 m respectively. Fig. 9(a) shows the shear strain distributions of the analyzed area with different mesh sizes. In the Fig. 9(a), the point x = 0.01 m is the interface between the two materials. It can be seen from Fig. 9(a) that as the mesh size decreases, the boundary layer thickness (the width of the region where the shear strain changes in the vicinity of the interface) of the material on the left side of the interface (Material 1) continuously reduces, and the thickness of the boundary layer of the material on the right side of the interface (Material 2) also decreases. But when the mesh size reduces to 5.00 × 10−4 m, the thickness of the boundary layer on the right side of the interface is basically no longer changed. According to the shear strain calculation results of the points on the x-axis (y = 0), the shear strain curves on the corresponding section (xaxis) of the analyzed area can be plotted, as shown in Fig. 9(b). Fig. 9(b) shows: (1) For material 1, the strain at the fixed end and the vicinity of the interface between the two materials (regions I and II) varies with the mesh size, which means there is a mesh dependence. It is because that the internal characteristic length of material 1 (l1x = 5 × 10-5 m) is smaller than the five mesh sizes (cmin = 1.25 × 10−4 m). (2) For material 2, when the mesh sizes are c = 2.00 × 10-3 m and c = 1.00 × 103 m, the strain value is different from those in the other three conditions within a certain range of interface attachments (regions III and IV); when c = 5.00 × 10-4 m, c = 2.50 × 10-4 m and c = 1.25 × 10-4 m, the curves in the Fig. 9(b) basically coincide, that is, there is no mesh dependency. This is because when c = 2.00 × 10-3 m and c = 1.00 × 10-3 m, the mesh sizes are larger than the internal characteristic length of material 2 (l2x = 5 × 10-4 m); when c = 5.00 × 10-4 m, c = 2.50 × 104 m and c = 1.25 × 10-4 m, the mesh sizes do not exceed the internal characteristic length of the material 2.
Mesh dependency analysis When the classical elastic theory is used to analyze the finite element calculation of the bi- material shear boundary layer problem, the calculation results have mesh dependence [20]. The mesh dependence means that the calculation results vary greatly with the size of the finite element mesh. Here, the mesh dependence of the finite element method 6
Results in Physics 13 (2019) 102231
J. Wang, et al.
Fig. 8. Finite element meshes: (a) c = 2.00 × 10−3 m; (b) c = 1.00 × 10−3 m; (c) c = 5.00 × 10−4 m; (d) c = 2.50 × 10−4 m; (e) c = 1.25 × 10−4 m.
Based on the above analysis, when the mesh size is larger than the internal characteristic length of the material, the calculation result has mesh dependence; when the mesh size does not exceed the internal characteristic length, the calculation result is independent of the mesh size. This is consistent with the findings of Song et al. [30]. Therefore, when calculating the bi-material shear boundary layer with the finite element method, it is recommended that the mesh size to be divided should not exceed the internal characteristic length value of the material.
Shu et al. [20], the internal characteristic length of the two materials is l1 = l2 = l = 5.0 × 10-4 m. Fig. 10 shows the shear strain results for the points on the x-axis (y = 0) calculated by the three methods. It can be seen from Fig. 10 that: (1) the shear strain obtained by the three methods is continuous at the interface of the two materials; (2) The calculation results of shear strain based on the MGE triangular element are basically consistent with the calculation results of analytical solutions deduced by Shu et al. [20] in the whole analyzed region, and there is small error within a certain range (0.008≦x≦0.015 m) near the interface of the two materials; (3) The shear strain calculation results based on MGE rectangular elements are significantly different from the calculation results of analytical solutions deduced by Shu et al. [20] in the whole analyzed region. In the areas where the shear strain evenly distributes (0.0025≦x≦0.085 m and 0.015≦x≦0.02 m), the shear strain value is about 1/2 of the analytical solution. It can be seen that the calculation accuracy of the triangular element constructed in this paper is higher than that of the rectangular element. Because the theory and the size of the two finite elements are both the same, the only difference between them is the displacement mode. The displacement mode of triangular element constructed in this paper is obtained by the first-order Hermite interpolation, while the displacement mode of rectangular element is obtained by the bidirectional Hermite interpolation. Therefore, the difference in precision between the two finite elements may be caused by the displacement mode.
Comparison of calculation methods In order to illustrate the rationality and accuracy of calculation results of the MGE-based triangular finite element constructed in this paper, the MGE triangular finite element, the MGE rectangular finite element [30] and the analytical solution proposed by Shu et al. [20] based on the Toupin-Mindlin coupled stress theory were used to compare and analyze the bi-material shear boundary layer in Fig. 4. Material parameters are shown in Table 2. When meshing, the size of the triangular element is taken as c = 1.25 × 10-4 m, and the size of the rectangular element is m = n = 1.25 × 10-4 m (m and n are the length and width of the rectangular element). In the MGE theory, the internal characteristic length of the two materials in the y direction is l1y = l2y = ly = 0, and the internal characteristic length in the × direction is l1x = l2x = lx = 5.0 × 10-4 m. In the analytical solutions deduced by c=2.00×10-3 m c=1.00×10-3 m c=5.00×10-4 m c=2.50×10-4 m c=1.25×10-4 m
x (×10-2 m)
(a)
(b)
Fig. 9. The shear strain calculation results with different mesh sizes: (a) Distributions of Shear strain; (b) Shear strain curves(x-axis). 7
Results in Physics 13 (2019) 102231
J. Wang, et al.
Fig. 10. Comparison of shear strain calculation results.
Fig. 11. Effect of ly on the shear strain curves (x-axis).
Discussions
the rectangular element, the triangular element constructed in this paper has the following advantages: (1) As can be seen from the comparison in Fig. 10, the shear strain results based on the triangular element constructed in this paper are in the better agreement with analytic solutions. (2) The triangular element has better adaptability to the boundary of the analyzed material than the rectangular element, especially when the material boundary is complicated (such as arc or irregular), it has more advantage in the meshing. (3) Since the triangular element has smaller number of degrees of freedom, its stiffness matrix is of the lower order when compared with the rectangular element, and its calculation is faster. Thus, the calculation efficiency of the triangular element is higher in the case of the same algorithm. Many scholars have carried out a lot of research on the MGE theory and its finite element implementation, but there are still some shortcomings in the current research, which are mainly in: (1) The internal characteristic length is an important parameter in the strain gradient theory, but its physical meaning is still not clearly and uniformly explained. (2) It is unclear whether the internal characteristic length parameter is the intrinsic property of the material or changing with the deformation of the material. (3) The existing method for inverting the internal characteristic length parameter of the material through experiments needs high experimental conditions, which is not universally applicable. (4) At present, the research on strain gradient theory is still in theoretical and experimental research, and the theory is not well applied to engineering practice. Therefore, the current research on the strain gradient theory and its programmatic is still not perfect, and there is still a lot of research work to be done in this area.
For the bi-material shear boundary layer problem, the classical elastic theory holds that the stress continuously and evenly distributes at the interface of the two materials, and the strain is abrupt, which is obviously inconsistent with the actual condition [19–21,37,38]. Therefore, the finite element analysis of the bi-material shear boundary layer problem based on the classical elasticity theory inevitably has mesh dependence. The MGE theory relates the macroscopic deformation of the material to its internal microstructure by introducing the internal characteristic length vector, so it can reflect the influence of the microstructure of the material on the mechanical behavior. The calculation results in this paper (Figs. 7 and 9) show that the shear strain at the interface of the two materials calculated by the MGE theory is continuous. Therefore, compared with the solution given by the classical elastic theory, the calculation results of the bi-material shear boundary layer given by the MGE theory are obviously more in line with the actual situation. In addition, Shu et al. [20] gave an analytical solution of the shear strain of the bi-material shear boundary layer based on the ToupinMindlin coupled stress theory. But since the internal characteristic length in the Toupin-Mindlin coupled stress theory is scalar, it makes the analytical solution fail to reflect the influence of the anisotropy of the microstructure inside the material on the mechanical behavior. For the MGE theory, although the shear calculation results (Fig. 7(b)) of the bi-material model with infinite length in the y direction in Fig. 4 show that the variation of the internal characteristic length in the y direction has no effect on the calculation results of the shear strain, the size of the material in one direction cannot be infinitely long. Take the length of the two materials in the y direction of Fig. 4 as 0.004 m (finite length), and the rest of the conditions are unchanged, the material parameters are shown in Table 2; take the internal characteristic lengths of the two materials in the × direction are both 5.0 × 10-4 m (l1x = l2x = 5.0 × 10-4 m), and the internal characteristic lengths of the two materials in the y direction are also the same (l1y = l2y = ly). Fig. 11 shows the effect of variation of the internal characteristic length in the y direction on shear strain calculation results of the points on the x axis. It can be seen that the variation of the internal characteristic length in the y direction has a certain influence on the shear strain calculation result of the bimaterial shear boundary layer, and the influence on the material with small elastic modulus (material 2) is more obvious. Since the internal characteristic length of the material in the MGE theory is vector, this theory can reflect the influence of the anisotropy of the microstructure inside the material on the mechanical behavior. The triangular element and the rectangular element [30] based on the MGE theory are both high-order finite elements, but compared with
Conclusions 1) In this paper, a triangular finite element is constructed based on the MGE theory, and its finite element format is derived. We also compiled the corresponding finite element program by the MATLAB software. The program was used to analyze the bi-material shear boundary layer problem with unidirectional infinite length. 2) The variation of the internal characteristic length of the material in the infinite length direction has little effect on the shear strain calculation results, and the variation of the internal characteristic length in the finite length direction has a certain influence on the shear strain calculation result, and the impact on the material with small elastic modulus in the bi-material is greater. 3) When the finite element mesh size is larger than the internal characteristic length of the material, the calculation result has mesh dependence, and the mesh dependence can be eliminated when the mesh size is smaller than the characteristic length. 4) The shear strain calculated by the MGE triangular element is in good agreement with the analytical solutions based on the Toupin8
Results in Physics 13 (2019) 102231
J. Wang, et al.
Mindlin coupled stress theory, and the calculation result of the MGE rectangular element is only about 1/2 of the analytical solutions in the uniformly distributed area of the shear strain, indicating that the MGE triangle element has higher calculation accuracy.
[18] Askes H, Aifantis EC. Gradient elasticity in statics and dynamics: an overview of formulations, length scale identification procedures, finite element implementations and new results. Int J Solids Struct 2011;48:1962–90. https://doi.org/10. 1016/j.ijsolstr.2011.03.006. [19] Fleck NA, Hutchinson JW. A phenomenological theory for strain gradient effects in plasticity. J Mech Phys Solids 2001;41:1825–57. https://doi.org/10.1016/00225096(93)90072-n. [20] Shu JY, King WE, Fleck NA. Finite elements for materials with strain gradient effects. Int J Numer Methods Eng 2015;44:373–91. https://doi.org/10.1002/(sici) 1097-0207(19990130)44:3<373::aid-nme508>3.0.co;2-7. [21] Zhang ZH, Zhuang Z, Gao Y, Liu ZL, Nie JF. Cyclic plastic behavior analysis based on the micromorphic mixed hardening plasticity model. Comput Mater Sci 2011;50:1136–44. https://doi.org/10.1016/j.commatsci.2010.11.013. [22] Yang F, Chong ACM, Lam DCC, Tong P. Couple stress based strain gradient theory for elasticity. Int J Solids Struct 2002;39:2731–43. https://doi.org/10.1016/S00207683(02)00152-X. [23] Mohammad-Abadi M, Daneshmehr AR. Modified couple stress theory applied to dynamic analysis of composite laminated beams by considering different beam theories. Int J Eng Sci 2015;87:83–102. https://doi.org/10.1016/j.ijengsci.2014. 11.003. [24] Lim CW, Zhang G, Reddy JN. A higher-order nonlocal elasticity and strain gradient theory and its applications in wave propagation. J Mech Phys Solids 2015;78:298–313. https://doi.org/10.1016/j.jmps.2015.02.001. [25] Saroukhani S, Vafadari R, Simone A. A simplified implementation of a gradientenhanced damage model with transient length scale effects. Comput Mech 2013;51:899–909. https://doi.org/10.1007/s00466-012-0769-8. [26] Askes H, Gutiérrez MA. Implicit gradient elasticity. Int J Numer Methods Eng 2010;67:400–16. https://doi.org/10.1002/nme.1640. [27] Soh AK, Wanji C. Finite element formulations of strain gradient theory for microstructures and the C0–1 patch test. Int J Numer Methods Eng 2004;61:433–54. https://doi.org/10.1002/nme.1075. [28] Chen Z. Assumed strain finite element method based on the theory of strain gradient. J Univ Sci Technol Beijing 2007;29:807–10. https://doi.org/10.1016/S10016058(07)60030-4. [29] Zhao J, Chen WJ, Bin JI. Formula derivation in axisymmetric strain gradient theory and finite element implementation. Chin J Computat Mechan 2009;26:209–14. https://doi.org/10.1109/CLEOE-EQEC.2009.5194697. [30] Song Z, Bing Z, He J, Zheng Y. Modified gradient elasticity and its finite element method for shear boundary layer analysis. Mech Res Commun 2014;62:146–54. https://doi.org/10.1016/j.mechrescom.2014.09.008. [31] Wunderlich W, Pilkey WD. Mechanics of structures. Variational and computational methods. Meccanica 2004;39:291–2. https://doi.org/10.1023/B:MECC. 0000023038.64148.bc. [32] Li X, Zhang J, Zheng Y. Static and free vibration analysis of laminated composite plates using isogeometric approach based on the third order shear deformation theory. 232019-232019 Adv Mech Eng 2015;6. https://doi.org/10.1155/2014/ 232019. [33] Zienkiewicz OC, Nakazawa S. On variational formulation and its modifications for numerical solution. Comput Struct 1984;19:303–13. https://doi.org/10.1016/ 0045-7949(84)90231-1. [34] Meleshkina AV. On the approximation of the derivatives of the Hermite interpolation polynomial on a triangle. Zh Vychisl Mat Mat Fiz 2010;50:201–10. https:// doi.org/10.1134/S0965542510020016. [35] García-Doñoro D, García-Castillo LE, Gómez-Revuelto I. An interface between an – adaptive finite element package and the pre- and post-processor GiD. Finite Elem Anal Des 2010;46:328–38. https://doi.org/10.1016/j.finel.2009.11.005. [36] Wang J, Huo Q, Song Z, Zhang Y. Study on adaptability of primary support arch cover method for large-span embedded tunnels in the upper-soft lower-hard stratum. Adv Mech Eng 2019;11:1–15. https://doi.org/10.1177/ 1687814018825375. [37] Liu X, Han Y, Li D, Tu Y, Deng Z, Yu C, et al. Anti-pull mechanisms and weak interlayer parameter sensitivity analysis of tunnel-type anchorages in soft rock with underlying weak interlayers. Eng Geol 2019;253:123–36. https://doi.org/10.1016/ j.enggeo.2019.03.012. [38] Lai J, Qiu J, Fan H, Zhang Q, Hu Z, Wang J, et al. Fiber bragg grating sensors-based in situ monitoring and safety assessment of loess tunnel. J Sens 2016;2016:1–10. https://doi.org/10.1155/2016/8658290.
Acknowledgments The authors gratefully acknowledge the support received from National Natural Science Foundation of China (Grant Nos. 51578447, 51404184), China, the Science and Technology Project of Ministry of Housing Urban-Rural Construction (Grant No. 2017-K4-032), China, the Youth Science and Technology Nova Program of Shaanxi Province (Grant No. 2018KJXX-061), China and Natural Science Basic Research Program of Shaanxi Province (Grant No. 2016JQ4009), China. References [1] Huang Y, Zhang F, Hwang KC, Nix WD, Pharr GM, Feng G. A model of size effects in nano-indentation. J Mech Phys Solids 2016;54:1668–86. 0.1016/ j.jmps.2006.02.002. [2] Xu J, Dai CY, Zhang B, Zhang GP. Strain-gradient dependent fatigue behavior of micron-thick copper single crystal foils. Comput Mater Sci 2014;85:223–9. https:// doi.org/10.1016/j.commatsci.2014.01.004. [3] Sahmani S, Aghdam MM. Nonlinear primary resonance of micro/nano-beams made of nanoporous biomaterials incorporating nonlocality and strain gradient size dependency. Results Phys 2018;8. DOI:10.1016/j.rinp.2018.01.002. [4] Wang J, Song Z, Zhao B, Liu X, Liu J, Lai J. A study on the mechanical behavior and statistical damage constitutive model of sandstone. Arab J Sci Eng 2018;43:5179–92. https://doi.org/10.1007/s13369-017-3016-y. [5] Hoover CG, Bažant ZP. Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests. Int J Fract 2014;187:133–43. https://doi.org/10.1007/s10704-013-9926-0. [6] Zervos A. Some finite elements for elasticity with microstructure and gradient elasticity. Int J Numer Methods Eng 2010;73:564–95. https://doi.org/10.1002/ nme.2093. [7] Lai JX, Mao S, Qiu JL, Fan HB, Zhang Q, Hu ZN, et al. Investigation progresses and applications of fractional derivative model in geotechnical engineering. Math Probl Eng 2016;3:1–15. https://doi.org/10.1155/2016/9183296. [8] Qiu J, Qin Y, Lai J, Wang K, Niu F, Wang H, et al. Structural response of the metro tunnel under local dynamic water environment in loess strata. Geofluids 2019;2019:16. https://doi.org/10.1155/2019/8541959. [9] Wang J, Ren Z, Song Z, Huo R, Yang T. Study of the effect of micro-pore characteristics and saturation degree on the longitudinal wave velocity of sandstone. Arab J Geosci 2019. in presss. [10] Cosserat E, Cosserat F. Théorie des corps déformables. Nature 1909;81:67. https:// doi.org/10.1038/081067a0. [11] Borg U. A strain gradient crystal plasticity analysis of grain size effects in polycrystals. Eur J Mech A-Solids 2007;26:313–24. https://doi.org/10.1016/j. euromechsol.2006.09.006. [12] Farajpour A, Rastgoo A. Influence of carbon nanotubes on the buckling of microtubule bundles in viscoelastic cytoplasm using nonlocal strain gradient theory. Results Phys 2017;7:1367–75. https://doi.org/10.1016/j.rinp.2017.03.038. [13] Lazopoulos KA, Lazopoulos AK. Fractional derivatives and strain gradient elasticity. Acta Mech 2016;227:823–35. https://doi.org/10.1007/s00707-015-1489-x. [14] Tsagrakis I, Konstantinidis A, Aifantis EC. Strain gradient and wavelet interpretation of size effects in yield and strength. Mech Mater 2003;35:733–45. https://doi. org/10.1016/s0167-6636(02)00205-3. [15] Mindlin RD. Micro-structure in linear elasticity. Arch Ration Mech Anal 1964;16:51–78. https://doi.org/10.1007/bf00248490. [16] Konstantopoulos I, Aifantis E. Gradient elasticity applied to a crack. J Mech Behav Biomed Mater 2013;22:193–201. https://doi.org/10.1515/jmbm-2013-0026. [17] Aifantis EC. Update on a class of gradient theories. Mech Mater 2003;35:259–80. https://doi.org/10.1016/s0167-6636(02)00278-8.
9
Contents lists available at ScienceDirect
Results in Physics journal homepage: www.elsevier.com/locate/rinp
Development and implementation of new triangular finite element based on MGE theory for bi-material analysis Junbao Wanga,b, Weiwei Lia, Zhanping Songa,b, a b
T
⁎
School of Civil Engineering, Xi’an University of Architecture and Technology, Xi’an 710055, China Institute of Tunnel and Underground Structure Engineering, Xi’an University of Architecture and Technology, Xi’an 710055, China
ARTICLE INFO
ABSTRACT
Keywords: MGE theory Triangular finite element Shear boundary layer Internal characteristic length
Compared with the classical continuum theory, the modified gradient elastic theory (MGE theory) better explains the size effect, strain and damage localization of materials. In this paper, we constructed a C1 planar triangular finite element based on the MGE theory. We also derived its finite element format and compiled the corresponding finite element calculation program by the MATLAB software. The program was used to analyze the bi-material shear boundary layer problem with unidirectional infinite length. The results show that the change of the internal characteristic length of the material in the infinite direction (the same direction as the shear stress) has little effect on calculation results of the shear strain, while the change of the internal characteristic length in the finite direction has a certain impact on calculation results of the shear strain, and its impact is even greater on the material with small elastic modulus in the bi-material; when the size of the finite element mesh does not exceed the internal characteristic length of the material, the mesh dependence of calculation results can be eliminated. Compared with the MGE rectangular element, the triangular element constructed in this paper has higher calculation accuracy and better adaptability.
Introduction Experiments indicated that the mechanical properties of metallic materials have distinct size effect at the micrometer and nanoscale [1–3], and the strain and damage localization of geomaterials also have strong size effect [4,5]. However, the classical continuum theory has limitations in explaining the size effect of the mechanical properties of materials, and the numerical simulation for the size effect based on the classical continuum theory excessively depends on the mesh [6]. This is because in the classical continuum theory, the interaction among the microstructures of the materials is not taken into consideration, and it can only describe the macroscopic mechanical behavior of the material, failing to reflect its microscopic nature. Every material has its microstructure, and the macroscopic mechanical behavior of the material is a macroscopic reflection of its microstructural changes. Therefore, if an appropriate theory can be used to explain the macroscopic mechanical behavior of the material on a microscopic level, it is of great significance for the prediction of material damage [7–9]. The emergence of the strain gradient theory provides a new idea for solving the size effect [10–14]. As a non-local theory, the strain gradient theory is introduced with a strain gradient term in its basic equation to elevate the interaction between material particles. Cosserat ⁎
et al. [10] assumed that the material microparticles were rigid particles, and that the particles both generated displacement and rotation during the deformation process. Correspondingly, the coupled stress and microrotation components were introduced into the equilibrium equation of the classical elastic theory, but only the rotation gradient was taken into consideration in the theory. Mindlin and Tiersten [15] introduced the constants related to the microstructure of the material in the constitutive equation of the classical elastic isotropic material to reflect the gradient effect, and solve the stress concentration factor of the infinite central plate with a circular hole affected by the far-field uniform force. This theory is complicated in the formula and does not define the physical meaning of the newly introduced material constant. However, it points out the concept of the internal length parameter, which positively affects the development of the strain gradient theory. Most scholars later have focused on the definition of the internal characteristic length and its physical meaning. Among them, the typical one is the strain gradient elastic theory proposed by Aifantis et al. [16–18], which assumes that the strain energy density is determined by the strain tensor and the strain gradient tensor. Combined with the classical elastic constitutive theory, they derive the relevant formula. Fleck and Hutchinson [19] further developed the theory of coupled stress, deduced the strain gradient theory with a focus on the plasticity of
Corresponding author at: School of Civil Engineering, Xi’an University of Architecture and Technology, Yanta Road 13, Xi’an, Shaanxi 710055, China. E-mail address: [email protected] (Z. Song).
https://doi.org/10.1016/j.rinp.2019.102231 Received 11 December 2018; Received in revised form 11 March 2019; Accepted 22 March 2019 Available online 29 March 2019 2211-3797/ © 2019 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/BY-NC-ND/4.0/).
Results in Physics 13 (2019) 102231
J. Wang, et al.
materials, and raised the analytical solution of the shear strain of the interface boundary layer. Shu and Fleck [20] studied the interface boundary layer with the finite element method, but the calculation results had mesh dependence. Huang and Zhang et al. [21] studied the interface boundary layer and proposed the CMSG theory. In recent years, scholars have proposed many gradient theories [21–25], but the numerical methods of these theories have the following shortcomings: (1) Some strain gradient theories have no numerical implementation methods; (2) Although some strain gradient theories have numerical implementation methods, their discrete elements are singular; (3) So far, there is no mature numerical simulation software for strain gradient theories. In the finite element numerical analysis of the gradient theory, Zervos [6] provided a finite element method based on the Mindlin microstructural linear elasticity theory, but this method does not eliminate finite element mesh dependence. Askes and Gutiérrez [26] proposed a numerical method based on the Eringen theory. However, it has a low displacement interpolation polynomial order and a slow convergence speed. Soh et al. [27] constructed a planar triangular element with 18 degrees of freedom based on the mesoscopic finite element elastic strain gradient theory. Based on the assumed strain finite element method of the strain gradient theory, Chen [28] analyzed the size effect of the micro-scale beam bending with the 8-node rectangular isoparametric element. Zhao et al. [29] established a strain gradient axisymmetric model and simulated the size effect of the material. Song et al. [30] proposed a modified strain gradient elastic theory (MGE theory), and constructed a C1 planar rectangular element with 4 nodes and 32 degrees of freedom based on the theory. Due to various gradient elasticity theories and diversified finite element structures, the corresponding finite element numerical analysis method needs further improvement and development. In this paper, we constructed a C1 planar triangular finite element (4 nodes and 20 degrees of freedom) based on the modified strain gradient elasticity theory [30], and derived its finite element format. We also compiled a corresponding finite element program by the MATLAB software. Finally, we verified the rationality of the program with numerical examples, analyzed the shear strain distribution of the bi-material shear boundary layer, and simulated the bi-material shear boundary layer effect.
ij
lke =
(k = x , y , z )
(3)
ijk
The strain energy density is developed using the Taylor series and retained to the quadratic term of ij and ijk . In the initial state of the material, ij = kl = mnq = 0 , so the strain energy density can be expressed as
= 0.5Dijkl
ij kl
+ Dijmn lqe (0.5lke
ijk mnq
+
ij mnq )
(4)
where Dijkl and Dijmn are elastic tensors. Substituting Eq. (4) into Eq. (2), we can obtain ij
= Dijmn (
ijk
+ lqe
mn
= Dijmn (
mn
mnq)
+ lqe
(5)
e mnq) lk
(6)
Eq. (5) is the constitutive equation of the MGE theory. When the internal characteristic length is lqe = 0 , the Eq. (5) degenerates into the classical elastic constitutive equation. Principle of virtual work According to Eqs. (1) and (2), the variation U of the total potential energy in the deformed body can be expressed as
U=
V
=
S
(
ij
ij
[nj ( (
V
+
ijk ) dV
ijk
ij, j ij, j
ijk , kj )
Ej (nk
ijk , kj )
ui dV +
ijk )
S
+ El (nl ) nk nj ijk nk nj E (
ijk ]
ui dS
ui ) dS
where V is the area occupied by the deformed body, S is the boundary of the deformed body, ui is the displacement component of the deformed body, nj is the outer normal unit vector on the boundary S of the deformed body, ij, j is the first-order Cauchy stress gradient, ijk, kj is the (·) is the tensor xk (·) nj nk ) x is a tensor k
second-order stress gradient, E (·) = nk
MGE theory
f V i
ui dV +
S
t¯i ui dS +
S
r¯i E ( ui ) dS
(9)
U= W
Substituting Eqs. (7) and (8) into equation (9), the equilibrium equation of the deformed body is
Constitutive equation
( ij,
ijk )
is
ijk ,
ij, j
(1)
nj (
ijk
= =
(10)
+ fi = 0
ij, j
ijk, kl )
Ej (nk
ijk )
+ El (nl ) nk nj
ijk
= t¯i
(11)
The high-order surface force boundary condition is ijk nk nj
= r¯i
(12)
The displacement boundary condition is
ij ijk
ijk, kj
The normal surface boundary condition is
In the isothermal and infinitesimal deformation process, it can be expressed as ij
(8)
where fi is the physical force of the deformed body, t¯i the normal surface force of the deformed body, and r¯i the high-order surface force of the deformed body. Since the temperature and speed changes are not considered, we can obtain from the virtual work principle of the deformed body [31,32] that
The MGE theory assumes that the strain energy density depends on the strain tensor and the strain gradient tensor [30]. It introduces the internal characteristic length vector related with the material microstructure into the constitutive equation to evaluate the influence of the material microstructure on the mechanical behavior of the material.
=
normal differ-
ential operator, and Ej (·) = ( jk tangential differential operator. The variation W of external force virtual work can be expressed as
W=
The MGE theory [30] assumes that the strain energy density determined by the strain tensor ij and the strain gradient tensor namely
(7)
ui = u¯ i
(2)
(13)
In Eqs. (11)–(13), u¯ i , r¯i and t¯i can be determined according to actual problems. Eqs. (11)–(13) constitute all boundary conditions of the deformed body.
where ij represents the Cauchy stress, and ijk the high-order stress. Define an internal characteristic length vectorlke , and 2
Results in Physics 13 (2019) 102231
J. Wang, et al.
Variational principle
L1 = L2 =
Define the potential energy functional of the strain gradient elastic (ui , ij, ijk ) , then material as
(u i ,
ij,
ijk ) = U
W
L3 =
(14)
ij,
ijk )
= U
W= S
V
dV
f V i
ui dV
(t¯i ui + r¯i ui, j ) dS = 0
(15)
According to Eq. (15), the total potential energy of the strain gradient elastic material is
(u i ,
ij,
ijk )
=
V
dV
f ui dV V i
S
(t¯i ui + r¯i ui, j ) dS
(17)
where A is the area of Δ123, and A1, A2 and A3 are the areas of ΔP23, ΔP13 and ΔP12, respectively. The ratios of L1, L2 and L3 are the area coordinates of the P point. In addition, the conversion relation between the area coordinates and the global coordinates is consistent with that in the definition of the classical constant strain triangle element [33]. The element's basic unknown parameter ae contains the displacement of each node and the displacement gradients at the three corners, and di (i = 1, 2, 3, 4) is the basic unknown parameter of node i.
According to Eq. (9), we can obtain
(u i ,
A1 A A2 A A3 A
di = {u, v, ui, x , vi, x , ui, y , vi, y } d 4 = {u , v } ae = {d1, d2, d3, d 4}
(16)
Obviously, Eqs. (10)–(12) are consistent with the standing-value (ui , ij, ijk ) = 0 ) of the potential energy functional condition ( (ui , ij, ijk ) required by the variational principle [33].
(18)
The conversion relation between the element’s basic unknown parameter ae and the overall basic unknown parameter a is as follows (19)
ae = Ga
where G is the transformation matrix. In order to meet the continuity of the C1 element shape function, the Hermite interpolation function [34] is used to determine the shape function of the triangular element. The triangular element shape function N constructed in this paper is expressed as follows
Triangular element programming based on the MGE theory In order to apply the MGE theory to finite element numerical calculations, Song et al. [30] constructed the corresponding C1 planar rectangular finite element (4 nodes and 32 degrees of freedom). Given that the rectangular element has poor adaptability to irregular boundaries, and in order to further enrich the finite element library based on the MGE theory, in this paper we developed a C1 planar triangular finite element based on the MGE theory (4 nodes and 20 degrees of freedom), and achieved its programming.
N = [ N1 I2 N1 I2
(20)
N1 I2 ]
where I2 is a second-order unit matrix, and
N1 = L12 (L1 + 3L 2 + 3L3)
Discretization and finite element format
N2 = L12 (c3 L2
c2 L3) + (c2
N3 = L12 (b2 L3
b3 L 2 ) + (b3
N4 = L22 (L2 + 3L3 + 3L1)
Fig. 1 shows a schematic diagram of the planar triangular element constructed herein. The 4 nodes of the element include the three corner points (points 1, 2, 3 in Fig. 1) and the centroid of the triangle (the point 4 in Fig. 1). Each corner point of the triangular element has 6 degrees of freedom (u , ui, x , ui, y , v , vi, x , vi, y ), and the centroid point has 2 degrees of freedom (u and v ) with a total degrees of freedom of 20. In Fig. 1, x and y are the global coordinates of the triangular element, L1, L2 and L3 are the area coordinates of the triangular element. Wherein, the area coordinates are defined as follows: in the triangular element of Fig. 1, the position of the any point P can be determined according to L1, L2 and L3 , and
N5 = L 22 (c1 L3
c3 L1) + (c3
N6 = L 22 (b3 L1
b1 L3) + (b1
N7 = L32 (L3 + 3L1 + 3L2 ) N8 = L32 (c2 L1
7L1 L 2 L3 c3) L1 L 2 L3 b2 ) L1 L 2 L3 7L1 L 2 L 3 c1) L1 L2 L3 b3) L1 L2 L3 7L1 L 2 L3
c1 L2) + (c1
c2) L1 L 2 L3
N9 = L32 (b1 L2 b2 L1) + (b2 N10 = 27L1 L2 L3
b1) L1 L 2 L3 (21)
Based on the above expression, the displacement of any point in the element can be expressed as follows (22)
u = Nae
With introduction of the differential operators L , L , B and B , the relation between the strain , the strain gradient and the element’s basic unknown parameter ae can be expressed as
3
= L u = B ae = L u = B ae
(23)
Only consider the two-dimensional plane, then
B =LN B =L N
4 1
(24)
where L and L can be expressed as
P
L =
2
x
0
Fig. 1. Schematic diagram of the triangular element. 3
y
y 2
L =
0
x2
0
(25a)
x
0
2
2
x y
x y
2
2
x y
x2
0
0
2
y2
2
2
y2
x y
(25b)
Results in Physics 13 (2019) 102231
J. Wang, et al.
F
F 1
3
2
4
5
7
9
8
10
H L (a)
6 (b)
Fig. 2. The beam structure : (a) Schematic diagram; (b) Discretization model.
Based on the definition above, the characteristic length lke is expressed as a matrix I, and the equation (15) can be expressed as (ae ) =
aT
GT ne
+ aT
GT ne
Ve
Ve
parameters; (Ⅳ) gnf: it stored structural load information; (Ⅴ) gbc: it stored structural boundary constraint information. 3) The finite element model was solved by the SolveModel subprogram. Firstly, the element stiffness matrix and the element load vector were solved to integrate the structural overall stiffness matrix K and the overall load vector P. On this basis, the displacement, strain and stress distribution of the structure under a given load were further obtained. 4) The solved results were post-processed by the Megpost subprogram to obtain the displacement, strain data and the corresponding cloud image, and analyze the finite element numerical of the structure.
(B T DB + B T IDI TB + 2B T DI TB ) dV G a NTfdV +
Se
(NTt¯ +
NTr¯ T ) dS = 0
(26)
where I is the matrix form of D is the elastic tensor matrix, and ne is the number of elements. Define the element stiffness matrix as K e and the element load vector Pe , we can obtain from Eq. (26) that
lke ,
Ke =
Pe =
Ve
Ve
(B T DB + B T IDI TB + 2B T DI TB ) dV
N TfdV +
Se
G TK e G a = ne
(N Tt¯ +
N Tr¯ T ) dS
G TPe ne
(27)
Case verification
(28)
Fig. 2(a) is a schematic view of a beam fixed at both ends. The beam span-height-ratio is 4, the span length (L) is 4000 mm, the height (H) is 1000 mm, and the thickness is 200 mm. We assumed the elastic modulus of the beam material is 2.1 × 105 MPa, the Poisson's ratio is 0.3, the internal characteristic length is 0, and its midspan is subjected to a vertical concentrated force F of 20 kN. MIDAS software is a finite element analysis software for structural design, which mainly includes Midas building, Midas Civil, Midas GTS and Midas XD. In order to verify the accuracy of the program developed in this paper, we analyzed the deflection of the beam shown in Fig. 2(a) using both the program developed in this paper and the Midas GTS software. Midas GTS software [36] is equipped to handle the entire range of geotechnical design applications including deep foundations, excavations, complex tunnel systems, seepage analysis, consolidation analysis, embankment design, dynamic and slope stability analysis. The discretization model of the beam structure and the node number are shown in Fig. 2(b) (the element internal nodes are not marked). Table 1 lists the deflection value of the each node calculated by the two methods, Fig. 3 shows the comparison of the values of each node. It can be seen from Table 1 and Fig. 3 that the results from both the softwares show that the deflection of the node at the middle of the span is the largest, while the deflection of the nodes at the end of the beam is
(29)
Thus, the expressions of the structural overall stiffness matrix K and the overall load vector P are
K= P=
ne
G TK e G
ne
G TPe
Ka = P
(30) (31)
Preparation of finite element program In section 3.1, the finite element format of the triangular element based on MGE theory was derived. Based on the finite element format, the corresponding calculation program was compiled by the MATLAB software. The process of programming includes the following steps: 1) The input file (the format is dat) including the node number, node coordinates, element number, element composition, material properties, node load, and constraints was obtained by pre-processing. The pre-processing phase is performed by the GiD software which is a universal, adaptive and user-friendly pre and post processor for numerical simulations in science and engineering [35]. The preprocessing specifically includes the following four steps: (Ⅰ) building the finite element model using the GiD software; (Ⅱ) exporting the finite element analysis of model information in text form by the GiD software; (III) modifying the derived model information to acquire the input file matching by the MATLAB software; (Ⅳ) reading the input file to the MATLAB software using the FemModel subprogram. 2) We used the FemModel subprogram to define the model and read the model information file (dat file), and obtained the following variables: (Ⅰ) gNode: it stored the node number and coordinate information; (Ⅱ) gElement: it stored the element number and its composition information; (III) gMaterial: it stored material
Table 1 The result of node deflection. Node number
Proposed program Y1 (×10−4 mm)
MIDAS Y2 (×10−4 mm)
Relative error (%)
1 2 3 4 5 6 7 8 9 10
0.00 0.00 1.63 1.64 3.47 3.04 1.63 1.64 0.00 0.00
0.00 0.00 1.65 1.61 3.57 2.96 1.65 1.61 0.00 0.00
0.00 0.00 1.76 2.17 2.75 2.65 1.76 2.17 0.00 0.00
Note: The relative error is calculated as |(Y2 − Y1)/Y2|. 4
Results in Physics 13 (2019) 102231
J. Wang, et al.
y
∞
Interface
2×10-2 m
σxy x
Material 1
Material 2
Analyzing domin
∞ Fig. 4. Bi-material shear boundary layer analysis model. Table 2 Material parameters.
Fig. 3. Comparison of calculated deflection of each node.
zero due to the displacement constraints. The closer the node is to the middle of the span, the greater the deflection. The highest error of 2.75% is obtained in node 5. If the size of the finite element is small enough, the deflection curve of the beam must be smooth with the center of the span symmetrically, and the deflection at the middle of the span is the largest. Therefore, the maximum error can occur only in the midspan of the beam, namely node 5. The errors are zero in nodes 1, 2, 9 and 10 because there are displacement constraints in these nodes (the deflection value of nodes 1,2,9 and 10 are all zero). Overall, the calculation results from the program developed in this paper are close to that of the Midas GTS. Therefore, the finite element program based on the MGE theory proposed in this paper is feasible.
Material number
Elasticity modulus E (MPa)
Poisson’s ratio μ
Characteristic length l
1 2
2 × 104 2 × 103
0.3 0.3
l1x , l1y l2x , l2y
c c c c
Fig. 5. Schematic diagram of the divided triangle mesh element.
Analysis of shear boundary layer calculation
c = 1.25 × 10−4 m. Take the internal characteristic lengths of the two materials in the y direction to be 0, that is, l1y = l2y = 0, and the internal characteristic lengths in the × direction are l1x = l2x = lx = 2.0 × 10-4 m, l1x = l2x = lx = 3.0 × 10-4 m, l1x = l2x = lx = 4.0 × 10-4 m and l1x = l2x = lx = 5.0 × 10-4 m. The influence of the internal characteristic length in the × direction on the calculation results is analyzed. Fig. 6 shows the shear strain distributions of the analyzed area at different internal characteristic lengths. In the figure, the point x = 0.01 m is the interface between the two materials. It can be
The bi-material shear boundary layer is a classic case study in the strain gradient theory. In order to further explain the rationality of the finite element program proposed in this paper, it was used to analyze the bi-material shear boundary layer issue. The effect of internal characteristic length on the calculation results The MGE theory relates the macroscopic deformation characteristics to the microstructure by introducing the internal characteristic length vector of the material, so it can reflect the influence of the material microstructure on the macroscopic mechanical behavior. In order to obtain the influence of the internal characteristic length of the material on the macroscopic mechanical properties, four different internal characteristic lengths were used to analyze the bi-material shear boundary layer. A schematic diagram of the selected analysis model is shown in Fig. 4. The model is a strip-like structure in which two strain gradient materials are bonded (regardless of gravity), with material 1 on the left and material 2 on the right. In the model, the y direction is infinitely long, the left side is a fixed end, and the right side is subjected to a uniform shear stress σxy = 50 MPa. The area between the two dashed lines in Fig. 4 was selected as the analysis area (0.02 × 0.004 m2). Material parameters [30] are shown in Table 2. Discretization (meshing) of the analyzed area is performed using the triangular element proposed in Section 3.1. A schematic diagram of the individual grids divided is shown in Fig. 5. The triangular mesh element in Fig. 5 is a right-angled isosceles triangle, and the size of the triangular element refers to the side length c of the right-angled side. Here, when dividing the mesh, the size of the triangular element is
lx = 2.0×10-4 m lx = 3.0×10-4 m
lx = 4.0×10-4 m
lx = 5.0×10-4 m
x(×10-2m) Fig. 6. Distributions of Shear strain with different lx. 5
Results in Physics 13 (2019) 102231
J. Wang, et al.
Fig. 7. The shear strain curves of the points on the x-axis: (a) Shear strain with different lx; (b) Shear strain curves with different ly.
seen from Fig. 6 that as the internal characteristic length increases (from 2.0 × 10−4 m to 5.0 × 10−4 m), the boundary layer thickness of the material on the left side of the interface (Material 1) changes a little, while the boundary layer thickness of the material on the right side of the interface (Material 2) is increasing. According to the shear strain calculation results of the points on the x-axis (y = 0), the shear strain curves with different lx can be plotted, as shown in Fig. 7(a). It can be seen from Fig. 7(a) that for material 1, the variation of the internal characteristic length has a certain influence in the area near the fixed end of the left side (area I). At the same position in the area, the shear strain decreases slightly as the length of the internal characteristic increases. For material 2, the internal characteristic length has a much greater influence on the strain value in the area close of the material interface (region II), and the material shear strain at the same position in the region decreases with the increase of the internal characteristic length. From the above analysis, we can see that when the MGE theory is used to calculate and analyze the bi-material shear boundary layer, the calculation result of the material with small elastic modulus is more sensitive to the value for the internal characteristic length. Assume the internal characteristic lengths of the two materials in the x direction are both 0, that is, l1x = l2x = 0, and the internal characteristic lengths in the y direction are as follows: l1y = l2y = ly = 2.0 × 10−4 m, l1y = l2y = ly = 3.0 × 10−4 m, l1y = l2y = ly = 4.0 × 10−4 m and l1y = l2y = ly = 5.0 × 10−4 m. Fig. 7(b) shows the shear strain curves with different ly. It can be seen that the shear strain curves for different internal characteristic lengths almost coincide. Therefore, the change in the characteristic length of the y-direction has no effect on the calculation result, which is mainly due to the infinite length of the material in the y direction. In summary, when the MGE theory is used to analyze the bi-material shear boundary layer problem with infinite length in one direction, the change of the internal characteristic length of the material in the infinite direction (in the same direction as the shear stress) has little effect on the shear strain calculation results, while the variation of the internal characteristic length in the finite length direction has a certain influence on the calculation results of the shear strain, and the impact on the material with small elastic modulus in the bi-material is greater.
is analyzed with the program proposed in the paper. The analytical model, material parameters, and analysis areas are shown in Fig. 4 and Table 2, respectively (same as Section 4.1). In order to analyze the mesh dependence of the calculation results, the analyzed regions were divided into five different sizes of meshes. The mesh shape is shown in Fig. 5. The mesh sizes (right-edge length is c) are 2.00 × 103 m (2 × 10 × 2, 40 elements), 1.00 × 10−3 m (4 × 20 × 2, 160 elements), 5.00 × 10−4 m (8 × 40 × 2, 640 elements), 2.50 × 10−4 m (16 × 80 × 2, 2560 elements) and 1.25 × 10-4m (32 × 160 × 2, 10,240 elements). The finite element meshes are shown in Fig. 8. According to the analysis in the Section 4.1, the variation of the characteristic length in the y direction has no effect on the calculation results of the bi-material shear boundary layer. Therefore, the internal characteristic lengths of the two materials in the y direction are both set to 0, that is, l1y = l2y = 0, and the internal characteristic lengths in the x direction are set to l1x = 5 × 10−5 m and l2x = 5 × 10−4 m respectively. Fig. 9(a) shows the shear strain distributions of the analyzed area with different mesh sizes. In the Fig. 9(a), the point x = 0.01 m is the interface between the two materials. It can be seen from Fig. 9(a) that as the mesh size decreases, the boundary layer thickness (the width of the region where the shear strain changes in the vicinity of the interface) of the material on the left side of the interface (Material 1) continuously reduces, and the thickness of the boundary layer of the material on the right side of the interface (Material 2) also decreases. But when the mesh size reduces to 5.00 × 10−4 m, the thickness of the boundary layer on the right side of the interface is basically no longer changed. According to the shear strain calculation results of the points on the x-axis (y = 0), the shear strain curves on the corresponding section (xaxis) of the analyzed area can be plotted, as shown in Fig. 9(b). Fig. 9(b) shows: (1) For material 1, the strain at the fixed end and the vicinity of the interface between the two materials (regions I and II) varies with the mesh size, which means there is a mesh dependence. It is because that the internal characteristic length of material 1 (l1x = 5 × 10-5 m) is smaller than the five mesh sizes (cmin = 1.25 × 10−4 m). (2) For material 2, when the mesh sizes are c = 2.00 × 10-3 m and c = 1.00 × 103 m, the strain value is different from those in the other three conditions within a certain range of interface attachments (regions III and IV); when c = 5.00 × 10-4 m, c = 2.50 × 10-4 m and c = 1.25 × 10-4 m, the curves in the Fig. 9(b) basically coincide, that is, there is no mesh dependency. This is because when c = 2.00 × 10-3 m and c = 1.00 × 10-3 m, the mesh sizes are larger than the internal characteristic length of material 2 (l2x = 5 × 10-4 m); when c = 5.00 × 10-4 m, c = 2.50 × 104 m and c = 1.25 × 10-4 m, the mesh sizes do not exceed the internal characteristic length of the material 2.
Mesh dependency analysis When the classical elastic theory is used to analyze the finite element calculation of the bi- material shear boundary layer problem, the calculation results have mesh dependence [20]. The mesh dependence means that the calculation results vary greatly with the size of the finite element mesh. Here, the mesh dependence of the finite element method 6
Results in Physics 13 (2019) 102231
J. Wang, et al.
Fig. 8. Finite element meshes: (a) c = 2.00 × 10−3 m; (b) c = 1.00 × 10−3 m; (c) c = 5.00 × 10−4 m; (d) c = 2.50 × 10−4 m; (e) c = 1.25 × 10−4 m.
Based on the above analysis, when the mesh size is larger than the internal characteristic length of the material, the calculation result has mesh dependence; when the mesh size does not exceed the internal characteristic length, the calculation result is independent of the mesh size. This is consistent with the findings of Song et al. [30]. Therefore, when calculating the bi-material shear boundary layer with the finite element method, it is recommended that the mesh size to be divided should not exceed the internal characteristic length value of the material.
Shu et al. [20], the internal characteristic length of the two materials is l1 = l2 = l = 5.0 × 10-4 m. Fig. 10 shows the shear strain results for the points on the x-axis (y = 0) calculated by the three methods. It can be seen from Fig. 10 that: (1) the shear strain obtained by the three methods is continuous at the interface of the two materials; (2) The calculation results of shear strain based on the MGE triangular element are basically consistent with the calculation results of analytical solutions deduced by Shu et al. [20] in the whole analyzed region, and there is small error within a certain range (0.008≦x≦0.015 m) near the interface of the two materials; (3) The shear strain calculation results based on MGE rectangular elements are significantly different from the calculation results of analytical solutions deduced by Shu et al. [20] in the whole analyzed region. In the areas where the shear strain evenly distributes (0.0025≦x≦0.085 m and 0.015≦x≦0.02 m), the shear strain value is about 1/2 of the analytical solution. It can be seen that the calculation accuracy of the triangular element constructed in this paper is higher than that of the rectangular element. Because the theory and the size of the two finite elements are both the same, the only difference between them is the displacement mode. The displacement mode of triangular element constructed in this paper is obtained by the first-order Hermite interpolation, while the displacement mode of rectangular element is obtained by the bidirectional Hermite interpolation. Therefore, the difference in precision between the two finite elements may be caused by the displacement mode.
Comparison of calculation methods In order to illustrate the rationality and accuracy of calculation results of the MGE-based triangular finite element constructed in this paper, the MGE triangular finite element, the MGE rectangular finite element [30] and the analytical solution proposed by Shu et al. [20] based on the Toupin-Mindlin coupled stress theory were used to compare and analyze the bi-material shear boundary layer in Fig. 4. Material parameters are shown in Table 2. When meshing, the size of the triangular element is taken as c = 1.25 × 10-4 m, and the size of the rectangular element is m = n = 1.25 × 10-4 m (m and n are the length and width of the rectangular element). In the MGE theory, the internal characteristic length of the two materials in the y direction is l1y = l2y = ly = 0, and the internal characteristic length in the × direction is l1x = l2x = lx = 5.0 × 10-4 m. In the analytical solutions deduced by c=2.00×10-3 m c=1.00×10-3 m c=5.00×10-4 m c=2.50×10-4 m c=1.25×10-4 m
x (×10-2 m)
(a)
(b)
Fig. 9. The shear strain calculation results with different mesh sizes: (a) Distributions of Shear strain; (b) Shear strain curves(x-axis). 7
Results in Physics 13 (2019) 102231
J. Wang, et al.
Fig. 10. Comparison of shear strain calculation results.
Fig. 11. Effect of ly on the shear strain curves (x-axis).
Discussions
the rectangular element, the triangular element constructed in this paper has the following advantages: (1) As can be seen from the comparison in Fig. 10, the shear strain results based on the triangular element constructed in this paper are in the better agreement with analytic solutions. (2) The triangular element has better adaptability to the boundary of the analyzed material than the rectangular element, especially when the material boundary is complicated (such as arc or irregular), it has more advantage in the meshing. (3) Since the triangular element has smaller number of degrees of freedom, its stiffness matrix is of the lower order when compared with the rectangular element, and its calculation is faster. Thus, the calculation efficiency of the triangular element is higher in the case of the same algorithm. Many scholars have carried out a lot of research on the MGE theory and its finite element implementation, but there are still some shortcomings in the current research, which are mainly in: (1) The internal characteristic length is an important parameter in the strain gradient theory, but its physical meaning is still not clearly and uniformly explained. (2) It is unclear whether the internal characteristic length parameter is the intrinsic property of the material or changing with the deformation of the material. (3) The existing method for inverting the internal characteristic length parameter of the material through experiments needs high experimental conditions, which is not universally applicable. (4) At present, the research on strain gradient theory is still in theoretical and experimental research, and the theory is not well applied to engineering practice. Therefore, the current research on the strain gradient theory and its programmatic is still not perfect, and there is still a lot of research work to be done in this area.
For the bi-material shear boundary layer problem, the classical elastic theory holds that the stress continuously and evenly distributes at the interface of the two materials, and the strain is abrupt, which is obviously inconsistent with the actual condition [19–21,37,38]. Therefore, the finite element analysis of the bi-material shear boundary layer problem based on the classical elasticity theory inevitably has mesh dependence. The MGE theory relates the macroscopic deformation of the material to its internal microstructure by introducing the internal characteristic length vector, so it can reflect the influence of the microstructure of the material on the mechanical behavior. The calculation results in this paper (Figs. 7 and 9) show that the shear strain at the interface of the two materials calculated by the MGE theory is continuous. Therefore, compared with the solution given by the classical elastic theory, the calculation results of the bi-material shear boundary layer given by the MGE theory are obviously more in line with the actual situation. In addition, Shu et al. [20] gave an analytical solution of the shear strain of the bi-material shear boundary layer based on the ToupinMindlin coupled stress theory. But since the internal characteristic length in the Toupin-Mindlin coupled stress theory is scalar, it makes the analytical solution fail to reflect the influence of the anisotropy of the microstructure inside the material on the mechanical behavior. For the MGE theory, although the shear calculation results (Fig. 7(b)) of the bi-material model with infinite length in the y direction in Fig. 4 show that the variation of the internal characteristic length in the y direction has no effect on the calculation results of the shear strain, the size of the material in one direction cannot be infinitely long. Take the length of the two materials in the y direction of Fig. 4 as 0.004 m (finite length), and the rest of the conditions are unchanged, the material parameters are shown in Table 2; take the internal characteristic lengths of the two materials in the × direction are both 5.0 × 10-4 m (l1x = l2x = 5.0 × 10-4 m), and the internal characteristic lengths of the two materials in the y direction are also the same (l1y = l2y = ly). Fig. 11 shows the effect of variation of the internal characteristic length in the y direction on shear strain calculation results of the points on the x axis. It can be seen that the variation of the internal characteristic length in the y direction has a certain influence on the shear strain calculation result of the bimaterial shear boundary layer, and the influence on the material with small elastic modulus (material 2) is more obvious. Since the internal characteristic length of the material in the MGE theory is vector, this theory can reflect the influence of the anisotropy of the microstructure inside the material on the mechanical behavior. The triangular element and the rectangular element [30] based on the MGE theory are both high-order finite elements, but compared with
Conclusions 1) In this paper, a triangular finite element is constructed based on the MGE theory, and its finite element format is derived. We also compiled the corresponding finite element program by the MATLAB software. The program was used to analyze the bi-material shear boundary layer problem with unidirectional infinite length. 2) The variation of the internal characteristic length of the material in the infinite length direction has little effect on the shear strain calculation results, and the variation of the internal characteristic length in the finite length direction has a certain influence on the shear strain calculation result, and the impact on the material with small elastic modulus in the bi-material is greater. 3) When the finite element mesh size is larger than the internal characteristic length of the material, the calculation result has mesh dependence, and the mesh dependence can be eliminated when the mesh size is smaller than the characteristic length. 4) The shear strain calculated by the MGE triangular element is in good agreement with the analytical solutions based on the Toupin8
Results in Physics 13 (2019) 102231
J. Wang, et al.
Mindlin coupled stress theory, and the calculation result of the MGE rectangular element is only about 1/2 of the analytical solutions in the uniformly distributed area of the shear strain, indicating that the MGE triangle element has higher calculation accuracy.
[18] Askes H, Aifantis EC. Gradient elasticity in statics and dynamics: an overview of formulations, length scale identification procedures, finite element implementations and new results. Int J Solids Struct 2011;48:1962–90. https://doi.org/10. 1016/j.ijsolstr.2011.03.006. [19] Fleck NA, Hutchinson JW. A phenomenological theory for strain gradient effects in plasticity. J Mech Phys Solids 2001;41:1825–57. https://doi.org/10.1016/00225096(93)90072-n. [20] Shu JY, King WE, Fleck NA. Finite elements for materials with strain gradient effects. Int J Numer Methods Eng 2015;44:373–91. https://doi.org/10.1002/(sici) 1097-0207(19990130)44:3<373::aid-nme508>3.0.co;2-7. [21] Zhang ZH, Zhuang Z, Gao Y, Liu ZL, Nie JF. Cyclic plastic behavior analysis based on the micromorphic mixed hardening plasticity model. Comput Mater Sci 2011;50:1136–44. https://doi.org/10.1016/j.commatsci.2010.11.013. [22] Yang F, Chong ACM, Lam DCC, Tong P. Couple stress based strain gradient theory for elasticity. Int J Solids Struct 2002;39:2731–43. https://doi.org/10.1016/S00207683(02)00152-X. [23] Mohammad-Abadi M, Daneshmehr AR. Modified couple stress theory applied to dynamic analysis of composite laminated beams by considering different beam theories. Int J Eng Sci 2015;87:83–102. https://doi.org/10.1016/j.ijengsci.2014. 11.003. [24] Lim CW, Zhang G, Reddy JN. A higher-order nonlocal elasticity and strain gradient theory and its applications in wave propagation. J Mech Phys Solids 2015;78:298–313. https://doi.org/10.1016/j.jmps.2015.02.001. [25] Saroukhani S, Vafadari R, Simone A. A simplified implementation of a gradientenhanced damage model with transient length scale effects. Comput Mech 2013;51:899–909. https://doi.org/10.1007/s00466-012-0769-8. [26] Askes H, Gutiérrez MA. Implicit gradient elasticity. Int J Numer Methods Eng 2010;67:400–16. https://doi.org/10.1002/nme.1640. [27] Soh AK, Wanji C. Finite element formulations of strain gradient theory for microstructures and the C0–1 patch test. Int J Numer Methods Eng 2004;61:433–54. https://doi.org/10.1002/nme.1075. [28] Chen Z. Assumed strain finite element method based on the theory of strain gradient. J Univ Sci Technol Beijing 2007;29:807–10. https://doi.org/10.1016/S10016058(07)60030-4. [29] Zhao J, Chen WJ, Bin JI. Formula derivation in axisymmetric strain gradient theory and finite element implementation. Chin J Computat Mechan 2009;26:209–14. https://doi.org/10.1109/CLEOE-EQEC.2009.5194697. [30] Song Z, Bing Z, He J, Zheng Y. Modified gradient elasticity and its finite element method for shear boundary layer analysis. Mech Res Commun 2014;62:146–54. https://doi.org/10.1016/j.mechrescom.2014.09.008. [31] Wunderlich W, Pilkey WD. Mechanics of structures. Variational and computational methods. Meccanica 2004;39:291–2. https://doi.org/10.1023/B:MECC. 0000023038.64148.bc. [32] Li X, Zhang J, Zheng Y. Static and free vibration analysis of laminated composite plates using isogeometric approach based on the third order shear deformation theory. 232019-232019 Adv Mech Eng 2015;6. https://doi.org/10.1155/2014/ 232019. [33] Zienkiewicz OC, Nakazawa S. On variational formulation and its modifications for numerical solution. Comput Struct 1984;19:303–13. https://doi.org/10.1016/ 0045-7949(84)90231-1. [34] Meleshkina AV. On the approximation of the derivatives of the Hermite interpolation polynomial on a triangle. Zh Vychisl Mat Mat Fiz 2010;50:201–10. https:// doi.org/10.1134/S0965542510020016. [35] García-Doñoro D, García-Castillo LE, Gómez-Revuelto I. An interface between an – adaptive finite element package and the pre- and post-processor GiD. Finite Elem Anal Des 2010;46:328–38. https://doi.org/10.1016/j.finel.2009.11.005. [36] Wang J, Huo Q, Song Z, Zhang Y. Study on adaptability of primary support arch cover method for large-span embedded tunnels in the upper-soft lower-hard stratum. Adv Mech Eng 2019;11:1–15. https://doi.org/10.1177/ 1687814018825375. [37] Liu X, Han Y, Li D, Tu Y, Deng Z, Yu C, et al. Anti-pull mechanisms and weak interlayer parameter sensitivity analysis of tunnel-type anchorages in soft rock with underlying weak interlayers. Eng Geol 2019;253:123–36. https://doi.org/10.1016/ j.enggeo.2019.03.012. [38] Lai J, Qiu J, Fan H, Zhang Q, Hu Z, Wang J, et al. Fiber bragg grating sensors-based in situ monitoring and safety assessment of loess tunnel. J Sens 2016;2016:1–10. https://doi.org/10.1155/2016/8658290.
Acknowledgments The authors gratefully acknowledge the support received from National Natural Science Foundation of China (Grant Nos. 51578447, 51404184), China, the Science and Technology Project of Ministry of Housing Urban-Rural Construction (Grant No. 2017-K4-032), China, the Youth Science and Technology Nova Program of Shaanxi Province (Grant No. 2018KJXX-061), China and Natural Science Basic Research Program of Shaanxi Province (Grant No. 2016JQ4009), China. References [1] Huang Y, Zhang F, Hwang KC, Nix WD, Pharr GM, Feng G. A model of size effects in nano-indentation. J Mech Phys Solids 2016;54:1668–86. 0.1016/ j.jmps.2006.02.002. [2] Xu J, Dai CY, Zhang B, Zhang GP. Strain-gradient dependent fatigue behavior of micron-thick copper single crystal foils. Comput Mater Sci 2014;85:223–9. https:// doi.org/10.1016/j.commatsci.2014.01.004. [3] Sahmani S, Aghdam MM. Nonlinear primary resonance of micro/nano-beams made of nanoporous biomaterials incorporating nonlocality and strain gradient size dependency. Results Phys 2018;8. DOI:10.1016/j.rinp.2018.01.002. [4] Wang J, Song Z, Zhao B, Liu X, Liu J, Lai J. A study on the mechanical behavior and statistical damage constitutive model of sandstone. Arab J Sci Eng 2018;43:5179–92. https://doi.org/10.1007/s13369-017-3016-y. [5] Hoover CG, Bažant ZP. Cohesive crack, size effect, crack band and work-of-fracture models compared to comprehensive concrete fracture tests. Int J Fract 2014;187:133–43. https://doi.org/10.1007/s10704-013-9926-0. [6] Zervos A. Some finite elements for elasticity with microstructure and gradient elasticity. Int J Numer Methods Eng 2010;73:564–95. https://doi.org/10.1002/ nme.2093. [7] Lai JX, Mao S, Qiu JL, Fan HB, Zhang Q, Hu ZN, et al. Investigation progresses and applications of fractional derivative model in geotechnical engineering. Math Probl Eng 2016;3:1–15. https://doi.org/10.1155/2016/9183296. [8] Qiu J, Qin Y, Lai J, Wang K, Niu F, Wang H, et al. Structural response of the metro tunnel under local dynamic water environment in loess strata. Geofluids 2019;2019:16. https://doi.org/10.1155/2019/8541959. [9] Wang J, Ren Z, Song Z, Huo R, Yang T. Study of the effect of micro-pore characteristics and saturation degree on the longitudinal wave velocity of sandstone. Arab J Geosci 2019. in presss. [10] Cosserat E, Cosserat F. Théorie des corps déformables. Nature 1909;81:67. https:// doi.org/10.1038/081067a0. [11] Borg U. A strain gradient crystal plasticity analysis of grain size effects in polycrystals. Eur J Mech A-Solids 2007;26:313–24. https://doi.org/10.1016/j. euromechsol.2006.09.006. [12] Farajpour A, Rastgoo A. Influence of carbon nanotubes on the buckling of microtubule bundles in viscoelastic cytoplasm using nonlocal strain gradient theory. Results Phys 2017;7:1367–75. https://doi.org/10.1016/j.rinp.2017.03.038. [13] Lazopoulos KA, Lazopoulos AK. Fractional derivatives and strain gradient elasticity. Acta Mech 2016;227:823–35. https://doi.org/10.1007/s00707-015-1489-x. [14] Tsagrakis I, Konstantinidis A, Aifantis EC. Strain gradient and wavelet interpretation of size effects in yield and strength. Mech Mater 2003;35:733–45. https://doi. org/10.1016/s0167-6636(02)00205-3. [15] Mindlin RD. Micro-structure in linear elasticity. Arch Ration Mech Anal 1964;16:51–78. https://doi.org/10.1007/bf00248490. [16] Konstantopoulos I, Aifantis E. Gradient elasticity applied to a crack. J Mech Behav Biomed Mater 2013;22:193–201. https://doi.org/10.1515/jmbm-2013-0026. [17] Aifantis EC. Update on a class of gradient theories. Mech Mater 2003;35:259–80. https://doi.org/10.1016/s0167-6636(02)00278-8.
9