An improved shear deformation theory for bending beams with symmetrically varying mechanical properties in the depth direction

The paper is devoted to simply supported beams under three-point bending. Their mechanical properties symmetrically vary in the depth direction. The individual shear deformation theory for beams of such features is proposed. Based on the principle of stationary total potential energy the differential equations of equilibrium are obtained. The system of the equations is analytically solved, and the shear coefficients and deflections of example beams are calculated. The solution is compared with other analytical results obtained with the use of another deformation function. Moreover, the bending problem of these beams is also numerically studied using the finite element method. Results of analytical and numerical studies are presented in Figures and Tables.


Introduction
The shearing effect occurring in bent constructions was noticed already in the nineteenth century, and studied in detail for homogeneous and layered constructions in the twentieth century. Assumption of a proper theory of deformation of the straight line normal to the neutral surface makes a basis for analytical modeling of heterogeneous structures, especially those with mechanical properties varying in the wall thickness direction.
Reddy [1] developed a theoretical model of bending of functionally graded rectangular plates considering the shearing effect. Detailed analysis is made taking into account the first and third order shear deformation theory. Zenkour [2] presented a generalized shear deformation theory and its application to the analysis of functionally graded rectangular plates subjected to uniformly distributed load. The transverse shear effect is studied in detail. Aydogdu [3] proposed a new shear deformation theory for laminated composite plates. This theory exactly meets the conditions for zeroing shear stresses on the upper and lower surface of the plate. Reddy [4] presented a reformulation of the classical and shear deformation beam and plate theories taking into account the nonlocal differential constitutive relations of Eringen and the von Kármán nonlinear strains. The equilibrium equations of the nonlocal beam theories and the classical and first-order shear deformation theories of plates are formulated. Carrera et al. [5] described in detail classic and advanced theories, including: the basics of the theory of deformable bodies, Euler-Bernoulli and Timoshenko theories of beams, non-linear theories, e.g. the parabolic, cubic, quartic, and n-order beam theories, as well as modeling of beams made of functionally graded materials. Meiche et al. [6] presented a new hyperbolic shear deformation theory on the example of buckling and free vibration analysis of thick functionally graded sandwich plates. This theory is more perfect in relation to the simple shear deformation theories of Mindlin and Reissner. Moreover, it provides parabolic variation of the transverse shear stresses across the thickness, and also their zeroing on external surfaces. Thai and Vo [7] developed various higher-order shear deformation theories for testing bending and free vibration of functionally graded beams. These theories account for higher-order variation of transverse shear strain in the depth direction of the beam, and satisfy the stress-free boundary conditions on the upper and lower surfaces of the beam. Thai and Vo [8] developed a new sinusoidal shear deformation theory for functionally graded rectangular plates. This theory describes the sinusoidal distribution of the transverse shear stress and meets the conditions for zeroing shear stress on the outer surfaces of the plate. Detailed tests concerning bending, buckling, and vibration of these plates have been performed.
Akgöz and Civalek [9] presented a new higher-order shear deformation analytical beam model with consideration of the strain gradient elasticity theory. This model describes the microstructural and shear deformation effects without the need of shear correction factors. The problems of static bending and free vibration of simply supported microbeams are investigated. Grover et al. [10] proposed a new inverse hyperbolic shear deformation theory of laminated composite and sandwich plates. This theory is formulated based on the shear strain shape function and validated by numerical studies of the bending and buckling problem of rectangular plates. Sahoo and Singh [11] proposed a new inverse trigonometric zig-zag theory for laminated composite and sandwich plates. This theory ensures the continuity conditions at the layer interfaces and zeroing shear stress on the outer surfaces of the plate. The effective finite element model is developed for numerical studies of static problems of these plates. Xiang [12] improved the n-order shear deformation theory with consideration of the condition for zeroing shear stress on the outer surfaces of the functionally graded beam. The free vibration problems of this beam are analyzed. Kumar and Chakraverty [13] proposed four new inverse trigonometric shear deformation theories allowing to study the free vibration of isotropic thick rectangular plates. The theories ensure meeting the transverse stress boundary conditions on both plate surfaces. A test of convergence and validation was carried out with the cases from the available literature. Mahi et al. [14] presented a new hyperbolic shear deformation theory describing bending and free vibration of isotropic, functionally graded, sandwich and laminated composite plates. The approach does not require a shear correction factor. Based on Hamilton's principle the energy functional of the system was obtained. The accuracy of the method was shown by comparisons with the numerical solution of the problem.
Darijani and Shahdadi [15] proposed a new deformation plate theory with consideration of the shear deformations. The transverse shear stresses vary across the plate thickness according to a power-law relationship. The upper and bottom surfaces of the plate are free of shear stress. The governing equations and boundary conditions of the plate are derived with the use of Hamilton's principle. The results are comparable to those obtained using higher-order theories. Lezgy-Nazargah [16] considered the thermo-mechanical phenomena in the beams made of a functionally graded material. A refined high order theory was used for this purpose, while the in-plane displacement field was depicted by polynomial and exponential expressions. The numerical results so obtained were compared with the solutions of other authors. Sobhy [17] used a new four-variable shear deformation plate theory to depict vibration and buckling of functionally graded sandwich plates supported by elastic foundations. The equations of motion were derived based on Hamilton's principle. The validity of the theory was verified by comparison of the obtained results with the previous ones. Sarangan and Singh [18] developed several new shear deformation theories applicable to analyzing the static, buckling, and free vibration behaviour of laminated composite and sandwich plates. The theories ensure zeroing of the transverse shear stresses at the plate's outer surfaces. The accuracy of the models was positively verified by comparison with the results of 3D elasticity solutions and existing theories. Chen et al. [19] investigated the free and forced vibration of functionally graded porous beams. The Timoshenko beam theory with consideration of the effect of transverse shear strain allowed to derive the equation of motion. The approach enabled the effective calculation of natural frequencies and transient dynamic deflections for the porous beams subjected to various loading conditions. Singh and Singh [20] dealt with laminated and three dimensional braided composite plates. The authors developed two new shear deformation theories for this purpose. The governing differential equations were formulated based on the virtual work principle. The results obtained with the use of the finite element method confirmed good effectiveness of both proposed theories. Shi et al. [21] formulated a new shear deformation theory applicable for free vibration and buckling analysis of laminated composite plates. The theory ensures disappearance of the shear stresses at the surfaces of the plates. Moreover, the shear correction factors are not required. The solutions available in the literature confirmed high accuracy and efficiency of the new method. Thai et al. [22] presented a simple beam theory used for the analysis of static bending and free vibration of isotropic nanobeams. The governing equation was derived based on the equilibrium equations of elasticity theory. Analytical solutions were obtained for nonlocal beams, imposing various types of boundary conditions. Verification has shown good accuracy and effectiveness of the theory. Pei et al. [23] worked out a modified higher-order theory of functionally graded beams using the principle of virtual work. The theory draws a distinction between the centroid and the neutral point of the cross section. In addition, the relation with the traditional higher-order theory is explained, which simplifies a comparative study on various higher-order beam theories. Kumar et al. [24] analyzed functionally graded material plates using two own new higher order transverse shear deformation theories. The energy principle was used to derive the governing differential equation of the plate. The obtained results of deflection and stresses were compared with other published data. The effects of various load types, span to thickness ratio, and grading index were investigated. Magnucki and Lewiński [25] considered simply supported beams with symmetrically varying mechanical properties in the depth direction, subjected to various load types-from uniformly distributed to concentrated one. The deformation of a beam planar cross section after bending was determined based on an own nonlinear "polynomial" hypothesis. The differential equation of equilibrium was formulated based on the definitions of bending moment and shear transverse force and then solved for several beam examples. Magnucki et al. [26] proposed a new formulation of the functions determining the variation of mechanical properties of a beam in the depth direction. The approach consists in a generalization enabling to describe homogeneous, nonlinearly variable and sandwich structures with the use of an universal analytical model. The equations of motion were derived based on Hamilton's principle and analytically solved. The results were verified by FEM computation. Katili et al. [27] proposed a higher-order two-node beam element developed to solve the static and free vibration problems. The Timoshenko beam theory was modified with a view to proper consideration of the transverse shear effect. Effectiveness of the approach was verified by comparison with other data published in the literature. Lezgy-Nazargah [28] developed a global-local shear deformation theory accurately predicting the static and dynamic behaviour of thin and thick layered curved beams. Variation of the shear stress in the beam thickness direction is approximated by a parabolic function. Zeroing of the shear stress on the beam boundary surfaces is ensured without the need for a shear correction coefficient. The results obtained from static and free vibration computations are positively validated by the ones calculated with FEM.
The main goal of the present paper consists in improving the shear deformation theory of bending in case of symmetrically varying mechanical properties of the material in the depth direction of the cross section. The individual nonlinear function of deformation of the planar cross section is proposed. The improved shear deformation theory is applied to the exemplary beams, the analytical model of which is developed. The analytical model of these beams is developed. The analytical results are compared with those obtained by a FEM numerical approach. The presented problem of bending beams with consideration of the shear effect is a continuation of the research submitted by Magnucki and Lewinski [25] and Magnucki et al. [26].

Analytical modelling of the beam bending with consideration of the shear effect
The subject of the study is a simply supported beam of symmetrically varying mechanical properties in the depth direction. The beam of length L and rectangular cross section of depth h and width b is subjected to three-point bending and situated in the Cartesian coordinate system xyz (Fig. 1).
Taking into account the papers by Magnucki and Lewinski [25] and Magnucki et al. [26], the symmetrically varying Young's modulus in the depth direction of the beam is assumed in the following form: where the dimensionless function of the symmetrically varying Young's modulus is An example graph of the dimensionless function of the symmetrically varying Young's modulus in the depth direction of the beam is shown in Fig. 2.
The deformation of a planar cross section after bending of the beam is shown in Fig. 3. The upper and lower surfaces of the beam are free from shear stresses, therefore, the line depicting the shape of this deformation is perpendicular to these surfaces.
The longitudinal displacement according to Fig. 3 takes the following form: where  The conditions for the DFD f d (η) according to Fig. 3 are as follows: Therefore, the strains are and consequently the stresses in accordance with Hooke's law are Taking into account the paper [18] the Poisson's ratio value ν of the beam is assumed to be constant. The elastic strain energy of the beam with consideration of the expressions (5.1), (5.2), (6.1), (6.2) and after integration in the depth direction is in the following form: where Therefore, the first variation of the elastic strain energy is The shear effect arising in the beam, especially in case of three-point bending, is significant (Fig. 4.) The work of the load according to Fig. 4 and its first variation are as follows: where T (x)-shear force-transverse force. Therefore, based on the principle of stationary total potential energy δ(U ε -W ) =0, the system of two differential equations of equilibrium of the considered beam is obtained in the following form: Taking into account the bending moment M b (x) = A yσ x (x, y)dA and the expression (6.1), after simple transformation one obtains the equation The differential equations of fourth order (10.1) and second order (11) are equivalent. Therefore, Eqs. (11) and (10.2) are governing equilibrium equations of the bending beams with consideration of the shear effect. This system, after simple transformation, is reduced to one differential equation in the following form: where: ξ = x/L-dimensionless coordinate, λ = L/ h-relative length of the beam, α = 1 2(1+ν) The shear-transverse force of the three-point bending is as follows: The solution of the Eq. (12) for the first interval of the shear-transverse force (0 ≤ ξ < 1/2) is in the following form: where the integration constants: C 1 = 0 result from the condition dψ/dξ | 0 = 0, and Thus, the dimensionless longitudinal displacement of the upper and lower surfaces of the beam for the first interval (0 ≤ ξ < 1/2) is as follows: and for the second interval (1/2 < ξ ≤ 1) The DFD of the planar cross section of the beams or DFD of the straight line normal to the neutral surface of the plates and shells are a basis for consideration of the shear effect. The DFD for bending beams with symmetrically varying mechanical properties in the depth direction of the cross section is assumed in the following form: where C 0 = 1/2 0 1−4η 2 ks f e (η) dη-coefficient, k s -unknown exponent (positive real number). This function satisfies the conditions (4). The analytical model of the beam bending is formulated with consideration of the shear effect pertains to three-point bending.
The results of the calculations, i.e. the values of the exponent k s of the deformation function (17), the shear coefficient C s (22), and the dimensionless maximum deflectionv max (21), for nine selected beams of relative length λ = 10 and Poisson's ratio ν = 0.3, are specified in Tables 1, 2, and 3.
The shapes of the deformations of planar cross sections-the function (17) of the selected beams are shown in Fig. 7.
The graph of the shear-transverse force (25) for the example beam B-2-2 is shown in Fig. 8. The shear force diagrams for the other eight beams are similar to the above diagram (Fig. 8).
It may be noticed that in the particular case of a homogeneous beam e 0 =1, and in the result the exponent k s =1, shear coefficient C s =3.085, dimensionless maximum deflectionv max = 25.771, and maximum dimensionless shear stress (24)τ max = 0.75.

Comparative analysis-analytical approach
Taking into account the paper by Magnucki and Lewinski [25] similar analytical studies are carried out for comparative purposes for example beams, with consideration of the deformation function developed in this paper in the form where: β =1/(1+k sc )-parameter, k sc -even exponent (2 ≤ k sc )-natural number. Consequently, the derivative of the function is The maximum relative deflection of the beam with consideration of the deformation function (26) is consistent with the expression (20), thereforeṽ    where the dimensionless maximum deflection is and the shear coefficient is The shear stress (6.2) with consideration of the functions (15) and (27) for the first interval (0 ≤ ξ < 1/2) is as follows: Therefore, the dimensionless shear stress for ξ = 0 reads The results of the calculations, i.e. the values of the exponent k sc of the deformation function (26), the shear coefficient C sc (30), and the dimensionless maximum deflectionv max (29), for nine selected beams of relative length λ = 10 and Poisson's ratio ν = 0.3, are specified in Tables 4, 5, and 6. The graphs of the dimensionless shear stresses (32) in the depth direction of the selected beams (B-1-1, B-2-2, B-3-3) are shown in Fig. 9.
Comparison of the dimensionless values of the maximum deflection of both series of the results (Tables 1, 2 Young's modulus of the middle part of B-3-1 is very small, and, in consequence, the deformation function (26) does not represent the beam behaviour correctly. This is the case since the value k sc must not be less than 2.
What concerns the shear stresses, their plots (Figs. 6 vs 9) differ significantly. This allows to conclude that the function (26) does not depict the shear stresses properly.

Numerical FEM study of bending of the beams with symmetrically varying mechanical properties
The beams are modelled using the SolidWorks software. The model of the B-2-2 beam is shown in Fig. 10. Symmetry of the beams and their loads allow to take into account a quarter of the whole structure. The model is divided into the layers distinguished by different Young's modulus values. In case of the B-2-2 beam its mechanical properties in the relatively large middle part remain unchanged. Young's modulus significantly varies only in the parts located in the vicinity of the beam's outer surfaces. Therefore, ten relatively thin layers are just there located, each of them having different Young's moduli.
The B-2-2 model is composed of over 5 million 3D tetrahedral finite elements with 4 Jacobian points. The number of the nodes reaches nearly 7 million. Additionally, a finer mesh was used with a view to find whether the computation accuracy is sufficient. The results so obtained were equal to the primary ones with the accuracy of five significant figures. These numbers of the finite elements and nodes are smaller for the other eight beams. A part of the B-2-2 meshes is shown in Fig. 10. The mesh being coarse in the middle part becomes much finer as it approaches the surfaces.
The beam is located in the Cartesian coordinate system, the origin of which is situated at the beginning of the beam neutral axis. The x-axis coincides with the neutral axis, the y-axis points down, the z-axis is normal to the longitudinal middle plane of the beam.
The following boundary conditions are imposed, with a view to ensure proper behaviour of one fourth of the beam: (i) The beam model is simply supported at its edge (i.e. for x =0), hence, the y displacements of the wall coinciding with the yz -plane are zero. The model is loaded with the force 1/4F downward directed and applied to the surface mentioned in the point 2 above (Fig. 11).
The FEM calculation has been carried out with the following data: E 1 = 2 · 10 5 MPa, ν = 0.3, L = 1000 mm, b = 60 mm, h = 100 mm.  The dimensionless shear stresses shown in Fig. 12 for three variants of the beam are plotted along the AB line (Fig. 10).
Maximum values of the shear stress (for η = 0) are very close to those of Fig. 6. The difference does not exceed 2.5%. Moreover, the respective graphs are very similar to each other. It should be noticed that in the FEM approach Young's modulus takes discrete values, different in each of the layers, instead of the continuous pattern depicted in Fig. 5a-c. This is the reason for these differences.
The values of the dimensionless maximum deflection determined FEM numerically for nine selected beams of relative length λ = 10 and Poisson's ratio ν = 0.3 are specified in Tables 7, 8, and 9.
The deflection results obtained analytically with consideration of the deformation function (17) (Tables 1, 2, 3) and numerically comply very well to each other. The maximum difference between them amounts to 1% in the case of the B-3-3 beam.

Conclusions
The presented analytical studies of bending beams with symmetrically variable mechanical properties in their depth direction allow to formulate the following conclusions: • The assumed function of symmetrically varying Young's modulus (1) describes a family of beams with the structures from homogeneous to approaching the three-layer. • The assumed function of deformation of the planar cross section (DFD) of the beams in the form (17) is original and generalizes the theory of beams. This function meets the condition of zeroing shear stresses on the upper and lower surfaces of beams. • This function may be used in analytical modeling of the plates and shells to describe the deformation of a straight line normal to the neutral surface.
The proposed individual shear deformation beam theory formulated based on the function of deformation of a planar cross section (17) accurately describes the distribution of shear stress in the depth direction of the beam.
Comparison with the FEM numerical study gives evidence of good compliance of the results.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.