# Accurate Solution for Orthotropic Beams with Continuously Variable Cross Sections

Authors: Weihua Fang; Zhi Chen; Yepeng Xu
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

Orthotropic materials and components are widely used in various engineering fields, such as aerospace, energy, infrastructure, and water conservancy projects. The orthotropic components such as slabs and beams cannot be investigated accurately by using traditional beam theory because it is unable to take more elastic constants into account. The precise stress distribution of orthotropic simply supported beams of any thickness is determined in this article. Meanwhile, all undetermined coefficients of the stress function are derived by performing a Fourier sine series expansion upon the bottom and top surfaces of the beam. It shows that the resolution corresponds well with the numerical analysis findings from the finite element method. It indicates rapid convergence and can obtain high-precision and high stress analysis results, which has a good application prospect in engineering.


---

## Body

## 1. Introduction

Civil and water conservancy projects such as dykes include a large number of slabs, beams, and column structures. With the development of science and technology, especially the continuous deepening of intelligent diagnosis and digital twin research on structural safety, it is of great significance to integrate structural analytical solutions to evaluate the safety status of structures and achieve accurate safety assessment of water projects. With the continuous application of anisotropic construction, composite materials, and composite structures, it is also necessary to consider the structural anisotropy and variable cross section characteristics in order to more closely approximate the hydraulic engineering components such as dykes. Due to the needs of mechanical, civil, and various engineering practices, accurate solutions to the plane stress (plane strain) problem have always been an important direction. The accurate solutions of thin isotropic beams and moderately thick beams can be generally determined with negligible error by the Euler–Bernoulli and Timoshenko beam theories, respectively. However, the transverse shear deformation in beams will have greater influence on the solutions when the thickness of a beam increases. In such instances, the conventional beam concept, founded on the Euler–Bernoulli theorem, would provide incorrect or even erroneous findings. Due to the difficulty in determining the shear coefficient, it is difficult to apply the classical Timoshenko beam theory to solve the problems in beams with variable cross sections. The higher-order beam theory has better accuracy in solving the problem of variable thickness beams. It is worth noting that traditional orthotropic beam theories frequently fail to account for all elastic constants in solution. The two-dimensional elasticity theory can provide more precise conclusions in this case.For isotropic beams with constant thickness, Timoshenko and Goodier [1] have systematically investigated. To evaluate the displacements and stresses in anisotropic beams, Hashin [2] employed polynomials in two coordinate variables to generate stress functions. Rao and Rao [3] investigated the post-buckling problem in an S-S beam under small strain conditions, ignoring shear and cross-sectional deformation. For beams fixed at both ends, Ding et al. gave corresponding analytical solutions for different material properties, in which the stress function is given in biharmonic form [4, 5]. Yaghoubi et al. [6, 7] studied the buckling of centrosymmetric anisotropic beams and conducted variational formulations and isogeometric analysis for the dynamics of anisotropic gradient beams. Recently, Niiranen et al. [8] gave an overview on the variational formulations, model comparisons, and numerical methods for Euler–Bernoulli micro and nanobeam models. Relevant research methods have also been applied to functional gradient beams in recent years. For instance, Ding et al. [9] gave the elastic solution of two-dimensional anisotropic functionally gradient beam. Ying et al. [10] achieved the analytical solution for free vibration of a functionally gradient beam resting on Winkeler–Pasternak foundation. Chu et al. [11] gave an elastic solution for two-dimensional elastic beams and functional gradient material beams under stress and bending conditions. Li et al. [12–14] studied the bending, vibration, and large deformation of functionally gradient beams with variable cross-sections using the two-dimensional elasticity theory. Lin and Muliana [15] and Su et al. [16] studied the electromechanical and vibrational properties of gradient beams. In terms of plane stress problems, Benguediab et al. [17] highlighted the analytical solution of cantilever beams of functional gradient materials loaded by uniform loads. Li et al. [18] studied the solution of the nonlinear bending of functional gradient beams based on the kinematic theory of Euler–Bernoulli beams.The theory of Timoshenko beams or the theory of Euler beams was used in most early studies to solve variable thickness beams. For example, Zenkour [19] investigated the elastic properties of orthotropic beams with uniform and variable thicknesses. Nonlinear analysis of any highly deflective composite Timoshenko beam under arbitrary boundary conditions was performed using the analog equation method [20]. Based on the Euler–Bernoulli theory, Grover [21] obtained a closed formula for the lateral vibration of a uniform, isotropic, heat-conducting Kelvin–Voigt type variable thickness viscoelastic thin beam. In Cazzani et al.’s overview paper [22], we can have an understanding into the contributions of Gustav Kirchhoff to the dynamics of tapered beams. Xu and Zhou established an elastic solution of multi-span variable thickness isotropic simple support beams [23]. Xu et al. [24] proposed the two-dimensional elasticity solution for bending of functionally graded beams with variable thickness. For the transient elastic dynamic response of a simple support beam of any thickness under arbitrary transverse load, a semi-analytical method has been proposed by Hasheminejad and Rafsanjani [25]. Boreyri et al. [26] used differential transformation techniques to study the vibration of gradient beams with exponentially increasing thickness on Winkler foundations. An analytical calculation of variable section cantilever beams under elastic constraints was performed by Akbarzade and Farshidianfar [27]. A two-dimensional elastic solution to the vibration problem of a simple support beam of arbitrary continuous thickness is given, and the nonlinear vibrations of an axial functional gradient Euler–Bernoulli beam with uneven cross-sectional shape were studied by Sinir et al. [28]. In recent studies, different numerical methods are employed to compute the static and dynamic behavior of kinds of beams including orthotropic beams and beams with variable thickness. Taken as typical examples, Cazzani et al. [29] conducted isogeometric analyses of Timoshenko beams and strongly curved beams. Balobanov and Niiranen [30] proposed locking-free formulations and isogeometric analysis for Timoshenko beam models. Giorgio [31] developed a three-dimensional nonlinear model to analyze the dynamic problem of a Kirchoff rod. Zhang et al. [32] analyzed the stresses of orthotropic laminated beams subjected to thermomechanical loads, and most recently, Kim [33] proposed asymptotic solutions including boundary layers for orthotropic beams. Furthermore, Li et al. [13] used the separation variable method and the Laplace transformation to analyze the free vibration of functionally gradient beams supported on the Pasternak elastic foundation, and in Li et al.’s work [34], the nonlocal method is employed to conduct thermoelastic analysis of functionally graded beams. However, the study of the elastic analytical solution of the variable thickness orthotropic beam has not yet been reported.The stress function of a continuous variable thickness orthotropic beam that meets the conditions of a continuously variable section that governs differential equations and simple constraints is analyzed. Based on Fourier series expansion along the beam’s bottom and top surfaces, the solutions’ unknown coefficients can be estimated. This approach may be used to solve the stress distribution of any continuous variable thickness orthotropic (including isotropic) beam. These findings can be used as standards or criteria for comparing or verifying the correctness of solutions derived from various approximate analysis and numerical methods related to classical beam theory.


## 2. Governing Equations and Methodology

As shown in Figure1, a simple two-dimensional simply supported (S-S) beam is taken as an example. The top surface is straight and the bottom surface is curved, denoted by fx. The beam’s total length is L, and the thickness is H at the point of x = 0. The top surface is subjected to uniform load qx.Figure 1
Simply supported beam with variable thickness of arbitrary shape.In thex-y plane, the beam’s constitutive equations made of orthotropic material can be given as(1)∂u∂x=s11σx+s12σy,∂v∂y=s12σx+s22σy,∂u∂y+∂v∂x=s66τxy,in which σx, σy and τxy are the normal stresses and shear stress. u and v are the displacement components along x and y directions, respectively. sij (i, j = 1, 2, …, 6) denote the elastic compliance constants.According to the classical elastic theory, the concept of stress functionϕ can be introduced to describe the stress as follows:(2)σx=∂2ϕ∂y2,σy=∂2ϕ∂x2,τxy=−∂2ϕ∂x∂y.Also, the compatibility equation can be accordingly written as(3)s22∂4ϕ∂x4+2s12+s66∂4ϕ∂x2∂y2+s11∂4ϕ∂y4=0.Suppose that the stress function can be written as(4)ϕ=∑m=1∞sinmπxLfmy.It is obvious that the expression of stress function in (4) satisfies the conditions of simply supported beams.Letam=mπ/L, and combine (3) and (4); one can obtain(5)s11fm4y−2s12+s66am2fm2y+am4s22fmy=0,where fmiy denotes the ith order derivative of fmy. The square roots of (5) can be denoted as(6)ε1m2=mπL22s12+s66+2s12+s662−4s11s222s11,ε2m2=mπL22s12+s66−2s12+s662−4s11s222s11.The solution form of (5) depends on the sign of the term 2s12+s662−4s11s22. Thus, three possible solutions are presented for (5).Solution 1.
When2s12+s662−4s11s22>0,(7)fmy=Amshλ1amy+Bmchλ1amy+Cmshλ2amy+Dmchλ2amy,(8)σx=π2L2∑m=1∞m2sinmπxLλ12shλ1amyAm+λ12chλ1amyBm+λ22shλ2amyCm+λ22chλ2amyDm,σy=−π2L2∑m=1∞m2sinmπxLshλ1amyAm+chλ1amyBm+shλ2amyCm+chλ2amyDm,τxy=−π2L2∑m=1∞m2cosmπxLλ1chλ1amyAm+λ1shλ1amyBm+λ2chλ2amyCm+λ2shλ2amyDm,in which(9)λ1=2s12+s66+2s12+s662−4s11s222s11,λ2=2s12+s66−2s12+s662−4s11s222s11.Solution 2.
When2s12+s662−4s11s22=0,(10)fmy=Amshλamy+Bmchλamy+Cmyshλamy+Dmychλamy,(11)σx=π2L2∑m=1∞m2sinmπxLλ2shλamyAm+λ2chλamyBm+λ2yshλamy+2λLmπchλamyCm+λ2ychλamy+2λLmπshλamyDm,σy=−π2L2∑m=1∞m2sinmπxLshλamyAm+chλamyBm+yshλamyCm+ychλamyDm,τxy=−π2L2∑m=1∞m2cosmπxLλchλamyAm+λshλamyBm+λychλamy+1amshλamyCm+λyshλamy+1amshλamy,Dmin which(12)λ=2s12+s662s11.
Ifλ=1, one can deduce that s11=s22=1/E, s66=1/G, and s12=−μ/E. It means that the material of the beam is isotropic.Solution 3.
When2s12+s662−4s11s22<0,(13)fmy=Amshλ1amycosλ2amy+Bmchλ1amycosλ2amy+Cmshλ1amysinλ2amy+Dmchλ1amysinλ2amy,(14)σx=π2L2∑m=1∞m2sinmπxLλ12−λ22shλ1amycosλ2amy−2λ1λ2chλ1amysinλ2amyAm+λ12−λ22chλ1amycosλ2amy−2λ1λ2shλ1amysinλ2amyBm+λ12−λ22shλ1amysinλ2amy+2λ1λ2chλ1amycosλ2amyCm+λ12−λ22chλ1amysinλ2amy+2λ1λ2shλ1amycosλ2amyDm,σy=−π2L2∑m=1∞m2sinmπxLshλ1amycosλ2amyAm+chλ1amycosλ2amyBm+shλ1amysinλ2amyCm+chλ1amysinλ2amyDm,τxy=−π2L2∑m=1∞m2cosmπxLλ1chλ1amycosλ2amy−λ2shλ1amysinλ2amyAm+λ1shλ1amycosλ2amy−λ2chλ1amysinλ2amyBm+λ1chλ1amysinλ2amy+λ2shλ1amycosλ2amyCm+λ1shλ1amysinλ2amy+λ2chλ1amycosλ2amyDm,in which(15)λ1=2s11s22+2s12+s664s11,λ2=2s11s22−2s12+s664s11.
Considering the change of beam thickness and the uniform load on the top surface, the unknown coefficients in (8), (11), and (14) can be solved.
Figure1 shows that the uniform load qx acts on the top surface of the beam, so the boundary conditions can be written as follows:(16)τxy=0,σy=−qx.
The tangent and normal stresses should be zero according to the curved lower surface. Thus, the boundary conditions are(17)lxσx+mxτxy=0,mxσy+lxτxy=0,where(18)lx=cosN,x=−dfx/dx1+dfx/dx2,mx=cosN,y=11+dfx/dx2.
Multiplying (16) and (17) by sinnπx/L and integrating them regarding the range of 0 to L yields(19)∫0LsinnπxLτxydx=0,∫0LsinnπxLσydx=−∫0LsinnπxLqxdx,∫0LsinnπxLlxσx+mxτxydx=0,∫0LsinnπxLmxσy+lxτxydx=0.
Combine (8), (11), and (14) and take the term from 1 to N+1 of the series to obtain the following equation:(20)A111⋯A1N1B111⋯B1N1C111⋯C1N1D111⋯D1N1⋮⋱⋮⋮⋱⋮⋮⋱⋮⋮⋱⋮AN11⋯ANN1BN11⋯BNN1CN11⋯CNN1DN11⋯DNN1A112⋯A1N2B112⋯B1N2C112⋯C1N2D112⋯D1N2⋮⋱⋮⋮⋱⋮⋮⋱⋮⋮⋱⋮AN12⋯ANN2BN12⋯BNN2CN12⋯CNN2DN12⋯DNN2A113⋯A1N3B113⋯B1N3C113⋯C1N3D113⋯D1N3⋮⋱⋮⋮⋱⋮⋮⋱⋮⋮⋱⋮AN13⋯ANN3BN13⋯BNN3CN13⋯CNN3DN13⋯DNN3A114⋯A1N4B114⋯B1N4C114⋯C1N4D114⋯D1N4⋮⋱⋮⋮⋱⋮⋮⋱⋮⋮⋱⋮AN14⋯ANN4BN14⋯BNN4CN14⋯CNN4DN14⋯DNN4⋅A1⋮ANB1⋮BNC1⋮CND1⋮DN=S11⋮SN1S12⋮SN2S13⋮SN3S14⋮SN4,in which Aijk, Bijk, Cijk, Dijk, and Sik can be determined using the integration. Then, A1… AN, B1… BN, C1… CN, D1… DN in (20) can be established. Combining (8), (11), and (14), the beam’s stress distribution can be determined.


## 3. Convergence and Comparison Studies

In the following studies, the material properties (unit: m2.N−1) of orthotropic beams are fixed at s11=8×10−11, s12=4×10−11, s22=8×10−11, and s66=9×10−11, unless otherwise stated. The matrix elements in (20) were numerically computed by the piecewise Gaussian quadrature. To verify the accuracy of the proposed method, the convergence of the present solutions for an S-S wedge-shaped orthotropic beam with linearly varying lower surface as shown in Figure 2 is studied firstly. The minimum value of H/L is 0.05, and a uniformly distributed load q acts on the top surface.Figure 2
Simply supported beam with linear thickness variation.As listed in Table1, solutions to four case terms have been examined. The stress values under different H1/H values are calculated. At x=L/4,y=H/4, the stresses are also given in Table 1. By comparing the results of N = 80 and N = 100, it can be found that the results have stabilized. The stress difference between N = 60 and N = 80 results is not more than 0.5% under certain significant digit accuracy. This validates the rapid convergence of the proposed method. Therefore, N = 80 is used in the following calculations. Special attention is required that the number of truncation terms available in the calculation is not infinite, which will affect the number of significant digits of the computer used. A large number of truncated terms leads to ill-conditioned solutions.Table 1
Stress convergence of simply supported orthotropic beams atx=L/4,y=H/4, and H/L=0.05.
H1/HN = 40N = 60N = 80N = 100σx/q1−112.3−112.4−112.4−112.41.25−105.3−105.4−105.4−105.41.50−98.5−98.6−98.6−98.6σy/q1−0.832−0.848−0.843−0.8431.25−0.822−0.838−0.833−0.8331.50−0.821−0.836−0.832−0.832τxy/q15.625.635.635.631.254.324.314.314.311.503.313.303.303.30It should be noted that the present elasticity solutions are also valid for the beam made of isotropic materials. Figure3 shows a simply supported isotropic beam with symmetrical linear thickness distribution.Figure 3
Simply supported beam with linear variation bottom surface in two directions.The boundary condition shown in Figure3 is that a uniformly distributed load q acts on the top surface, and the two ends are simply supported, and the maximum value of H/L is 0.1. In the finite element analysis, Quad 8-Node-82 entity element is adopted. The surface stresses of different H1/H values are listed in Table 2. By comparing the elastic solution with the ANSYS results, the accuracy of the proposed method is verified.Table 2
The stresses ony=0.2H for the beam with linearly changed thicknesses, in which H/L=0.05.
x/LH1/H=1.5H1/H=1.75H1/H=2.0PresentFEPresentFEPresentFE0.000.01700.01500.0140.1−56.4−56.5−52.9−53.0−49.6−49.70.2−88.7−88.8−78.5−78.6−69.9−70.00.3−103−103−86.9−86.8−74.0−74.00.4−105−105−84.5−84.5−69.4−69.40.5−99.9−100−78.3−78.3−63.0−63.0


## 4. Numerical Examples

The method presented in this paper can be applied to one class of problems. Several examples are given in this section to illustrate the applicability of the method to other problems. Among them, the first is the example of a simply supported beam with a parabola bottom, as shown in Figure4, which is subjected to a uniform distributed load q, where H/L = 0.05.Figure 4
Concave parabolic simply supported beam with variable section.By comparing the normal stressσx on the top surface, it can be found that the distribution trend of normal stress is the same with different H1/H values, but the values are different. As shown in Figure 5, the greater the depth ratio, the smaller the absolute value of stress. In addition, it can be found that the stress value at x = L/2 increases gradually as the depth ratio increases, that is, the change of beam thickness decreases. In the limit case of H1/H = 1, the stress distribution should be the same as that of the constant section.Figure 5
The evenly distributed stress ofσx/q on the upper surface.The second study is a simply supported wedge-shaped orthotropic beam with linear variation bottom surface, as shown in Figure2. The maximum depth of the beam is H = 0.1 L. H2 is the beam’s thickness at x=L/2. As shown in Figure 6, different values of H/H1=1.5,2 correspond to different stress distribution. In Figure 6(a), the stress σx is zero at y/H2=0.5 and attains the maximum absolute value at y/H2=0 or y/H2=1. Figure 6(b) shows that the maximum absolute value of stress σy is located at y/H2=0. Moreover, Figure 6(c) shows that the shear stress τxy is zero at y/H2=0 and attains the maximum absolute value at y/H2=1.Figure 6
Cross-sectional stress distribution at locationx=L/2: (a)σx/q, (b)σy/q, and (c)τxy/q.
(a)(b)(c)Finally, we give an example of the convex shape parabolic variable thickness simply supported beam, as shown in Figure7. The thickness of the left and right ends of the beam is H = 0.05 L. As in the previous examples, the uniform load q is applied to the beam’s upper surface.Figure 7
Convex parabolic simply supported beam with variable section.By comparing the shear stress distribution on the surface whose coordinatey = 0.5H, it can be found that under different H1/H values, the shear stress τxy decreases with an increase of H1/H values, as shown in Figure 8. Moreover, in the limit case of H1/H = 1, the shear stress follows an exact linear distribution, and as expected, the maximum absolute value of the shear stress occurs at the left and right segments of the beam.Figure 8
The shear stress distribution ofτxy/q on y=H/2.


## 5. Conclusions

An analytical expression of the stress function for orthotropic simply supported beams is given and the unknown coefficients are determined by the upper and lower boundary conditions. The two-dimensional analytical solution obtained in this paper can be used for stress analysis of orthotropic beams of arbitrary continuous varying thickness. Of course, isotropy cases can also be analyzed using this method. Verification indicates that the elasticity solution derived in this paper has high accuracy. When compared to the standard beam theory, the two-dimensional elastic theory used in this study can reflect all material properties, ensuring that the findings achieved are more practical. With these features, the proposed method is suitable for various situations requiring high-precision stress analysis, such as a range of high-precision sensor designs [35–38].



---

*Source: 4777608-2022-10-14.xml*