Acta Polytechnica CTU Proceedings doi:10.14311/APP.2017.7.0033 Acta Polytechnica CTU Proceedings 7:33–37, 2017 © Czech Technical University in Prague, 2017 available online at http://ojs.cvut.cz/ojs/index.php/app 1D HERMITE ELEMENTS FOR C1-CONTINUOUS SOLUTIONS IN SECOND GRADIENT ELASTICITY Christian Liebold∗, Wolfgang H. Müller Berlin University of Technology, Chair of Continuum Mechanics and Materials Theory, Einsteinufer 5, 10587 Berlin, Germany ∗ corresponding author: christian.liebold@tu-berlin.de Abstract. We present a modified strain gradient theory of elasticity for linear isotropic materials in order to account for the so-called size effect. Additional material length scale parameters are introduced and the problem of static beam bending is analyzed. A numerical solution is derived by means of a finite element approach. A global C1-continuous displacement field is applied in finite element solutions because the higher-order strain energy density additionally depends on second gradients of displacements. So-called Hermite finite elements are used that allow for merging gradients between elements. The element stiffness matrix as well as the global stiffness matrix of the problem is developed. Convergence, C1-continuity and the size effect in the numerical solution is shown. Experiments on bending stiffnesses of different sized micro beams made of the polymer SU-8 are performed by using an atomic force microscope and the results are compared to the numerical solution. Keywords: second gradient continuum, size effect, Hermite finite elements. 1. Introduction Materials with intrinsic micro- or nano-structure can show size-dependent material behavior, which is re- flected, e.g., by a stiffer elastic response to external forces when the size of the material body is reduced. A quantitative understanding of a possible size effect in micro- and submicrostructures is of great importance when modeling Micro- and Nano-Electro-Mechanical Systems (MEMS/NEMS). Driven by miniaturization, the need for saving materials and improving the per- formance, the requirement of reliability in simulation technologies increases. A size effect is reflected in a different elastic response to external loads if the size of the material body is reduced. Experimental validation therefore is given in, e.g., [1–6]. Second Gradient (SG) continuum theories are used if, in the absence of strain gradients (in uniaxial tensile tests), the elastic behavior of the material is independent of the thickness of the sample [3]. They are applicable to so-called non-simple materials of the gradient type, for example, polymers at a small scale. Since con- ventional continuum theories based on the Cauchy continuum are not able to predict size effects, the present work deals with the Modified Strain Gradient theory (MSG) developed in, e.g., [7–9]. By employ- ing the macroscopic rotation vector and second order derivatives of the displacement vector, MSG-theory is able to account for quantities like rotation and curva- ture. A generalization of strain gradient continua has been developed by Eringen, who proposed “nonsim- ple materials of gradient type [10]” in order to derive the corresponding higher-order material dependencies in a rational manner. The application of conventional Finite Element (FE) strategies may lead to inaccurate results, if finite element formulations are used which only fullfill global C0–continuity. The scope of this work is to use FE formulations based on Hermite polynomials in order to fulfill global C1–continuity for solutions in the modified strain gradient theory. Section 2 presents the strain energy function of the modified strain gradient theory and a one dimensional differential equation is derived by making use of the Euler-Bernoulli assumptions for slender beams. The finite element formulation based on Hermite polynomials is presented in Sect. 3 as well as the structure of the element stiffness- and system matrix. A numerical example for a loaded cantilever beam is then given in Sect. 4. 2. Modified strain gradient theory 2.1. Strain energy function The present study starts with one of the three re- duced forms of the strain energy densities for small deformations, uSG, as postulated by Mindlin (1962) [8]. The Einstein summation convention on repeated indices is used and spatial partial derivatives in the Cartesian coordinate system are denoted by comma- separated indices. Mindlin’s second form of linear isotropic strain energy density reads: uSG = α1εij εij + α2εkk εmm + β1ηijk ηijk + β2ηiik ηjjk + + β3ηiik ηkjj + β4ηijj ηikk + β5ηijk ηkji , (1) where εij = 1/2(ui,j + uj,i ) denotes the small strain tensor and ηijk = 1/2(uk,ij + uj,ki ) = εkj,i the gradient 33 http://dx.doi.org/10.14311/APP.2017.7.0033 http://ojs.cvut.cz/ojs/index.php/app Christian Liebold, Wolfgang H. Müller Acta Polytechnica CTU Proceedings of strain. α1 and α2 lead to Lamé’s constants and β1,...,5 are five additional higher-order material con- stants [11]. σij and µijk denote the Cauchy stress tensor and the higher-order stress tensor, respectively: σij = ∂uSG ∂εij , µijk = ∂uSG ∂ηijk . (2) The inclusion of the macroscopic rotation vector ϕi = 1/2�ijk uk,j (�ijk being the Levi-Civita symbol) leads to a reduction of the additional material constants from five to three. Based on Fleck & Hutchinson (1997) [12] we introduce independent expressions of ηijk and decompose the second order displacement gradient into a symmetric part ηSijk and a remaining part ηRijk (which is not necessarily anti-symmetric): Hj Hj Hj Z~ Hj �� �� �� �� ηijk ηSijk η R ijk η̄ij χSijχ A ij η (0) ijk η (1) ijkχ A ijεmm,i Figure 1. Scheme of decomposition ηijk = ηSijk + η R ijk , η S ijk = 1 3 (uk,ij + ui,jk + uj,ki ) , ηRijk = 2 3 (�ikl η̄lj + �jkl η̄li ) + �kjl η̄li , (3) where η̄ij = ϕj,i is the gradient of rotation, which we decompose into its symmetric and anti-symmetric parts, χSij and χAij , respectively: χAij = 1 2 (ϕi,j −ϕj,i ) , χ S ij = 1 2 (ϕi,j + ϕj,i ) . (4) The tensor ηSijk is further decomposed into its spheri- cal and deviatoric part, η(0)ijk and η (1) ijk , cf., Fig. 1 . The quantity η(0)ijk is related to χ A ij and the dilatation gra- dient εmm,i in the following manner: η (0) ijk = 1 5 ( δij η S mmk + δjk η S mmi + δki η S mmj ) ηSmmi = εmm,i + 2 3�iln χ A ln , η (1) ijk = η S ijk −η (0) ijk . (5) As shown in [13, 14] χAij does not influence the strain energy if symmetry of the couple stress tensor µij is assumed. As a consequence the linear isotropic strain energy density for non-simple materials of gradient type reads: uMSG = 2Gεij εij + λεkk εii + 2G`20εmm,i εkk,i + + 2G`21η (1) ijk η (1) ijk + 2G` 2 2χ S ij χ S ij , (6) where λ and G are Lamé’s constants, whereas `0,`1 and `2 denote additional material length scale parame- ters. These parameters are chosen to be of the dimen- sion of a length and squared in order to guarantee a positive definite problem in the energy minimization. The normalization of the higher-order terms by G is arbitrary. 2.2. Governing Euler-Bernoulli differential equation By using the displacement field of an Euler- Bernoulli beam of length L, thickness t, cross- section A, Young’s modulus E, and second moment of inertia I cf., Fig. 2, ux = −z dw(x) dx , uy = 0 , uz = w(x), (7) L, EI, A w(x) dw(x) dx F xy z Figure 2. Bending line of a cantilever where w(x) characterizes the bending line and F de- notes the force at the free end, Eqs. (3)–(6) lead to the strain energy density of the body. By applying the principle of virtual work to the variation of the func- tion w(x), a triple integration by parts leads to the Ordinary Differential Equation (ODE) of the problem [15]: S d4w dx4 −K d6w dx6 = 0 , ∀ x ∈ [0,L] . (8) The resulting ODE is of rank six, where S = EI(1 + 4308225 `2 t2 ) and K = 75EI` 2 (9) are given constants, if a rectangular cross-section and `0 = `1 = `2 = ` is assumed. 3. Finite element implementation 3.1. Weak form and discretization A multiplication of Eq. (8) with a test function v(x) and triple integration by parts yields in the following weak integral form of the ODE: L∫ 0 ( S d2w dx2 d2v dx2 + K d3w dx3 d3v dx3 ) dx = 0 . (10) Now the Galerkin method is applied to discretize the interval [0,L] by using linear combinations of indepen- dent base functions φα for each element e = 1, . . . ,N, which are uniformly defined on the element’s coordi- nate 0 ≤ ζ ≤ 1: ve(ζ) = ∞∑ β=1 φβ , w e(ζ) = ∞∑ δ=1 ceδφδ . (11) The general element function of this method is derived by inserting Eq. (11) in Eq. (10), evaluated on the domain ζ: 1∫ 0 ∞∑ δ=1 ∞∑ β=1 ( Sceδ d2φδ dζ2 d2φβ dζ2 + Kceδ d3φδ dζ3 d3φβ dζ3 ) dζ = 0 . (12) 34 vol. 7/2017 C1-Continuous Solutions in Second Gradient Elasticity 3.2. Hermite elements for global c1–continuity In order to solve the present weak form of the problem, the base functions need to be differentiable three times at least. In order to guarantee that information of second gradients of the trial function w(x) is included, we ask for global continuity of the first derivative of w(x), which is called C1–continuity of the solution. Therefore, we choose a linear combination of four different Hermite polynomials as follows: H1 = 2ζ3 − 3ζ2 + 1 , H2 = ζ3 − 2ζ2 + ζ H3 = −2ζ3 + 3ζ2 , H4 = ζ3 − ζ2, (13) H1 H2 H3 H4 0 ≤ ζ ≤ 1 Figure 3. Four Hermite polynoms per element which allow to influence the first derivative at each boundary of the element independently. The unit ele- ment now is a linear combination of the four Hermite polynomials: ve(ζ) = 4∑ β=1 Hβ , w e(ζ) = 4∑ δ=1 ceδHδ (14) where ceδ are four constants per element. 3.3. Structure of the element stiffness- and system matrix By using Hermite elements as base functions Eq. (12) can be written in the following manner: 1∫ 0 4∑ δ=1 4∑ β=1 ( Sceδ d2Hδ dζ2 d2Hβ dζ2 +Kceδ d3Hδ dζ3 d3Hβ dζ3 ) dζ = 0. (15) The element stiffness matrix Ke reads: Keδβ = 1∫ 0 ( S d2Hδ dζ2 d2Hβ dζ2 + K d3Hδ dζ3 d3Hβ dζ3 ) dζ =   Ke11 K e 12 K e 13 K e 14 Ke22 K e 23 K e 24 Sym. Ke33 Ke34 Ke44   = const. (16) and is a symmetric 4×4 matrix with constant compo- nents that depend on the material data, beam geome- try, and on the length of the element. In this approach we choose equidistantly distributed elements. By spe- cial assembling of the system matrix it is now possible to connect the boundary values of an element to its neighbors and, in addition, to bring together the first derivatives in between. This will be achieved by the following structure of the system matrix K: K =   K1δβ + K2δβ + K3δβ + K···δβ + KNδβ   . (17) Due to the summation of the last entries, the specific components of neighboring elements are treated as equivalent. This refers to both, the components that represent direct node values as well as to the corre- sponding derivatives. The solution of the entire beam bending problem is now approximated by a linear system of equations: Kc = f ⇒ K\f = c , (18) where c is the vector of element coefficients and f is a right hand side, in which the external force F is imple- mented (here, situated at the very last node N). The system is solved using the ”backslash” operator of the commercial program MATLAB (LU decomposition, etc.). 4. Numerical example As an example we study bending of a cantilever beam, fixed at the left hand side and loaded with a single force at the free end, as presented in Fig. 2 . Material and geometry data is given in Tab. 1 . The material length scale parameter ` is chosen according to a literature value [3] and to measurements performed by the authors with an atomic force microscope [16, 17]. t W L E `0 = `1 = `2 F 100 µm 2t 20t 1.44 GPa 17.6 µm 1 µN Table 1. Dimensions, loads and material parameters The boundary condition for the fixation has been implemented by canceling the first and second row and column of the system matrix Eq. (17), since, regarding the first node of the first element, the current values of deflection and the derivative of the deflection line are equal to zero. 35 Christian Liebold, Wolfgang H. Müller Acta Polytechnica CTU Proceedings er r ×10−5 0 5 10 15 20 0 40 80 120 160 number of elements N Figure 4. Convergence The test for convergence with respect to higher num- bers of elements N yield positive results (cf. Fig. 4). The relative error of tip deflections err= wN/w1 −1 per increase of elements achieves values below 0.001 h. LKMBerlin University of TechnologyFaculty V – Mechanical Engineering and Transport Systems Institute of Mechanics Chair of Continuum Mechanics and Materials Theory Prof. Dr. rer. nat. W. H. Müller Table: Dimensions, loads and material parameters 𝑡 𝑊 𝐿 𝐸 ℓ0 = ℓ1 = ℓ2 𝐹 100 µm 2𝑡 20𝑡 1.44 GPa 17.6 µm 1 µN 𝑤 [m m ] 𝑥/𝐿 d𝑤 d𝑥 𝑥/𝐿 d2 𝑤 d𝑥2 𝑥/𝐿 d3 𝑤 d𝑥3 𝑥/𝐿 Dr.-Ing. C. Liebold | ExNum 2016, Liblice | 1D example 9/10 Figure 5. Deflection line The deflection line of the beam solution is presented in Fig. 5 . Five elements are shown so that a closer look at the solution of each element becomes possi- ble, which consists of reconstructed trial functions, cf. Eq. (14)2. Figure 6 presents a plot of the first derivative of the reconstructed trial functions per ele- ment and clearly demonstrates continuity in between neighboring elements and therefore C1–continuity. LKMBerlin University of TechnologyFaculty V – Mechanical Engineering and Transport Systems Institute of Mechanics Chair of Continuum Mechanics and Materials Theory Prof. Dr. rer. nat. W. H. Müller Table: Dimensions, loads and material parameters 𝑡 𝑊 𝐿 𝐸 ℓ0 = ℓ1 = ℓ2 𝐹 100 µm 2𝑡 20𝑡 1.44 GPa 17.6 µm 1 µN 𝑤 [m m ] 𝑥/𝐿 d𝑤 d𝑥 𝑥/𝐿 d2 𝑤 d𝑥2 𝑥/𝐿 d3 𝑤 d𝑥3 𝑥/𝐿 Dr.-Ing. C. Liebold | ExNum 2016, Liblice | 1D example 9/10Figure 6. First derivative of the solution As expected the global solution is not C2–continuous. This is confirmed by the jumps between the elements in their second derivative in Fig. 7. Such jumps are also present in the third derivative in Fig. 8, but invisibly small. LKMBerlin University of TechnologyFaculty V – Mechanical Engineering and Transport Systems Institute of Mechanics Chair of Continuum Mechanics and Materials Theory Prof. Dr. rer. nat. W. H. Müller Table: Dimensions, loads and material parameters 𝑡 𝑊 𝐿 𝐸 ℓ0 = ℓ1 = ℓ2 𝐹 100 µm 2𝑡 20𝑡 1.44 GPa 17.6 µm 1 µN 𝑤 [m m ] 𝑥/𝐿 d𝑤 d𝑥 𝑥/𝐿 d2 𝑤 d𝑥2 𝑥/𝐿 d3 𝑤 d𝑥3 𝑥/𝐿 Dr.-Ing. C. Liebold | ExNum 2016, Liblice | 1D example 9/10 Figure 7. Second derivative of the solution LKMBerlin University of TechnologyFaculty V – Mechanical Engineering and Transport Systems Institute of Mechanics Chair of Continuum Mechanics and Materials Theory Prof. Dr. rer. nat. W. H. Müller Table: Dimensions, loads and material parameters 𝑡 𝑊 𝐿 𝐸 ℓ0 = ℓ1 = ℓ2 𝐹 100 µm 2𝑡 20𝑡 1.44 GPa 17.6 µm 1 µN 𝑤 [m m ] 𝑥/𝐿 d𝑤 d𝑥 𝑥/𝐿 d2 𝑤 d𝑥2 𝑥/𝐿 d3 𝑤 d𝑥3 𝑥/𝐿 Dr.-Ing. C. Liebold | ExNum 2016, Liblice | 1D example 9/10Figure 8. Third derivative of the solution E ∗ [G P a] thickness t [mm] Figure 9. Size effect for decreasing thicknesses The sensitivity of the implemented numerical model of the modified second gradient continuum was tested for incorporation of size effect by decreasing the thick- ness values, while all other dimensions were scaled in parallel according to Tab. 1 . The solution of the tip deflection was used in a calculation of an equivalent Young’s modulus, E∗, with a Euler-Bernoulli analytical beam model of the conventional Cauchy continuum. As plotted in Fig. 9, a size effect is clearly present. For small thicknesses the equivalent Young’s modulus is calculated to be about three times larger than the usual one. 36 vol. 7/2017 C1-Continuous Solutions in Second Gradient Elasticity 5. Conclusions Based on a modified strain gradient theory for lin- ear isotropic materials we have evaluated the Euler- Bernoulli beam assumptions to derive the ordi- nary differential equation of the problem. The weak form was derived using test functions. Third order Hermite polynomials have been summed up for an element formulation: These are able to access and influence first derivatives between the elements. The element stiffness and system matrix have been for- mulated for the problem of a cantilever beam. The resulting system of linear equations was solved nu- merically to determine the coefficients of the trial functions for each element. The global solution was reconstructed element-wise and first, second and third order derivatives have been calculated. It was demon- strated, that this procedure guarantees C1–continuity of the solution. This is considered to be an important quality in second gradient theories. Furthermore, the size effect, measured by the authors in previous work, could be reproduced by using similar material as well as geometry data. This provides the impetus for the development of a three dimensional finite element for- mulation of a strain gradient continuum on the basis of Hermite elements. References [1] W. J. Poole, M. F. Ashby, N. A. Fleck. Micro-hardness of annealed and work-hardened copper polycrystals. Scripta Materialia 34(4):559–564, 1996. [2] X. H. Guo, D. N. Fang, X. D. Li. Measurement of deformation of pure ni foils by speckle pattern interferometry. Mechanical Engineering 27(2):21–25, 2005. [3] D. C. C. Lam, F. Yang, C. M. Chong, et al. Experiments and theory in strain gradient elasticity. J Mech Phys Sol 51(8):1477–1508, 2003. [4] A. W. McFarland, J. S. Colton. Role of material microstructure in plate stiffness with relevance to microcantilever sensors. Journal of Micromechanics and Microengineering 15(5):1060–1067, 2005. [5] S. Cuenot, C. Fretigny, S. Demoustier-Champagne, B. Nysten. Surface tension effect on the mechanical properties of nanomaterials measured by atomic force microscopy. Physical Review B 69:01–05, 2004. [6] X.-F. Li, B.-L. Wang, K. Y. Lee. Size effect in the mechanical response of nanobeams. J of Adv Research in Mech Eng 1(1):4–16, 2010. [7] R. Toupin. Elastic materials with couple-stresses. ARMA 11:385–414, 1962. [8] R. Mindlin. Effects of couple-stresses in linear elasticity. ARMA 11:415–448, 1962. [9] R. D. Mindlin, N. N. Eshel. On first strain-gradient theories in linear elasticity. International Journal of Solids and Structures 4:109–124, 1968. [10] A. Eringen. Nonlocal continuum field theories. Springer-Verlag, New York 2010. [11] M. Lazar. Irreducible decomposition of strain gradient tensor in isotropic strain gradient elasticity. ZAMM Z Angew Math Mech 2016. doi:10.1002/zamm.201500278. [12] N. A. Fleck, J. W. Hutchinson. Strain gradient plasticity. Advances in Applied Mechanics, Academic Press, New York 33:295–361, 1997. [13] C. Liebold, W. H. Müller. Measuring material coefficients of higher gradient elasticity by using afm techniques and raman-spectroscopy. Advanced Structured Materials, Springer Berlin Heidelberg 22:255–271, 2013. [14] F. Yang, C. M. Chong, D. C. C. Lam, P. Tong. Couple stress based strain gradient theory for elasticity. International Journal of Solids and Structures 39(10):2731–2743, 2010. [15] S. Kong, S. Zhou, Z. Nie, K. Wang. Static and dynamic analysis of micro beams based on strain gradient elasticity theory. International Journal of Engineering Science 47:487–498, 2009. [16] C. Liebold, W. H. Müller. Applications of strain gradient theories to the size effect in submicro- structures incl. experimental analysis of elastic material parameters. Bulletin of TICMI 19(1):45–55, 2015. [17] C. Liebold, W. H. Müller. Comparison of gradient elasticity models for the bending of micromaterials. Elsevier, Computational Materials Science 116:52–61, 2016. 37 http://dx.doi.org/10.1002/zamm.201500278 Acta Polytechnica CTU Proceedings 7:33–37, 2017 1 Introduction 2 Modified strain gradient theory 2.1 Strain energy function 2.2 Governing Euler-Bernoulli differential equation 3 Finite element implementation 3.1 Weak form and discretization 3.2 Hermite elements for global c1–continuity 3.3 Structure of the element stiffness- and system matrix 4 Numerical example 5 Conclusions References