Transverse free vibration of Euler-Bernoulli beam with pre-axial pressure resting on a variable Pasternak elastic foundation under arbitrary boundary conditions

With consideration of pre-axial pressure and two-parameter elastic foundation (Pasternak), a new method was put forward for analysis of transverse free vibration of a finite-length Euler-Bernoulli beam resting on a variable Pasternak elastic foundation. Matrices and determinants corresponding to arbitrary boundary conditions were provided for engineers and researchers according to their own demand. In derivation process of new proposed method, Schwarz distribution was adopted for simplifying the partial derivative of Dirac function and compound trapezoidal integral formula was adopted for discretizing the shape function of beam. For the symmetrical boundary conditions, it was concluded that natural frequency of transverse free vibration obtained by FEM highly agreed with the new proposed method. In contrary, for the asymmetrical cases, the calculation results were different from each other. For solving the ordinary differential equation with nonlinear partial derivative terms of shape function, the key point of new proposed method was to establish stiffness equation set composed of obtained matrices, rather than a single equation on the basis of classical theory. This point should be treated as a great advantage. New proposed method can be generalized to solve more complicated problems, which were illustrated in conclusion and prospect.


Introduction
The model of beams resting on elastic foundation have a wide range of applications in civil engineering, such as underground pipelines, rails, strip foundations, pile foundations, etc. [1][2][3][4][5] Winkler elastic foundation is treated as a single-parameter model and commonly used by engineers and researchers, but this model ignores the shearing interaction between soil elements. Previous studies have shown that [6,7], the assumption always results in a huge deviation. Therefore two-parameter model is introduced to overcome the shortcoming of Winkler foundation model, defined as Pasternak elastic foundation. In addition, structure components in engineering, such as railway or steel reinforced concrete, is often in the pre-axial load state. For example, application of pre-stress can improve the compressive strength of concrete. Axial force has a great influence on dynamic response of Euler-Bernoulli beam, and it should be considered in structural dynamic analysis.
Theoretical analysis method was the most commonly used method, with the development of computer technology in recent years, finite element method and other numerical methods have become more and more popular, and were applied for solving problems, such as uniform beams resting on a nonlinear elastic foundation. But most of these results were obtained on the assumption of space-invariable elastic foundation [8][9][10]. As we know, elastic coefficient of foundation always varied with the longitudinal axis of beam, namely variable elastic foundation. Variable elastic foundation was firstly proposed and solved by the use of Taylor series on the assumption of linearly variable elastic coefficient of foundation along longitudinal axis of beam. The same process was confirmed by Hetényi [11,12]. On the basis of compound trapezoidal integral formula, natural frequency of a finite-length Euler-Bernoulli beam resting on variable elastic Winkler foundation were derived under six boundary conditions [13], but the influence of axial force and shearing interaction between soil elements was neglected. Assuming that elastic coefficient of Winkler foundation varied with the negative fourth power of a polynomial function, analytical solution of a finite-length beam under the action of concentrated and uniform load was derived [14], but the whole process was a static analysis, and the assumption of special space-variable elastic coefficient was aimed to set up Euler equation. In case while elastic coefficient varied arbitrarily, the method failed. In the hypothesis that elastic coefficient of Winkler foundation varied in Taylor series along the longitudinal axis of beam, power series method was applied for structural dynamic analysis of a finite-length Euler-Bernoulli beam [15], but the application of power series method to solve differential equations always brought a large amount of calculation because the analytical solution was generally composed of multiple higher power functions.
In summary, none of the above methods could reflect the influence of variable second-parameter of Pasternak elastic foundation, either directly ignored the influence of second-parameter, or simply assumed that second-parameter was constant. It was well known that the second-parameter must also be variable if foundation stiffness varied along the longitudinal axis of beam. In the present paper this problem would be discussed.
With consideration of pre-axial pressure, transverse free vibration response of a finite-length Euler-Bernoulli beam resting on variable Pasternak elastic foundation was solved by a new proposed method in the present paper, which can't be achieved by previous methods. Separation variable method was used to separate frequency and space domain, and variation of parameters was applied to solve fourth-order ordinary differential equation to obtain analytical expression of structural shape function, and then compound trapezoidal integral formula was introduced to discretize structural shape function and establish the stiffness equations set composed of corresponding matrices and determinants to obtain structural natural frequency under arbitrary boundary conditions. In the derivation process, Schwarz distribution theory was used to simplify the partial derivative of Dirac function, and Schwarz distribution theory was illustrated in Appendix B. The new proposed method can be used to carry out research on transverse free vibration response of a finite-length Euler-Bernoulli beam resting on arbitrarily variable Pasternak elastic foundations; By comparing with finite element method(FEM), it was verified that the new proposed method was extended with good applicability for continuous variable foundation and piece-wise homogeneous foundation. In numerical calculation, the whole process of new proposed method was demonstrated, and the influence of continuous variable coefficients of Pasternak elastic foundation on mode of transverse free vibration under different boundary conditions was discussed. New proposed method was rigorously theoretical and suitable for further popularization.
system, x was the longitudinal axis of beam, , y z was coordinates of beam cross-section, the downward direction of y was positive.

The solution of structural motion equation
Transverse free vibration equation of a finite-length Euler-Bernoulli beam with pre-axial pressure resting on variable Pasternak elastic foundation, represented respectively fourth and second-order partial derivatives of ( , was second-order partial derivatives of ( , ) v x t to t , z EI represented flexural rigidity and A area of structural cross-section, According to principle of mass distribution m A ρ = ; ω represented structural natural frequency.
Transverse free vibration equation in the non-dimensional form was simplified and expressed as, According to the method of separation of variables ( , ) ( If the coefficients of elastic foundation were invariable, analytical solution of Eq(3) can be solved by Fourier integral and Laplace transformation under specified boundary conditions, but in present paper it aimed to solve motion equation of Euler-Bernoulli beam resting on variable Pasternak elastic foundation, so theoretical analysis method failed.
It was assumed that the solution of Eq(3) was composed of ( ) . The left side of Eq(3) was homogeneous differential equations, and general solution ( ) g V x was simplified to obtain analytical expression. In contrast, the right side of Eq(3) was much more complicated, variation of parameters was used to obtain ( ) p V x , theoretical basis of method and whole process of solution was illustrated in Appendix A.

Undetermined coefficients under specified boundary condition
Undetermined coefficients in Eq(4) were assumed as unknown real numbers,and can be determined under specified boundary condition. Herein Pinned-Pinned beam (P-P) was chosen as an example and its corresponding boundary condition was expressed as, Boundary condition of P-P beam was substituted into Eq(4), then undetermined coefficients m B were obtained and expressed as,

Compound trapezoidal integral formula
Due to the complexity of the integrand, well-known Newton-Leibniz integration formula was too difficult for solving the integration in many practical problems. Therefore, many numerical integral formulas were put forward, such as Gaussian integral formula, compound trapezoidal integral formula, and Simpson integral formula. Herein compound trapezoidal integral method was used and integral interval [0,1] was discretized into n equal parts, coordinates of each equal parts was with its corresponding length 1/ n . The approximation of Eq(4) and (6) was expressed as, In order to eliminate inconvenience of second-order partial derivative term in Eq (8) to the subsequent calculation, second-order partial derivative of Eq(8) to x was performed,

Structural natural frequency and mode of transverse free vibration
Eq (8) and (9) were respectively multiplied by Schwarz distribution theory was used for simplifying the partial derivative of Dirac function ( ) (11), and corresponding theory was illustrated in Appendix B. Eq (10) and (11) can be denoted in matrix form. 11 12   Through introducing identity matrix [I],Eq (12) was simplified as, Eq.(13) was defined as stiffness equation. According to classical theory of structural dynamics [16], stiffness equation of transverse free vibration was a single equation. However, Eq(13) was an equation set because second-order partial derivative term of was ingeniously handled in new proposed method and it could be extended to solve ordinary differential equations with first or third-order partial derivative terms. Therefore the method could be generalized to solve more complicated problems. This point should be treated as a great advantage.
In Eq (13) , and its determinant was recognized not equal to 0. The above determinant H was set as 0 for obtaining natural frequency of a finite-length Euler-Bernoulli beam with pre-axial pressure resting on variable Pasternak elastic foundation with consideration of axial force and bending. The key of new proposed method was to find structural natural frequency. In Chapter 3, an example would be chosen to illustrate the whole process of new proposed method. Then each-order natural frequency was respectively substituted into Eq (16) for corresponding mode of transverse free vibration { } d .

Undetermined coefficients and Matrices corresponded to solution of structural natural frequency and mode under arbitrary boundary conditions
In many practical problems, there always exists various boundary conditions. Therefore In the section, for the case of arbitrary boundary conditions, undetermined coefficients m B and corresponding matrices ij C , ij D , ij C  , ij D  have been deduced and provided for engineers and researchers. The derivation process was similar as Pinned-Pinned beam and was not described once more. 6 specified boundary conditions were taken into consideration and divided into symmetrical and asymmetrical conditions. Clamped-Clamped(C-C), Pinned-Pinned(P-P), Free-Free(F-F) were recognized as symmetrical, while Clamped-Free(C-F), Pinned-Free(P-F), Pinned-Clamp(P-C) asymmetrical.

Numerical discussion
For continuous variable Winkler elastic foundation, some analytical solutions have been proposed for certain specific situations. However, for most of cases or piece-wise homogeneous elastic foundation, the most commonly used solution was finite element method (FEM). There existed many corresponding studies. The following section was mainly aimed to demonstrate the whole process of new proposed method to solve structural natural frequency and corresponding mode under arbitrary boundary conditions.
A finite-length Euler-Bernoulli beam was taken into consideration, structural material was steel, and the values of parameters were denoted in Table 1, where all symbols were in accordance with previous definition.

The whole process of new proposed method
In finite element method (FEM), structural natural frequency and corresponding mode can finally be determined by solving the eigenvalues and eigenvectors of structural global stiffness matrix and mass matrix. In present paper, new proposed method was different from FEM.
The first step of new proposed method was to determine the range of structural natural frequency corresponded to minimum and maximum coefficients of Pasternak elastic foundation by comparison with analytical solution of natural frequency. The values in the range was defined as possible values. Then these possible values were substituted into Eq.(16) by small-step increment to find the structural natural frequency point corresponded to determinant 0 H = . In Fig.2 the curve was the values change of H in the whole process of seeking from 1 st to 4 th order structural natural frequency under P-F boundary condition. It was shown that the value of H should be increased and decreased repeatedly. When the value of H decreased to 0, the corresponding abscissa was structural natural frequency. The following notes should be paid attention in calculation process of new proposed method: 1) 0 was a small quantity that may be positive or negative, and it was generally taken as the absolute value. 2) For various-order natural frequencies, the magnitude of a small quantity was different (4 marked points as shown in Fig.2), therefore a constant small quantity cannot be set.
3) The value of determinant H may be very large (between marked points 3 and 4 in Fig.2) or very small (between marked points 2 and 3) in the calculation process, but it did not affect the use of new proposed method, there Latin American Journal of Solids and Structures, 2020, 17(7), e305 11/17 definitely exist structural natural frequency points corresponded to 0 H = in case that values of determinant H increased and decreased repeatedly. 4) In Fig.2 there existed a breakpoint between the 3 rd and 4 th -order natural frequency, please noted that ordinate of the breakpoint hadn't been drawn because of a large value at the breakpoint compared to other points.

The comparison between FEM and new proposed method
The comparison of 1 st to 6 th structural natural frequency obtained by FEM and new proposed method under specified symmetrical and asymmetrical boundary conditions was shown in Table 2-5. At the end of each table's title, symbols in bracket represented values of various parameters. Such as, from ( p G , k and 50F ) 50F indicated that 50 times axis pressure. In present paper maximum difference ratio κ was introduced to demonstrate difference ratio of structural natural frequency obtained by FEM and new proposed method. For example, for the case of C-C boundary condition in Table 2 By comparing Table 2 and 3, it was indicated that axial pressure could suppress structural free vibration; By comparing Table 2 and 4 it was indicated that the second-parameter of Pasternak elastic foundation could increase structural natural frequency. It was concluded that axial pressure and second-parameter of Pasternak elastic foundation should be considered in structural dynamic analysis. By comparing Table 2 and 5 it was indicated that structural natural frequency became smaller as stiffness of elastic foundation decreases. In Table 2-5 structural natural frequency was respectively calculated under two specified symmetrical (C-C and P-P) and asymmetrical boundary conditions(P-F and C-F). It was clearly demonstrated that maximum difference ratio κ was a small quantity under specified symmetrical boundary conditions, and in this case structural natural frequency obtained by FEM agreed with new proposed method. In contrast, for the case of asymmetrical boundary conditions, maximum difference ratio κ has reached to 12.09% with consideration of smaller axial pressure and larger foundation stiffness. It should be paid more attention that for the case of asymmetrical boundary conditions, there existed a certain difference between two methods. The comment can't be made on which method was more accurate. It was mentioned that high-order structural natural frequency of slender beam obtained by FEM had a large error with theoretical solution [17]. In present article another solution was put forward.
It was shown in Fig.3-5 that influence of variable foundation coefficient 0 b on 1 th and 2 nd -order modes of transverse free vibration under 3 specified boundary conditions(C-F, C-C, and F-F). For the case of C-F, it was indicated in Fig.3 that as elastic coefficients of Pasternak foundation gradually decreased along the longitudinal axis, right end of Euler-Bernoulli beam seemed to be subjected to a clockwise bending moment, meanwhile amplitude of mode was moving to the right direction because left end of Euler-Bernoulli beam was fixed; For the case of C-C, in Fig.4. it was clearly shown that shape change of mode was small because both ends of beam were fixed, but the amplitude of mode would shift toward the direction of foundation stiffness reduction. At this time, the mode should be treated as superposition of multiple sinusoidal and cosine functions rather than a single sinusoidal function; For the case of F-F, influence of variable foundation coefficient 0 b on modes of F-F was observed to be much greater than that of C-F and Latin American Journal of Solids and Structures, 2020, 17 (7), e305 15/17 C-C. It seemed to be subjected to simultaneous action of bending moments at both ends of beam. In summary variable elastic coefficients of Pasternak foundation was recognized as an indispensable part in the structural dynamic analysis.

Conclusion and prospect
In present paper a new method was put forward for obtaining natural frequency and corresponding mode of a finite-length Euler-Bernoulli beam with pre-axial pressure resting on variable Pasternak elastic foundation. New proposed method was applicable for continuous variable foundation and piece-wise homogeneous foundation. The matrices and determinants corresponding to arbitrary boundary conditions have been deduced and provided. These formulas could provide guidance for engineers and researchers according to their own demand. 6 Boundary conditions was divided into symmetrical and asymmetrical. For the case of symmetrical boundary conditions, structural natural frequency obtained by new proposed method was highly consistent with FEM. In contrast, for the case of asymmetrical boundary conditions, structural natural frequency between two methods was different., and more attentions should be paid on the asymmetrical conditions. Variable elastic coefficients of Pasternak foundation had a great influence on the shape and amplitude of 1st and 2nd-order transverse free vibration mode under specified C-F, C-C, and F-F boundary conditions, especially for the case of F-F condition. The vibration modes of a finite-length Euler-Bernoulli beam resting on variable Pasternak elastic foundation were no longer a single trigonometric function, but superposition of multiple trigonometric and cosine functions.
Compared with FEM, the new proposed method has many advantages. Firstly, FEM was on the basis of corresponding shape function, which may result in a huge deviation for slender beam [17]. In contrast, the shape function of new proposed method was obtained by variation of parameters. Secondly, for those who only concerned the lowerorder structural natural frequency and corresponding vibration mode, or even first-order [10], new proposed method was a good choice. For example, the critical speed of the train traveling on the track was often related to first-order natural frequency, and for the problem of a finite length Euler-Bernoulli beam with open edge cracks [18]. As we know, the cracks always decreased bending stiffness of beam, and were modeled as the sum of several Dirac functions in governing equation of a finite Euler-Bernoulli beam. New proposed method could be to treated as an analytical solution, because the number of cracks determined the size of determinant. For the case of Euler-Bernoulli beam with several open edge cracks, structural natural frequency could be directly solve as an equation of multiple degree with one unknown by MATLAB. Finally, in the new proposed method the treatment of second -order partial derivative term of shape function provided a reference for solving ordinary differential equations with several nonlinear partial derivatives of shape function.