#### 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.

#### 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 Figure 1, 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 . 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 .

In the *x-y* plane, the beam’s constitutive equations made of orthotropic material can be given asin which , and are the normal stresses and shear stress. and are the displacement components along *x* and *y* directions, respectively. (*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:

Also, the compatibility equation can be accordingly written as

Suppose that the stress function can be written as

It is obvious that the expression of stress function in (4) satisfies the conditions of simply supported beams.

Let , and combine (3) and (4); one can obtainwhere denotes the *i*th order derivative of . The square roots of (5) can be denoted as

The solution form of (5) depends on the sign of the term . Thus, three possible solutions are presented for (5).

*Solution 1. *When ,in which

*Solution 2. *When ,in whichIf , one can deduce that , , and . It means that the material of the beam is isotropic.

*Solution 3. *When ,in whichConsidering the change of beam thickness and the uniform load on the top surface, the unknown coefficients in (8), (11), and (14) can be solved.

Figure 1 shows that the uniform load acts on the top surface of the beam, so the boundary conditions can be written as follows:The tangent and normal stresses should be zero according to the curved lower surface. Thus, the boundary conditions arewhereMultiplying (16) and (17) by and integrating them regarding the range of 0 to *L* yieldsCombine (8), (11), and (14) and take the term from 1 to *N*+1 of the series to obtain the following equation:in which , , , , and can be determined using the integration. Then, … , … , … , … 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: m^{2}.N^{−1}) of orthotropic beams are fixed at , , , and , 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.

As listed in Table 1, solutions to four case terms have been examined. The stress values under different *H*_{1}/*H* values are calculated. At , 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.

It should be noted that the present elasticity solutions are also valid for the beam made of isotropic materials. Figure 3 shows a simply supported isotropic beam with symmetrical linear thickness distribution.

The boundary condition shown in Figure 3 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 *H*1/*H* values are listed in Table 2. By comparing the elastic solution with the ANSYS results, the accuracy of the proposed method is verified.

#### 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 Figure 4, which is subjected to a uniform distributed load *q*, where *H/L* = 0.05.

By comparing the normal stress on the top surface, it can be found that the distribution trend of normal stress is the same with different *H*_{1}/*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 *H*_{1}/*H* = 1, the stress distribution should be the same as that of the constant section.

The second study is a simply supported wedge-shaped orthotropic beam with linear variation bottom surface, as shown in Figure 2. The maximum depth of the beam is *H* = 0.1 *L*. is the beam’s thickness at . As shown in Figure 6, different values of correspond to different stress distribution. In Figure 6(a), the stress is zero at and attains the maximum absolute value at or . Figure 6(b) shows that the maximum absolute value of stress is located at . Moreover, Figure 6(c) shows that the shear stress is zero at and attains the maximum absolute value at .

**(a)**

**(b)**

**(c)**

Finally, we give an example of the convex shape parabolic variable thickness simply supported beam, as shown in Figure 7. 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.

By comparing the shear stress distribution on the surface whose coordinate *y* = 0.5*H*, it can be found that under different *H*1/*H* values, the shear stress decreases with an increase of *H*1/*H* values, as shown in Figure 8. Moreover, in the limit case of *H*1/*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.

#### 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].

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

This study was supported by the project of Water Conservancy Science and Technology Foundation of Jiangsu Province (no. 2021073).