# The Application of Mesh-Free Method in the Numerical Simulation of Beams with the Size Effect

Authors: Yuan Tian; Caili Zhang; Yuzhou Sun
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

The mesh-free method is employed to implement the numerical simulation of the bending behavior of beams with the size effect. On the basis of the classical Bernoulli-Euler beam theory, two higher-order strain components are involved in the beam model. The intrinsic bulk length and the directional surface length components are introduced into the constitutive relationship to describe the size effect, and the variation of the total potential is provided. The moving-least square approximation is used to construct the shape function and its second- and third-order derivatives, and the choice of the scaling factor is discussed in detail. A mesh-free scheme is built to implement numerical simulation, in which the higher-order strains are directly approximated with the nodal components due to the higher-order continuity of the shape function. The convergence of method is illustrated in virtue of an example of the simply supported beam, and the effect of the intrinsic bulk length and the directional surface length components are studied.


---

## Body

## 1. Introduction

The size effect has been found in the areas of micromechanics and nanomechanics, for example thin films, microelectromechanical systems, and nanoelectromechanical systems [1–3]. The scale effect can be successfully modeled by employing the higher-order continuum theory, in which the constitutive equations introduce some material length scale parameters, in addition to the classical material parameters [4–13]. Mostly the generally known higher-order theory includes the nonlocal theory of Eringen [14], strain gradient elasticity, and the couple stress theory [15]. The beam model is very efficient for studying static and dynamic behaviors of slender nanostructures such as carbon nanotubes and microtubules. Zhang et al. studied the small length effect on a nonlocal cantilever beam and carried out a detailed analysis of bending, buckling, and vibration of nanosized beams [4]. Lim gives detailed statements for static response and bending moments of nanosized beams via the nonlocal theory [5, 6]. Gao and Lei study the small scale effects on mechanical behaviors of protein microtubules based on the nonlocal elasticity theory [7]. The classical strain gradient elasticity theory and its modified version [8] have been applied for analysis of bending, vibration, and buckling of micro-sized beams under various boundary conditions by some researchers [9–11]. A microstructure-dependent model for the Timoshenko beam has been developed by using a modified couple stress theory and Hamilton’s principle to account for the microstructural effect [12]. Kahrobaiyan et al. has also investigated the size-dependent dynamic characteristics of atomic force microscopic cantilevers based on the modified couple stress theory [13].However, in all above studies, the microscale structural information is not incorporated, and the scale parameters need to be assumed to study their effect on the nanostructures’ behavior. On the other hand, most of studies are mainly on the theoretical analysis for some simple support and boundary conditions, and the numerical simulation is very rare due to the complexity from the introduction of the higher-order strain. Specifically, the inclusion of the higher-order gradient largely complicates the solving of the boundary value problems. For example, in the finite element method, the interpolation requiresC1-continuity [16, 17], and this leads to difficulty in the establishment of elements and the construction of the interpolation functions. Recently, researchers [17] have applied the mesh-free method [16–21] to simulate the materials with strain-gradients effects. The mesh-free method is a newly developed technique that has some distinct advantages. In particular, the moving-least square (MLS) approximations possess nonlocal properties and automatically satisfy the higher-order continuity requirement. This intrinsic nonlocal property leads to real rotation-free approximation [17], and displacements can thus be used as the only nodal freedoms. Motivated by the above fact, the present work attempts to employ the mesh-free method to build the numerical scheme for beams with size effect. To produce the continuum third-order derivatives of shape function, the cubic basis is needed, and the choice of scaling factor is thus needed to be studied. After the convergence of the established method is tested, the effect of the buck length and surface length components is studied and discussed.


## 2. Higher-Order Bernoulli-Euler Beam Model

As shown in Figure1, a beam is considered, in which the x-axis is the axis of the beam, and the y-axis is the deflection axis. The elastic line lies on the x-y plane. Following the works of Akgöz and Civalek [1] and Civalek and Demir [2], we treat the strain energy density of the beam as a function of the classical strain and the double strains. In Bernoulli-Euler principle, the classical strain is defined as
(1)εxx=-y∂2w∂x2
in which w denotes the displacement of the elastic line, that is, deflection.Figure 1
A beam model and related markers.The double strains are defined as [3]
(2)εxxx=-y∂3w∂x3,εyxx=-∂2w∂x2.
The classical stress becomes
(3)σxx=Eεxx+lxεxxx,
and the double stresses are given by
(4)τxxx=g2Eεxxx+lxEεxx,τyxx=g2Eεyxx,
where E is the elastic Young’s modulus, and g and lx are the size effect parameters which are the intrinsic bulk length and the directional surface length, respectively.The principle of virtual work requires(5)δV=δUin-δUex,
in which Uex is the work of the external forces, Uin is the strain energy of the beam, and V is the total potential. Moreover, the variation of the strain energy of the beam is defined by [3]
(6)δUin=∫0L∫A(σxxδεxx+τxxxδεxxx+τyxxδεyxx)dAdx,
in which L denotes the length of beam, and A denotes the area of the cross-section.In addition, the variation of the work of the external forces is given by(7)δUex=-∫0Lqδwdx-[Vnδw]0L+[Mδw,x]0L+[mδw,xx]0L,
where q denotes the continuously distributed loading, Vn the shear forces at the ends of the beam, and M and m the bending moments and double moments of the beam.Substituting (2)–(4) into (6), and calculating the inner integral, we can get
(8)δUin=∫0L{[E(I+g2A)∂2w∂x2+lxEI∂3w∂x3]δw,xxhhhhhh+[EIg2∂3w∂x3+lxEI∂2w∂x2]δw,xxx}dx,
where I is the moment of inertia of the cross-section with respect to the y-axis.


## 3. The MLS Shape Functions and Their Derivatives

In the present research, the MLS approximation is used to build the numerical discretization scheme. Since the MLS approximation satisfies the higher-order continuity automatically, the componentsw,xx and w,xxx in (8) can directly be approximated with the nodal components. The calculation of the MLS shape functions and their derivatives are illustrated as follows [18].The deflectw(x) defined on the segment [0,L] is set as
(9)wh(x)=∑i=1mpi(x)ai(x)=pT(x)a(x),
where pi(x) is the monomial basis function, m is the number of terms in the basis, and ai(x) are the coefficients of the basis function. To construct the third-order derivative, the following cubic basis function is used:
(10)pT=(1,x,x2,x3).
The unknown coefficients ai(x) in (9) can be determined by the minimization of the weighted discrete J norm
(11)J=∑I=1NPW(x-xI)[PT(xI)a(xI)-wI]2,
where W(x-xI) is the weight function with compact support, and  NP is the number of nodes with W(x-xI)>0. The minimum of J in (11) with respect to a(x) leads to a set of linear equations:
(12)A(x)a(x)=B(x)w,
where w=(w1,w2,…,wNP)T, and
(13)A(x)=∑I=1NPW(x-xI)p(xI)pT(xI),B(x)=[W(x-x1)p(x1),W(x-x2)p(x2),hhh…,W(x-xNP)p(xNP)].
Thus, the unknown coefficients a(x) can be obtained from (12) as
(14)a(x)=A-1(x)B(x)w.
Substituting (14) into (9), we can get
(15)wh(x)=PT(x)A-1(x)B(x)w=∑I=1NPϕI(x)wI(x),
where the MLS shape function ϕI(x) is defined as
(16)ϕI(x)=∑j=1mpj(x)[A-1(x)B(x)]jI.
If we set
(17)H(x)=pT(x)A-1(x)=∑j=1mpj(x)A-1(x),
the partial derivatives of ϕI(x) with respect with x can be obtained as(18a)ϕI,x(x)=[H,x(x)B(x)+H(x)B,x(x)]I,(18b)ϕI,xx(x)=[H,xx(x)B(x)+2H,x(x)B,x(x)+H(x)B,xx(x)]I,(18c)ϕI,xxx(x)=[H,xxx(x)B(x)+3H,x(x)B,xx(x)hhh+3H,xx(x)B,x(x)+H(x)B,xxx(x)]I,in which(19a)B,x(x)=[W1,x(x)P(x1),W2,x(x)P(x2),hhhh…,WN,x(x)P(xN)](19b)B,xx(x)=[W1,xx(x)P(x1),W2,xx(x)P(x2),hhh…,WN,xx(x)P(xN)](19c)B,xxx(x)=[W1,xxx(x)P(x1),W2,xxx(x)P(x2),hhh…,WN,xxx(x)P(xN)]with(20a)H,x(x)=p,xT(x)A-1(x)-A-1(x)A,x(x)H,(20b)H,xx(x)=P,xxT(x)A-1(x)-A-1(x)×(A,xx(x)H(x)+2A,x(x)H,x),(20c)H,xxx(x)=P,xxxT(x)A-1(x)-A-1(x)×[A,xxx(x)H(x)hhhhh+3A,xx(x)H,x(x)+3A,x(x)H,xx].In the present work, the weight function is chosen as(21)W(x-xI)=W(|x-xI|C),
where C is the support size for the node, and W(d) is the cubic spline function [18]:
(22)W(d)={23-4d2+4d3,d≤1243+4d+4d2-43d3,12<d≤10,d>1.


## 4. The Mesh-Free Computational Scheme

In our mesh-free simulation, some nodes are allocated uniformly along the beam axis. To implement the numerical integrals, the interval(0,L) is divided into a series of subdomains that are called integral cells in some papers [18–21]. For the calculation convenience, the integral cells are chosen to coincide with the nodal arrangements; that is, the N-1 integral cells are used if the N nodes are allocated.The second and third derivatives of the deflectionw can be approximated as
(23)w,xx=∑I=1NPϕI,xx(x)wI(x),w,xxx=∑I=1NPϕI,xxx(x)wI(x).
Their institution into (8) gives the stiffness matrix
(24)K=∫0L{ϕ,xxTE(I+g2A)ϕ,xx+ϕ,xxxTlxEIϕ,xxhhhhh+ϕ,xxxTEIg2ϕ,xxx+ϕ,xxTlxEIϕ,xxxφ,xxT}dx,
in which
(25)ϕ,xx=(ϕ1,xx,ϕ2,xx,…,ϕN,xx)T,ϕ,xxx=(ϕ1,xxx,ϕ2,xxx,…,ϕN,xxx)T,
with N is the total number of nodes.The Gaussian quadrature rule is used for each cell. In the practical calculation, the Gaussian points can be globally marked. To avoid the repeating calculation for each Gaussian point, the information that includes the nodes whose supports cover the Gaussian points, the corresponding shape functions, and their derivatives is calculated and deposited in advance.The solution can be found by solving the equation(26)KU=F,
where U is the total nodal component vector, andF is the nodal force vector. It is noted that the unknown vector does not involve the rotation angle; the essential boundary condition can be imposed with the penalty function method [18, 20, 21].


## 5. Numerical Results

The efficiency of the proposed method is first tested. A thin beam of lengthL is supported at two points A andB with an applied vertical F at its midpoint (Figure 2). The dimensions are fixed asthe height h=1, and thickness b=0.01. If both lx and g are chosen as zero, the deflection δ at the midpoint is FL3/48EI according to the classical beam theory. The support size of nodes C can be expressed as C=Dmax⁡(xj+1-xj) with Dmax⁡ denoting the scaling factor. Three point Gauss integral formula is used for each integral cell. Table 1 shows the midpoint deflection while the different numbers of nodes are used in the case of Dmax⁡ being fixed as 3.0. It can be seen that a stable convergence is achieved, and very good result can be obtained when the node number is larger than 15 (the exact solution is 0.25). Table 2 gives the results while the varying scaling factor is used if the node number is fixed as 61, which shows that Dmax⁡ has a large effect on the numerical results when the cubic basis is used, and a good choice for Dmax⁡ should be 3.0. For comparison, the results in the case of quadratic basis are also provided in Table 2. It should be pointed out that εxxx and εyxx are not necessary since lx andg are fixed as zero in these tests.Table 1
The midpoint deflectionδ(×FL3/EI) for different numbers of nodes (Dmax⁡ is fixed as 3.0).
Node number
5
9
13
17
21
25
29
33
61
δ
(
×
F
L
3
/
E
I
)
0.2716
0.2545
0.2516
0.2508
0.2505
0.2503
0.2502
0.2501
0.2500Table 2
The midpoint deflectionδ(×FL3/EI) for varying scaling factor (61 nodes are used).
D
max
⁡
2.0
2.5
3.0
3.5
4.0
4.5
5.0
Cubic basis
0.0564
0.0634
0.2500
0.2502
0.2817
0.2605
0.3544
Quadratic basis
0.2499
0.2501
0.2500
0.2502
0.3382
0.2601
0.2165Figure 2
A simply supported beam and its deformation under a vertical forceF.Next, the size effect is studied. The total 61 nodes are used, andDmax⁡ is chosen as 3.0. The cross -section is fixed as b=h, and the boundary conditions are same with the above tests. lx is first fixed as 0.00125b to test the effect of g on the midpoint deflection. Figure 3 plots the relatively increasing amount (%) in comparison with the results of those from the classical beam theory (FL3/48EI). An increasing effect on the results can be found with the increasingg.To test the effect of lx on the deflection,g is fixed as 0.1875b. Figure 4 plots the result comparison with the classical beam theory, which shows that an obvious effect is only associated with the value of lxbeing in the certain range.Figure 3
The relatively increasing amount (%) of the midpoint deflection for the varying value ofg in comparison with the classical beam theory.Figure 4
The relatively increasing amount (%) of the midpoint deflection for the varying value oflx in comparison with the classical beam theory.


## 6. Conclusions

The higher-order continuum theory is used to study the bending behavior of beams, which can reflect the size effect associated with some special structures. The derived formulas provide a foundation for the theoretical analysis and numerical simulation, and the current work is concerned mainly with the numerical simulation. A mesh-free scheme has been built to implement the numerical simulation. Specifically, the second- and third-order derivatives of deflection are directly approximated with the nodal components, which largely make the discretization easy. It is proven that the MLS approximation and mesh-free method are very efficient for the problems in which the higher-order strains need to be involved.The current research employs the cubic basis to produce the continuum third-order derivatives of shape function, and it is found that the better convergence can be obtained, the stricter the scaling factor needs to be used in comparison with the quadratic basis. For the beams with size effect, the size effect becomes larger with the increasing buck length component if the surface length component is fixed. But, only the value of surface length component in certain range has a large effect if the buck length component is fixed.



---

*Source: 590271-2014-03-30.xml*