Analysis of flexural behavior of a simple supported thin beam under concentrated multiloads using trigonometric series method

For accurately predicting the local stresses around the loading points, a stress function in the form of trigonometric series (TS) was applied to obtain the displacements and stresses for simply supported thin beams under concentrated multiloads. The convergences of TS solutions were testified for both displacements and stresses, where the convergence criteria were established and the proper iteration numbers were given. Besides, the accuracy of TS solutions was verified by finite element analysis, where the stress concentration effect is obvious for the shear stresses around the loading sections, not following the parabolic distribution proposed in literatures. Finally, taking thin beams under multi-wheels loads as an example, parameter studies were then performed to examine the effects of wheels' numbers, distances and locations on the flexural response of beams. Numerical results were summarized into a series of curves indicating the distribution of displacements and stresses for various parameters.


Introduction
Simply supported thin beams subjected to concentrated loads are a class of mechanic problems generally encountered in practical engineering, such as the web of I-shaped girders or off-track box girders under wheel loads, as shown in Fig. 1. Unlike the uniform load, the stress concentration around the loading point is generally obvious, which possibly induces the local buckling before overall yielding, so that more attentions need to be paid in design.
Due to the practical importance, the concentrated loading case has been focused by many researchers since the 1903's. The stress of an infinite long beam under two equal and opposite concentrated loads was initially analyzed by Filon using Fourier series [1], and he found that the stresses around the loading point were remarkable and diminished rapidly with the increment of the distance from the loading point. However, the convergence property around the loading point is not clear and the proper iteration numbers are not clearly specified for stresses and displacements. Afterwards, Seewald [2] solved the problems on the infinite long beam loaded by a concentrated force based on the semi-infinite-plane theory, and found out that the vertical stress (normal to the axis) has a good agreement with the exact solution from the Fourier series solution [1], but the warping one (along the neutral axis) at the bottom edge of the beam has a large deviation up to 90.9%. Similarly, Wang [3] solved the warping stress of the deep beam subjected to a concentrated load by using the semi-infinite-plane theory and the superposition principle, however, the error for the warping stress at the top edge reaches to 26.3% in comparison with the finite element analysis (FEA) result. In FE analysis, gradually refined grid will lead to the large increment of local stress at the loading point, which well corresponds to the infinite stress at the loading point obtained from the Fourier series method [1] and the semi-infinite-plane theory [2,3].
According to the simplification of the deformation, many beam theories have been established, including the classical beam theory (CBT), the first-order beam theory (FBT), and the high-order beam theory (HBT). The CBT known as the Euler-Bernoulli beam [4] is only applicable to slender beams. While for thick or deep beams, the CBT underestimates deflection and overestimate natural frequency and buckling load due to the ignorance of the transverse shear deformation effect [5]. The FBT known as the Timoshenko beam theory [6,7] is proposed to overcome the drawback of CBT by accounting for the transverse shear deformation effect for deep beams. In FBT, a shear correction factor (SCF) is needed to compensate the discrepancy between the actual stress state and the assumed constant stress state since the FBT violates the zero shear stress boundary conditions on beam edges [8]. The SCF depends on various parameters, such as geometry configurations, material properties, boundary conditions [9]; and, therefore, it needs a further research [10].
Besides, in order to avoid the use of SCF and get a better estimation of behavior of deep beams, the HBT, represented by various shape functions for the shear stress, has been developed, including the thirdorder theory [11], the trigonometric theory [12], the hyperbolic theory [13], the exponential theory [14,15], the mixed theory [16][17][18][19]. Although these HBTs were initially proposed for plates and shells, application of the shape functions to beams is immediate. Based on the assumption of a high-order variation of axial displacement and a constant transverse displacement, most of HBTs comply with the zero shear stress boundary conditions and produce a non-linear (generally parabolic-shaped) distribution for the transverse shear stress through the beam height. Applying a high-order variation to both axial and transversal displacements [17,18], Carrera [20][21][22][23][24] proposed the Carrera Unified Formulation (CUF). Then, Demasi [25] provides a hierarchical formulation leading to very accurate FE models for beams, plates and shells, in which the stretching effect is automatically taken into account.
A detailed observation on the literature reveals that the transverse shear strain generally varies in the form of parabolic function in most of HBTs. However, this parabolic variation may not be applicable to the shear strain around the loading sections for beams under concentrated loads, due to the stress concentration. Therefore, this paper deals with the flexural behavior of thin beams under concentrated multiloads by using the trigonometric series (TS) method. The emphasis is placed on the convergence property of TS solutions for both stresses and displacements, especially for those around the loading points, where two convergence criteria were established depending on the iteration type, and the proper iteration numbers were given. Also, the accuracy of TS solutions is verified by FEA for both displacements and stresses, especially for the transverse shear stresses around the loading sections, where the stress concentration is obvious, not following the parabolicshaped distribution proposed in literatures. Finally, taking thin beams under multi-wheels loads in parameter study, the effects of wheels' numbers, distances and locations on the displacements and stresses are investigated respectively.

Trigonometric series method
A thin beam with rectangular cross section subjected to a concentrated load is investigated firstly in Fig. 2a. Due to the small variability of stresses within the beam thickness, the original 3D beam model can be simplified into a 2D one with a unit thickness, as shown in Fig. 2, and the load is correspondingly transferred to P/t. For analysis, the coordinate system O-yz is established in Fig. 2b with its original point O set on left top of the beam. The beam is made of a homogeneous, isotropic and linearly elastic material with Young's and shear moduli E and G, respectively. The span is l and the height is h. The concentrated load P/t is applied at z P on the upper edge. The y-/z-axial displacements are v and w, respectively.
Considering the discontinuity for concentrated load, the stress function φ needs to be expanded in the form of trigonometric series (TS), given by where i is the iteration number. α i =iπ/l. A i and B i are the constant coefficients to be determined by boundary conditions. Substituting Eq. (1) into the compatibility equation, a four-order differentiate equation is obtained, where the subscripts (2) and (4) indicate the second-and four-order differentiates over variable y, respectively. The solution then is where c ki (k=1, 2, 3, 4) are constant coefficients. Consequently, the stress function φ is Based on the relations between the stress function φ and the stresses, we have In order to obtain the coefficients c ki , the concentrated load P/t is also expanded in the form of TS. To do this, we assume that the load P/t is acted within a infinitesimal domain [z P -ε/2, z P +ε/2], given by where ε is the size of domain. Substitute the Eqs. (6)-(8) into the stress boundary conditions, we have Furthermore, based on the constitutive equations, (17) the z-/y-axial displacements w and v can be obtained by integrals, given by where functions f(z) and g(y) are to be determined by Eq. (17). Substituting Eqs. (12)- (14) into Eqs. (18)- (19), then to Eq. (17), we then have It's evident that both f(y) and g(z) are linear functions, given by where M, M 1 and M 2 are constants to be determined by boundary conditions. Therefore, the displacements w and v can be obtained by substituting Eqs. (12), (13), (21) and (22) into Eqs. (18) and (19), given by where G=E/2(1+υ), F i =1-2E i /(1+υ). υ is the Poisson's ratio. For the simple supported boundary, with one end of the beam being hinged and the other being allowed to slide freely on a frictionless support, the restraints on displacements are Substitute Eqs. (23) and (24) Therefore, the z/y-axial displacements are It is obvious the y-axial displacements v at both ending sections (z=0 and z=l) keep zero. Therefore, the restraints Eq. (25) on displacement v at points are expanded to the whole ending sections, which corresponds to the restraints in the finite element model in Section 4.
For the multi-loads beam shown in Fig. 3, the loads P j /t are located at z Pj (j=1, 2, …, n). The stresses and displacements can be obtained by superposition.

Convergence test
A simply supported deep beam, subjected to two concentrated loads, is investigated to testify the convergence property of TS solutions. As shown in Fig. 4, the span l is 10 m, the height h is 2.5 m and the thickness t is 2 cm. The Young's modulus E is 210 GPa and Poisson's ratio υ is 0.3. Two loads P 1 and P 2 are 20 kN and 30 kN, acted at 4 m and 6 m, respectively. The tested points include K 1~K3 at mid span and J 1~J3 and L 1~L3 at loading sections, for three heights, respectively at the upper and lower edges and on the neutral axis.
The convergence results from TS method are comprehensively shown in Fig. 5 to Fig. 9 for the z/y-axial displacements w and v, stresses σ z and σ y and the in-plane shear stress τ zy , at the tested points from the upper edge to the lower one, respectively.
Some findings are drawn below (1) All displacements and stresses converge to one certain limit at infinity, except for the y-axial displacements v and stresses σ z and σ y at the upper edge. For the displacement v at the points J 1 and L 1 , as shown in Fig. 6a, the curves increase with a gradual descending slope ratio and cannot converge to a certain limit within 200 iterations. While for the z/y-axial stresses at the point K 1 , as shown in Fig. 7a and Fig. 8a, the curves vary periodically in the form of trigonometric function and have no convergence, and those at the points J 1 and L 1 generally increase linearly to infinity along with the iteration numbers.
(2) The shear stresses at the upper edge, as shown in Fig. 9a, constantly keep zero with the increment of iteration numbers, which obeys the stress restraints in Eq. (10). While for the y-axial stresses and shear stresses at the lower edge, as shown in Fig. 8c and Fig. 9c, they converge to small values, closed to zero, compared with those on the neutral axis.
For practical purpose, convergence criteria should be established to find the proper iteration numbers (PINs) in TS solutions for the displacements and stresses, which largely depends on the type of iteration results.
(1) For the z-axial displacement w, due to the progressive descending increment, it is stipulated that the summation of two adjacent iteration errors (IEs) should be less than the given error bound (EB), where the IE is defined as the relative error between two adjacent iteration results. (2) For the y-axial displacement, z/y-axial stresses and the shear stress, the calculation stability number (CSN), defined as the time for which the summation of two adjacent IEs is successively less than the given EB, is involved besides the former stipulation on IEs. This is due to the appearance of short platform (horizontal line) before the real convergence. For instance, there exists a short platform within 20 iterations for y-axial stresses at the lower edge in Fig. 8c and shear stresses in Fig. 9c, possibly resulting in a fake convergence before the real one coming. Considering the calculation accuracy and efficiency, the CSN is uniformly taken by 10.
Accordingly, the comparisons of PINs from different EBs or CSNs are tabulated in Table 1 for the displacements and stresses at the tested points on the neutral axis and at the lower edge, where the relative errors (REs) between the TS solutions corresponding to two PINs are calculated. The EBs for displacement w are set to 0.01 and 0.001 respectively, and the results with and without considering the CSNs under the same EB (EB=0.001) are compared for others.
Some findings are drawn from Table 1 (1) For the displacement w, the PINs increase largely from 3 to 14 for the point J 3 and from 4 to 24 for others when the EB reduces from 0.01 to 0.001, but all REs are less than 0.5%. This implies the zaxial displacement has a fast convergence speed, so that it can be obtained from the former 3 or 4 items by manual computation. (2) For the displacement v, the PIN is 33 for the point K 3 under EB=0.001 and CSN=10, which means that the summation of two adjacent IEs is successively less than 0.001 for 10 times from the   33th iteration for the displacement v at the point K 3 . Obviously, the REs are smaller than 0.5% at all chosen tested points for the displacement v, therefore, whether or not considering the CSNs, almost makes no difference for the displacement v. This infers that the former 7 items are adequate for the y-axial displacements at the points J 2 and L 2 , even 4 items for those at K 2 and K 3 , which needs no additional stipulation on the calculation stability. (3) For the z/y-axial stresses on the neutral axis, there exists a large deviation between the results with and without considering CSNs at the point K 2 , where RE=150.3% for the stress σ z and RE=44.65% for σ y , both of which are marked by pink highlights. This demonstrates that there occurs a short platform for the fake convergence      before the real one coming for the z/y-axial stresses at point K 2 , as shown in the zoomed pictures in Fig. 7b and Fig. 8b. So the calculation stability should be mandatorily involved when calculating the stresses on the neutral axis, and the IPNs are recommended as 47. While for the y-axial stresses and shear stresses at the lower edge, similar large deviations happen between the results with and without considering CSNs, and the IPNs are at least 26.
Additionally, a scenario is plotted in Fig. 10, where the points J 11 , J 12 and J 13 , at a distance of h/100, h/50 and h/10 respectively, beneath the point J 1 , are selected to analyze the convergence of TS solutions for the z/y-axial stresses around the point J 1 , so are the points K 1 and L 1 .
The convergence results of TS solutions are shown in Figs. 11 and 12 for the y/z-axial stresses around the points J 1 , K 1 and L 1 , respectively.
Some findings are drawn as follows (1) Obviously, the convergence speed for the y-axial stress becomes more and more slower when approaching to the points J 1 and L 1 .
Simultaneously, the y-axial compressive stresses at i=1000 increase significantly from 2.54 MPa to 25.42 MPa for the point J 1 and from 3.81 to 38.12 MPa for L 1 , as shown in Fig. 11a and c. This infers that the y-axial stress tends to get infinity at the loading points, due to the stress concentration effect. (2) For the z-axial compressive stresses in Fig. 12a and c, gradually slowness for convergence speed is obvious near the loading points J 1 and L 1 . The values at i=1000 increase slightly from 3.20 to 3.98 MPa for the point J 1 and from 3.3 Mpa to 4.12 MPa for L 1 , and tend to converge when approaching to the loading points. This means that the z-axial stresses at the loading points can be approximately estimated by the adjacent points. Therefore, the z-axial stresses are about −3.98 MPa for the point J 3 and −4.12 MPa for L 3 . (3) For the stresses at the point K 1 , as shown in Fig. 11b and Fig. 12b, they show a damped trigonometric-function variability, and those at i=1000 tend to converge when getting closed to the loading points. Similarly, the stresses at the point K 1 can be estimated by the  Table 4 Comparisons between TS method and FEA for z-axial displacements at z=l/2 (PIN=24, unit: μm).  Table 5 Comparisons between TS method and FEA for z-axial displacements at z=z P2 (PIN=24, unit: μm).

Accuracy test
Finite element analysis (FEA) is applied to verify the accuracy of TS solutions by software package ANSYS using the element of Solid45. The meshing grid is shown in Fig. 13a, where the meshing size is 0.2 m and a total of 1400 elements are adequate. Fig. 13b to f show the contours of displacements and stresses for the simple supported thin beam under two concentrated loads, respectively, where the measurements and loads are referred to those in the Section 3 and the boundary restraints to Eq. (25).
Comparison results between TS method and FEA are shown in Table 2 for y-axial displacement at the height z=l/2 and in Tables 3-5 for z-axial displacement for three cross sections (z=z P1 , l/2, z P2 ), respectively, where the PIN is referred to Table 1. It is seen that relative errors (RE) between the TS method and FEA are very small. So the TS method has a very high accuracy in calculating the y-axial (vertical) and z-axial (warping) displacements.
For stresses, the comparison results are plotted in Figs. 14-16. In order to guarantee the accuracy of stresses around the loading points, the PIN is uniformly enhanced to 100 for the y-axial stresses and shear ones, not the recommended 27 or 47 iterations in Table 1. It is seen that the TS method has a very high accuracy in calculating axial and shear stresses. Unlike the Euler beam theory or Timoshenko beam theory, the TS method is also capable of estimating the y/z-axial stresses around the loading points, known as the stress concentration effect, as shown in Figs. 14 and 15.
Besides, the shear stress distributions at sections z=z P1 and z=z P2 are approximately parabolic as shown in Fig. 16a, which satisfies the zero shear stress boundary condition on the upper and lower edge of the beam (Eq. (10)). However, as shown in Fig. 13f, the shear stresses around the loading sections may not follow the general parabolic distribution determined by the shape function, proposed by Levinson [11], Touratier [12], Soldatos [13], Karama [14], respectively.
In order to investigate the shear stresses around the loading sections z=z P1 and z=z P2 , the meshing size in FEA model is further diminished to 0.0125 m and the PINs in TS solutions are uniformly enhanced to 200. The comparison results are shown in Fig. 17, where the z c is the distance from the loading points. It is seen that the shear stresses obtained from TS method have a great agreement with the FEA ones, even for those around the loading points J 1 and L 1 . Obviously, the shear stresses increase significantly when gradually approaching to the   loading points J 1 and L 1 , due to the stress concentration, which no longer follows the former parabolic-shaped transverse distribution. So for the concentrated loading case, the TS method is capable of estimating the shear stresses around the loading points, which cannot be done by the existing parabolic-shaped functions.

Parametric study
Based on the TS method, three thin beams subjected to two/four/ eight-wheels loads respectively, as shown in Fig. 18, are applied to investigate the effects of the wheels' numbers, distances and locations on the displacements and stresses. The total force P is 160 kN. The hook is symmetrical with respect to the trolley, resulting in the uniform distribution of wheel loads in each cases. For beams, the span l=30 m, the height h=3 m, the thickness t=0.04 m, the elastic module E=210 GPa and the Poisson's ratio υ=0.3. For the trolley, the distances b=6 m, r=3 m and s=1 m, respectively.

Number of wheels
The influence of wheel numbers on the stresses and displacements are shown in Fig. 19.
Some findings are drawn as follows (1) In Fig. 19a, the z-axial tensile stresses on the lower edge (y=h) almost vary in the form of parabolic function, and get their maximum at mid span for the four-and eight-wheels cases or around the mid span for the two-wheels case. The differences between the three cases are very small for the z-axial stresses, so are the y-axial displacements v on the neutral axis in Fig. 19c and the shear stresses at the cross section z=(l-b)/2 in Fig. 19d. (2) Fig. 19b gives the distribution of y-axial compressive stresses at the height of y=h/4. It is seen that the maximum reduces from 1.5MPa to 0.42 MPa when the wheel number increases from 2 to 8. This infers that more wheels will reduce the maximum of y-axial stresses under the same total force P, but increase the fabrication cost on the trolley. Plus, the reduction on the maximum y-axial stress is approximately proportional to that on one single wheel load.

Distance between wheels
Taking the two-wheels case as an example, relations between the yaxial stresses and the calculation heights are analyzed in Fig. 20 in terms of the distance b between wheels. Some findings are drawn as follows (1) The y-axial stress increases significantly when the calculation height y approaches to the upper edge from h/2 to h/10, due to the stress concentration, which well demonstrates the well-known Saint-Venant principle. (2) For the cases in Fig. 20c and d, the distance between wheels is so close that the two valleys have superposed with each other at mid span, resulting in the local increase of the peak value. For instance, the peak values for y=h/2 increase from 0.61 to 1.14 MPa when the distance b decreases from 5 m to 0.5 m, and those for y=h/4 from 1.6 to 2.6 MPa and those for y=h/10 from 4.17 to 4.48 MPa. (3) In Fig. 20d, the two valleys have totally superposed to just one for the stresses at y=h/4 and y=h/2, but it still has two valleys for those at y=h/10. This implies that the superposition of two valleys is related to not only the distance between two wheels, but also the Fig. 20. Y-axial stresses at different calculation heights for various distance b for the thin beams under two-wheels loads. minimum distance z 0 between zero stresses on both sides of the symmetrical axis, as illustrated in the zoomed picture in Fig. 21a.
The distance z 0 is probably related to both the calculation height and the ratio of height to span of the beam. Fig. 21b gives the relations between the normalized distance z 0 and the calculation height y in terms of the height-to-span ratio. It is seen that both curves get their maximums around the height y=0.6 h, and those with larger ratio h/l are much higher. This implies that the y-axial stresses of deep beams have larger distance z 0 between zero stresses and are more inclined to superpose than those of slender beams for the same calculation height.

Location of wheels
Based on TS method, Fig. 22 gives the displacements and stresses at some specific points varying with the trolley hook's location z t , known as the influence line, for the trolley moving from the left to right for the simply supported thin beams under multi-wheels loads.
Some findings are drawn as follows (1) For the y-axial displacement at y=h/2 in Fig. 22a, it almost makes no difference between influence lines in the two-and four-wheels cases. The maximum is 4.6 mm for z=l/2 at mid span, and 3.25 mm for z=l/4 and 1.41 mm for z=l/10. Similarly, the z-axial stress at y=h in Fig. 22b shows its reduction on the maximum when the calculation section changes from 'z=l/2′ to 'z=l/10'. So attentions should be paid primarily to the mid span in design. (2) The influence lines of shear stresses at central point (h/2, l/2) are plotted in Fig. 22c for three cases. A group of anti-symmetrical distributions, with respect to the central point, are obvious for all three cases. It is seen that the up-and-down fluctuation damps with the increment of wheel numbers. (3) The influence lines of the y-axial stresses at the point (h/10, l/2) are depicted in Fig. 22d for three cases, where all lines are symmetrical about the mid span. For the two-wheels case, the distance between the valley bottom A1 and mid span is 3 m (=b/2), which equals to that between the wheels and the hook, so is the bottom A2. This implies that when the trolley hook is located at the valley bottoms, one of the wheels is just at the mid span, and the y-axial stress at mid span gets the most. Similar situations can be observed for the bottoms B1 and B2 for the four-wheels case and the C1, C2, C3 and C4 for the eight-wheels one. Furthermore, the influence lines for other points at the height y=h/10 can be obtained by figure translation, and similar conclusions can be drawn.

Conclusions
In this paper, the flexural behavior of a simple supported thin beam subjected to concentrated multiloads is investigated using the trigonometric series (TS) method, where the axial displacements, stresses and shear stress were all obtained in the form of TS according to the boundary conditions. The main conclusions are drawn as follows (1) The TS method has a good convergence property for the displacements and stresses for points away from the loading points. However, the convergence speed for the y/z-axial stresses are getting more and more slower when approaching to the loading points, and cannot converge at the loading points due to the stress concentration. (2) The TS method has a very high accuracy in calculating the axial displacements and stresses, even for those around the loading points, which has been well verified by FEA. Besides, the TS method is capable of estimating the shear stresses around the loading sections, where the stress concentration effect is obvious, not following the existing parabolic-shaped distribution in literatures. (3) Taking a multi-wheels loading thin beam as an example, parametric studies show that the increment of wheel numbers will effectively reduce the maximum of y-axial stresses under the same total force. Also, the superposition between y-axial stresses is related to not only the wheels' distance, but also the minimum distance between zero stresses on both sides of the symmetrical axis. Additionally, for the y-axial stress at mid span, the influence lines for all cases will get the most when one of the wheels is just located at mid span.