Geometrical Nonlinearity for a Timoshenko Beam with Flexoelectricity

The Timoshenko beam model is applied to the analysis of the flexoelectric effect for a cantilever beam under large deformations. The geometric nonlinearity with von Kármán strains is considered. The nonlinear system of ordinary differential equations (ODE) for beam deflection and rotation are derived. Moreover, this nonlinear system is linearized for each load increment, where it is solved iteratively. For the vanishing flexoelectric coefficient, the governing equations lead to the classical Timoshenko beam model. Furthermore, the influence of the flexoelectricity coefficient and the microstructural length-scale parameter on the beam deflection and the induced electric intensity is investigated.


Introduction
The size effect on micro/nano structural elements is observed due to the comparable sizes of these elements and the material microstructural length-scale. Classical continuum models include only the first-order gradients of primary fields and do not reflect the material microstructure. Therefore, these models are inapplicable of describing and explaining the phenomena (e.g., flexoelectricity), where the second-order gradients play a role. On the other hand, higher-grade continuum theories contain higher gradients of primary fields. In addition, the additional coefficients in constitutive relationships are determined by the material microstructure, which is usually accomplished through the microstructural length-scale parameter. Buckling of centrosymmetric anisotropic beams was analyzed in Reference [1] with the strain-gradient theory. Then, the geometric nonlinearity with von Kármán strains [2] are taken into account. Wang et al. [3] analyzed the nonlinear free vibration of electrically actuated nanobeams, where the surface energy, temperature change, geometrical nonlinear deformation, and intermolecular Casimir force are considered.
Thin beams are commonly used in nanoelectromechanical devices and for energy harvesting [4]. Traditionally, the piezoelectricity is utilized to convert mechanical energy to electrical energy or vice versa [5]. In piezoelectric materials, a uniform mechanical strain can induce an electric polarization. This conversion is observed only in non-centrosymmetric crystal structures. Numerous crystalline materials are not piezoelectric since their structure is centrosymmetric. However, a non-uniform strain or the presence of strain-gradients may potentially break the inversion symmetry and induce polarization, even in centrosymmetric crystals [6][7][8]. In the literature, this phenomenon is called the flexoelectric effect [9,10]. Flexoelectric energy harvesters are intensively studied in [11][12][13]. A rectangular beam is the most common model used for this purpose. The Euler-Bernoulli and Timoshenko beam theories are applied to analyze the beams. In the Euler-Bernoulli theory, there are vanishing shear strain deformations. The Timoshenko beam model based on the strain-gradient elasticity theory and the couple stress elasticity theory [14,15] has been developed only recently [16]. The microstructure-dependent Timoshenko beam model contains a material length-scale parameter and can capture the size effect. A flexoelectric Euler-Bernoulli model for energy harvesting is proposed in [17]. The first attempt to consider the geometric nonlinearity deformation in flexoelectricity has been developed by Wang and Wang [18] with a beam described by the Euler-Bernoulli theory. The authors used simplification, where the gradient of normal strain along the beam vanishes. This leads to the simplification of the nonlinear system of differential equations. An exact solution is proposed for the nonlinear forced vibration analysis of nanobeams made of functionally graded materials (FGMs) subjected to a thermal environment, including the effect of surface stress in [19]. Thai et al. [20] proposed an isogeometric approach for general problems of flexoelectricity in soft dielectric materials at finite deformations, accounting for Maxwell stresses. A nonlinear vibration of flexoelectric nanobeams with surface and thermal effects is investigated by the surface elasticity theory [21]. A similar nonlinear bending problem with Kármán nonlinear deformations for the coupled piezomagnetic-flexomagnetic nanosized Euler-Bernoulli beam is analyzed in [22]. Recently, they extended an early paper to the geometrically nonlinear vibration of the piezo-flexomagnetic nanotube [23]. Free vibrations of a visco-piezo-flexoelectric nanobeam are given in [24]. Sahmani and Aghdam [25] applied the nonlocal strain-gradient beam model with the third-order distribution of shear deformation to explore the nonlinear vibration of axially-loaded multilayer functionally graded nanobeams in both the pre-buckling and post-buckling domains. A similar nanoplate problem for the electro-mechanical shear buckling analysis by the modified couple stress theory with various boundary conditions is given in [26].
In the present paper, von Kármán large deformations in the direct flexoelectricity are considered for a cantilever beam without simplification. The Timoshenko model with geometrical nonlinearity is applied to derive the nonlinear system of ordinary differential equations (ODEs) for the beam deflection and rotation. The mechanical load in the nonlinear system is considered incrementally. In each increment, the system of ODEs should be linearized, where nonlinear terms are considered iteratively.

A Linear Theory of Direct Flexoelectricity
The electric enthalpy density for piezoelectric solids with the direct flexoelectricity can be written as [8,27]: where the symbols a and c are used for the second-order permittivity and the fourth-order elastic constant tensors, respectively. The piezoelectric coefficient is denoted by symbol e, and symbol f is the direct flexoelectric coefficient. The tensor g is used for higher-order elastic coefficients representing the strain-gradient elasticity. The strain tensor ε ij and the electric field vector E j are defined as [28]: where u i and φ are the displacements and electric potential, respectively. The strain-gradient tensor η ijk is given by: Under the infinitesimal deformations, the constitutive equations can be obtained from the electric enthalpy density expression (1) [27,29]: where σ ij , D k , and τ jkl are the stress tensor, electric displacements, and higher-order stress tensor, respectively.
The size scale of higher-order elastic parameters g jklmni is expressed by a proportionality of the conventional elastic stiffness coefficients c klmn and the microstructural length-scale parameter l [1,30]: with δ li as the Kronecker delta. Deng et al. [31] considered two independent components f 1 and f 2 for the direct Then, the electric enthalpy density has the following form: If the poling direction is along the x 3 -axis in the piezoelectric material, the electric enthalpy has the following form: The constitutive Equation (4) for orthotropic materials, can be rewritten into a matrix form as in [32]:  The governing equations for the piezoelectric solid with direct flexoelectric effects are given in [33]: Moreover, one can find the form of essential and natural boundary conditions (b.c.): 1. Essential b.c.: 2.
Natural b.c.: where the traction vector and the electric charge are defined as follows: with ρ i := n k π j τ ijk , δ(x) is the Dirac delta function, and π i is the Cartesian component of the unit tangent vector on Γ.
The jump at a corner (x c ) on the oriented boundary contour Γ is defined as: Regarding the electric boundary conditions, the prescription of the electric potential (as given by (12)) also includes the short-circuit case, when the electric potential is constant on the whole boundary (vanishing). On the other hand, in the case of open-circuit conditions, the electric potential is unknown and a certain value of the free electric charge on the boundary is prescribed. If the case of the short-circuit condition is applied to the beam problem considered in the next section, the electric field E 3 is equal to zero. Therefore, the open-circuit conditions are applied.

Timoshenko Beam with Flexoelectric Effect and Nonlinear Strains
Next, we analyze a cantilever beam by taking into account the direct flexoelectric effect and large von Kármán strains (geometric nonlinear deformation). The x 1 -axis is the neutral axis of the undeformed beam, and the x 3 -axis is along the thickness direction. A uniform transverse load q is applied on the upper surface of the cantilever beam ( Figure 1).
Regarding the electric boundary conditions, the prescription of the electric potentia (as given by (12)) also includes the short-circuit case, when the electric potential is constan on the whole boundary (vanishing). On the other hand, in the case of open-circuit condi tions, the electric potential is unknown and a certain value of the free electric charge on the boundary is prescribed. If the case of the short-circuit condition is applied to the beam problem considered in the next section, the electric field E3 is equal to zero. Therefore, the open-circuit conditions are applied.

Timoshenko Beam with Flexoelectric Effect and Nonlinear Strains
Next, we analyze a cantilever beam by taking into account the direct flexoelectri effect and large von Kármán strains (geometric nonlinear deformation). The 1 x -axis is the neutral axis of the undeformed beam, and the 3 x -axis is along the thickness direction. A uniform transverse load q is applied on the upper surface of the cantilever beam ( Figure  1). In the traditional Timoshenko beam theory with finite shear stresses, the displace ment of the beam can be expressed as:  In the traditional Timoshenko beam theory with finite shear stresses, the displacement of the beam can be expressed as: where w(x 1 ) is the transverse displacement of the neutral axis and φ(x 1 ) is the rotation of the cross-section. With the assumption of very small slopes in the beam after defor-mation but a possible finite transverse deflection, the von Kármán nonlinear strain of the beam [1,2] is: The only non-zero strain-gradients are given as: Wang and Wang [18] simplified their solution neglecting η 111 = ε 11,1 , as compared to η 113 . All of the three components of strain gradients in (20) are considered in this work.
For the beam, we assume that only the E 3 electric intensity component is nonvanishing. In the case of open-circuit condition, the electric displacement on the surface is vanishing. Therefore, D 3 = 0 on x 3 = ±H/2, and from the Maxwell equation for the electric displacement D 3,3 = 0, we obtain D 3 = 0 inside the beam. Then, from the constitutive Equation (9), we directly obtain the expression for the electric intensity component E 3 .
In the case of the short-circuit condition, the electric intensity component E 3 is equal to zero. Then, the stress components are dependent only on the strains. In this case, finite values of higher-order stress tensors are found in (4), only due to the non-vanishing higherorder elastic parameters g jklmni . In addition, the governing Equation (11) for mechanical quantities (beam deflection) is dependent on the electrical fields induced by flexoelectricity. If g jklmni are vanishing, the governing equation for the beam deflection is similar to the classical elasticity without electro-mechanical coupling [34].
Next, the open-circuit conditions are considered. Substituting the electric intensity component E 3 (21) into the stress (8) and higher-order stress tensors (10), we obtain: The variational principle for this special case with a vanishing electric displacement has the following form: thus with the use of strain (19) and strain-gradient (20) expressions, we obtain: where W is the work of the external transverse load q(x) and L is the beam length.
Defining the bending moment, shear and normal forces, as well as their higher-order counterparts, we obtain: Applying the integration in parts to the above equation, we obtain: Since variations δφ and δw can be arbitrary in the interval [0, L], we obtain two governing equations: Furthermore, the possible boundary conditions are obtained from Equation (27).
(i) Total shear force or the beam deflection vanishes at the end of the beam: (ii) Total bending moment or rotation of the cross-section vanishes at the end of the beam: (iii) Higher-order shear force or the deflection slope vanishes at the end of the beam: (iv) Higher-order bending moment or gradient of rotation vanishes at the end of the beam: Substituting Equations (19), (20), (22) and (23) into (25), we obtain: Iφ ,1 = −Sφ ,1 , Next, notations are used for the second moment of the cross-section area: Substituting (30) into (28), a nonlinear system of ordinary differential equations is obtained: where Λ = Z + S + Ψ. For the beam subjected to the uniform transverse load and clamped at the end x 1 = 0, while for the free end at x 1 = L, the following boundary conditions are prescribed: In order to solve the above nonlinear boundary value problem, we employ the weak formulation: Using the direct iteration technique (Picard iteration method), the solution at the k-th iteration is determined from the linearized equation: where the coefficient matrix is evaluated using the known solution from the (k − 1)-st iteration. The initial guess vector {U} (0) is taken as the solution of equation where K is the coefficient matrix obtained from K({U}) by neglecting the nonlinear terms. The iterative procedure stopped when the difference between the two consecutive iterations, measured by the Euclidean norm, is less than the tolerance ε:
The length of the beam is L = 500 nm, the width is B = H, and the thickness is H = 20 nm. The microstructural length-scale parameter is chosen as l = 1 × 10 −8 m.
The weak form of the linearized differential Equation (34) has been implemented into the commercial software Comsol. The cubic Hermitian shape function has been used as interpolation functions and equidistantly distributed 1D finite elements.
The computational procedure for the 1D Euler-Bernoulli model is first verified for a linear model with the vanishing von Kármán strains. The convergence of the numerical results for deflection with respect to the increasing number of discretization elements is achieved for N > 50 finite elements, as shown in Table 1. In the following sections, we will use the finite elements with a length of 10 nm, when the convergence is guaranteed. 10 elements 6.353723 × 10 −9 20 elements 6.750830 × 10 −9 50 elements 6.750842 × 10 −9 100 elements 6.750842 × 10 −9 Moreover, the variation of the beam deflection along x 1 for various flexoelectric coefficients is obtained with the 2D analysis by the mixed FEM program [36], where 2000 rectangular elements are used. An opposite verification of the new computer code for the 2D problems by the Euler-Bernoulli theory was applied in [36] for a body force load, where the analytical solution is available.
The solid line in Figure 2 corresponds to the 1D Euler-Bernoulli model and the symbols are valid for the 2D FEM analysis. Here, one can observe the excellent agreement between the present numerical and analytical results. The induced electric field represented by the electric intensity E 3 on the neutral beam axis is presented in Figure 3. In addition, one can observe that the induced electric intensity increases with the increasing flexoelectric coefficients. It is in agreement with the results recently presented by the authors in [36]. The largest values of the electric intensity are at the clamped end, where the bending moment and the strain-gradients reach their maximum values. The vanishing electric potential is prescribed at the free end of the cantilever beam.    [37][38][39]. Similar to the linear elastic case, the increasing flexoelectricity reduces the beam deflection in the nonlinear case, as well. The variation of the electric intensity 3 E on the neutral beam axis is presented in Figure 5, when the load In Figure 4, one can see that at this load intensity, the nonlinearity is significantly large. From a comparison with Figure 5 (nonlinear case) and Figure  3 (linear case), one can see a significantly lower quantity of the induced electric intensity in the nonlinear case. bending moment and the strain-gradients reach their maximum values. The vanishing electric potential is prescribed at the free end of the cantilever beam.   [37][38][39]. Similar to the linear elastic case, the increasing flexoelectricity reduces the beam deflection in the nonlinear case, as well. The variation of the electric intensity 3 E on the neutral beam axis is presented in Figure 5, when the load In Figure 4, one can see that at this load intensity, the nonlinearity is significantly large. From a comparison with Figure 5 (nonlinear case) and Figure  3 (linear case), one can see a significantly lower quantity of the induced electric intensity in the nonlinear case. Next, the effect of nonlinearity can be investigated. The iteration procedure is applied to consider the finite values of the von Kármán strains. For the iteration procedure, we have chosen the tolerance ε = 0.005, which has been achieved after five iterations. The variation of the beam maximal deflection w with the load q for various flexoelectric coefficients and at the microstructural length-scale parameter l = 1 × 10 −8 m is shown in Figure 4. The dashed line represents the deflection corresponding to the linear model. Here, one can observe that the decrease of nonlinear deflections in comparison with the linear ones is enhanced with the increasing load. It is in agreement with the observations in the classical theory of elasticity [37][38][39]. Similar to the linear elastic case, the increasing flexoelectricity reduces the beam deflection in the nonlinear case, as well. The variation of the electric intensity E 3 on the neutral beam axis is presented in Figure 5, when the load intensity is q = 0.02 N/m. In Figure 4, one can see that at this load intensity, the nonlinearity is significantly large. From a comparison with Figure 5 (nonlinear case) and Figure 3 (linear case), one can see a significantly lower quantity of the induced electric intensity in the nonlinear case.
Next, the influence of the microstructural length-scale parameter l on the beam deflection is investigated. If the value of this parameter grows by 50% with respect to the previous case in Figure 4, the beam deflection is further reduced as shown in Figures 4 and 6. Evidently, the effect of nonlinearity is weaker if the intrinsic material parameter is larger.  Next, the influence of the microstructural length-scale parameter l on the beam deflection is investigated. If the value of this parameter grows by 50% with respect to the previous case in Figure 4, the beam deflection is further reduced as shown in Figures 4  and 6. Evidently, the effect of nonlinearity is weaker if the intrinsic material parameter is larger.  Next, the influence of the microstructural length-scale parameter l on the beam deflection is investigated. If the value of this parameter grows by 50% with respect to the previous case in Figure 4, the beam deflection is further reduced as shown in Figures 4  and 6. Evidently, the effect of nonlinearity is weaker if the intrinsic material parameter is larger.  The variation of maximal deflection w with the intrinsic parameter l for two different loading levels and at the flexoelectric coefficient f 1 = 5 × 10 −8 C/m is shown in Figure 7. Both of the curves for the two different load intensities are nearly parallel.
The influence of the flexoelectric coefficient on the maximal beam deflection is clearly visible in Figures 8 and 9 for two different values of the microstructural length-scale parameter. Here, one can observe that the deflection response is lower for the beam with the larger value of the microstructural length-scale parameter. The beam deflection decreases with the increasing value of the flexoelectric coefficient.  The influence of the flexoelectric coefficient on the maximal beam deflection is clearly visible in Figures 8 and 9 for two different values of the microstructural length-scale parameter. Here, one can observe that the deflection response is lower for the beam with the larger value of the microstructural length-scale parameter. The beam deflection decreases with the increasing value of the flexoelectric coefficient.

Conclusions
In nano-sized structures, higher-order derivatives of strains require consideration due to the large gradients of deformations. The mechanical and electrical response of the

Conclusions
In nano-sized structures, higher-order derivatives of strains require consideration due to the large gradients of deformations. The mechanical and electrical response of the cantilever beam on mechanical loading is investigated within the higher-grade electroe-

Conclusions
In nano-sized structures, higher-order derivatives of strains require consideration due to the large gradients of deformations. The mechanical and electrical response of the cantilever beam on mechanical loading is investigated within the higher-grade electroelasticity, including the direct flexoelectricity and von Kármán nonlinear strain deformations. In addition, the bending of the beam is studied within the Timoshenko beam theory. The principle of virtual work is applied to derive a nonlinear system of ordinary differential equations (ODEs) for the coupled fields of the electric potential, beam deflection, and rotation. The FEM scheme is employed for discretization and approximation of spatial field variations with the resulting system of nonlinear algebraic equations, as solved by the iteration procedure.
From the numerical results obtained for the clamped beam, one can observe a reduction of the mechanical as well as electrical response of the beam within the nonlinear model, as compared with the linear one. Furthermore, the response decreases with the increasing level of the mechanical load. The flexoelectric effect consumes a portion of the energy of the external forces and effectively hinders the bending of the beam under a pure mechanical load. The intensity of the induced electric field rapidly grows with the increasing value of the flexoelectric coefficient.