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

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 beammodel.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.Themoving-least square approximation is used to construct the shape function and its secondand 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.


Introduction
The size effect has been found in the areas of micromechanics and nanomechanics, for example thin films, microelectromechanical systems, and nanoelectromechanical systems [1][2][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][5][6][7][8][9][10][11][12][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][10][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 sizedependent 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 requires C 1 -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][17][18][19][20][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 rotationfree 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.

Higher-Order Bernoulli-Euler Beam Model
As shown in Figure 1, a beam is considered, in which the xaxis 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 in which  denotes the displacement of the elastic line, that is, deflection.
The double strains are defined as [3] The classical stress becomes and the double stresses are given by where  is the elastic Young's modulus, and g and   are the size effect parameters which are the intrinsic bulk length and the directional surface length, respectively.The principle of virtual work requires in which  ex is the work of the external forces,  in is the strain energy of the beam, and  is the total potential.Moreover, the variation of the strain energy of the beam is defined by [3]  in which  denotes the length of beam, and  denotes the area of the cross-section.
In addition, the variation of the work of the external forces is given by where  denotes the continuously distributed loading,   the shear forces at the ends of the beam, and  and  the bending moments and double moments of the beam.
Substituting ( 2)-( 4) into ( 6), and calculating the inner integral, we can get where I is the moment of inertia of the cross-section with respect to the y-axis.

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 components  , and  , 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 deflect () defined on the segment [0, ] is set as where   () is the monomial basis function,  is the number of terms in the basis, and   () are the coefficients of the basis function.To construct the third-order derivative, the following cubic basis function is used: The unknown coefficients   () in ( 9) can be determined by the minimization of the weighted discrete  norm where (−  ) is the weight function with compact support, and NP is the number of nodes with ( −   ) > 0. The minimum of J in (11) with respect to a() leads to a set of linear equations: where w = ( 1 ,  2 , . . .,  NP )  , and Thus, the unknown coefficients a() can be obtained from (12) as Substituting ( 14) into ( 9), we can get where the MLS shape function   () is defined as If we set the partial derivatives of   () with respect with  can be obtained as in which with In the present work, the weight function is chosen as where  is the support size for the node, and () is the cubic spline function [18]: (22)

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, ) is divided into a series of subdomains that are called integral cells in some papers [18][19][20][21].
For the calculation convenience, the integral cells are chosen to coincide with the nodal arrangements; that is, the  − 1 integral cells are used if the  nodes are allocated.
The second and third derivatives of the deflection w can be approximated as Their institution into (8) with  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 where U is the total nodal component vector, and F 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].

Numerical Results
The efficiency of the proposed method is first tested.A thin beam of length L is supported at two points A and B with an applied vertical F at its midpoint (Figure 2).The dimensions are fixed asthe height ℎ = 1, and thickness  = 0.01.If both   and  are chosen as zero, the deflection  at the midpoint is  3 /48 according to the classical beam theory.The support size of nodes C can be expressed as  =  max ( +1 −   ) with  max 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  max 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  max has a large effect on the numerical results when the cubic basis is used, and a good choice for  max 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   and   are not necessary since   and g are fixed as zero in these tests.
Next, the size effect is studied.The total 61 nodes are used, and  max is chosen as 3.0.The cross -section is fixed as  = ℎ, and the boundary conditions are same with the above tests.  is first fixed as 0.00125b to test the effect of  on the midpoint deflection.Figure 3 plots the relatively increasing amount (%) in comparison with the results of those from the classical beam theory ( 3 /48).An increasing effect on the results can be found with the increasing g.To test the effect of   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 l x being in the certain range.

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.

Figure 1 :
Figure 1: A beam model and related markers.

Figure 2 :Figure 3 :
Figure 2: A simply supported beam and its deformation under a vertical force F.

Figure 4 :
Figure 4: The relatively increasing amount (%) of the midpoint deflection for the varying value of l x in comparison with the classical beam theory.

Table 1 :
The midpoint deflection (× 3 /) for different numbers of nodes ( max is fixed as 3.0).