Nonlinear Vibration of Axially Loaded Railway Track Systems Using Analytical Approach

Abstract

In this paper, the nonlinear vibration of railway track systems resting on elastic foundation has been studied. An axially loaded simply supported Euler–Bernoulli beam resting on a flexible foundation has been considered to provide a mathematical representation of the railway track system. Winkler springs have been used to model the elastic foundation. Nonlinear partial differential equation of the system has been presented and solved. A new approximate analytical solution called Improved Amplitude–Frequency Formulation (IAFF) is proposed to obtain nonlinear frequency of the system and an accurate analytical solution for the whole domain. The first iteration of the IAFF leads to a highly accurate solution in comparison with the exact frequency of the problem. The exact frequency of the problem is also presented, and the results of IAFF are compared and verified. Sensitive analysis of the soil stiffness and loading condition is studied for different parameters. A full comparison of the IAFF and exact solution results are illustrated. It has been proved that the IAFF can be potentially extended to highly nonlinear conservative problems.

Introduction

Rail stress condition is a crucial parameter driving rail safety. Effective rail thermal stress management is crucial for the prevention of rail buckles and pull-a-parts. According to the Federal Railroad Administration, there have been 5977 rail accidents in the United States over the last five years. Around 1% (58) of these are attributed to rail buckling. ¹ In railway engineering, one of the common approaches to describe the railway track system under different loading conditions is to represent it as a beam resting on an elastic foundation. Various approaches have been developed to simulate the soil foundation behavior and to consider soil–structure interaction (SSI) such as Pasternak, Winkler, or Vlasov, Flonenko–Borodich foundations. ² The Winkler model is one of the most common approaches to consider the normal displacement of the structure. To model the railway track systems, the rail is modeled as an infinitely long beam, and the soil medium is considered as a series of parallel linear spring elements. ³ Shamalta and Metrikine ⁴ studied the response of the railway track system and modeled the soil foundation by using Winkler springs. Timoshenko beam theory on the Winkler foundation was developed by Steele ⁵ and Chen and Huang^6,7 under different loading conditions. Kenney ⁸ studied the effect of soil–structure interaction of the railway track system. He developed closed-form solutions for dynamic response of an infinite Euler–Bernoulli beam on the Winkler spring foundation. Zhou ⁹ provided a general solution for the frequency analysis of beams on a variable Winkler elastic foundation. Auersch ¹⁰ investigated the dynamic behavior of the soil-structure interaction of beam-type structures. The soil was modeled as soil–finite and infinite, half-space, and Winkler models in his research. Ruge and Birk ¹¹ extended the work on Timoshenko and Euler–Bernoulli beam models on the Winkler foundation. Approximate dynamic analysis was conducted for high-frequency beam vibrations for both models.

Al Rjoub and Hamad ¹² extended the work on the analytical study of the multi-cracked, axially loaded beam with different boundary conditions. They implemented the Transfer Matrix Method to achieve the mode shape of the beam vibration.

Free vibration of a cracked Euler–Bernoulli beam was studied by Attar. ¹³ Attar developed a new approach to solve the inverse problem Euler–Bernoulli beam with cracks.

Boudaa et al. ¹⁴ presented finite element modeling of a beam resting soil layer. The soil was modeled as a linear and homogeneous isotropic continuum. Shear strain of the beam element and soil foundation were considered simultaneously.

Partial differential equations are used for mathematical representation of beams resting on the linear or nonlinear elastic foundation.^15–20 Generally, finding an exact solution for nonlinear vibration equations is very difficult. Javanmard et al. ²¹ implemented Energy Balance Method (EBM) to the Euler–Bernoulli beam resting on an elastic foundation. They modeled the elastic foundation with Winkler springs and verified their results with numerical solutions. Bayat et al. ²² developed the Variational Approach on the nonlinear vibration of an electrostatically actuated microbeam. Wu et al. ²³ studied on the application of the Homotopy Perturbation Method for nonlinear oscillators with coordinate dependent mass. Cha et al. ²⁴ extended the Perturbation Methods to provide first- and second-order solution for symmetric and asymmetric vibratory systems. Many scientists have been working on other asymptotic methods such as: Parameter Expansion Method, ²⁵ Differential Transform Method, ²⁶ Variational Iteration Method, ²⁷ Homotopy Perturbation Method, ²⁸ Max–Min Approach, ²⁹ and other analytical approaches.^30–32

The main purpose of this paper is to implement a new approximate method called Improved Amplitude–Frequency Formulation (IAFF) to consider the nonlinear vibration of axially loaded Euler–Bernoulli beams resting on Winkler foundations. Exact frequencies of the problem are presented, and the results of IAFF are verified. The effects of the important parameters on the nonlinear frequency of the problem are considered. The results have been shown that the IAFF solution is valid for the whole domain.

Description of the problem

Figure 1 represents a rail track resting on an elastic foundation that models with Winkler springs. A straight rail (beam element) with length $L$ , a cross-section $A$ , a mass per unit length $M$ , moment of inertia $I$ , and modulus of elasticity $E$ that subjected to an axial force of magnitude $\overline{P}$ . It is assumed that the section area is uniform, and the material is homogenous. The beam is modeled according to the Euler–Bernoulli beam theory. A series of parallel springs are used to model the Winkler foundation. Linear stress–strain relationship is considered to describe the soil behavior.

Figure 1.

Schematic representation of an axially loaded simply supported Euler–Bernoulli beam resting on Winkler foundation. ²¹

The basic assumptions of the Euler–Bernoulli beam theory are considered, such as ²⁹ :

The beam is isotropic and elastic.

The beam deformation is dominated by bending, and the distribution and rotation are negligible.

The beam is along as slender with a constant section along the axis.

A mathematical formulation of the axially loaded Euler–Bernoulli beam by considering the mid-plane stretching effect is ²⁹

$$EI\frac{{\partial }^4{W}^{\prime }}{\partial {X^{\prime }}^4}+M\frac{{\partial }^2{W}^{\prime }}{\partial {t^{\prime }}^2}+\overline{P}\frac{{\partial }^2{W}^{\prime }}{\partial {X^{\prime }}^2}+K^{\prime }(X)W^{\prime }-\frac{EA}{2L}\frac{{\partial }^2{W}^{\prime }}{\partial {X^{\prime }}^2}{{\displaystyle {\int }_0^L\left(\frac{\partial W^{\prime }}{\partial X^{\prime }}\right)}}^{2}d{X}^{\prime }=U(X^{\prime },t^{\prime })$$

(1)

where $K^{\prime }$ is an elastic foundation modulus and $U$ is a distributed load in the transverse direction.

Considering the non-conservative forces are equal to zero. Therefore, equation (1) can be written as follows

$$EI\frac{{\partial }^4{W}^{\prime }}{\partial {X^{\prime }}^4}+M\frac{{\partial }^2{W}^{\prime }}{\partial {t^{\prime }}^2}+\overline{P}\frac{{\partial }^2{W}^{\prime }}{\partial {X^{\prime }}^2}+K^{\prime }(X)W^{\prime }-\frac{EA}{2L}\frac{{\partial }^2{W}^{\prime }}{\partial {X^{\prime }}^2}{{\displaystyle {\int }_0^L\left(\frac{\partial W^{\prime }}{\partial X^{\prime }}\right)}}^2\text{d}{X}^{\prime }=0$$

(2)

Following non-dimensional variables are introduced

$$X=\frac{X^{\prime }}{L},W=\frac{W^{\prime }}{R},t=t^{\prime }\sqrt{\frac{EI}{M{L}^4}},P=\frac{\overline{P}{L}^2}{EI},K=\frac{K^{\prime }{L}^4}{EI}$$

(3)

where $R=\sqrt{(I/A)}$ is the radius of gyration of the cross-section. It is assumed that the elastic coefficient of the Winkler foundation is constant $K^{\prime }(X)=K_0$ . Then, equation (1) can be written as follows

$$\frac{{\partial }^{4}W}{\partial W^4}+\frac{{\partial }^{2}W}{\partial t^2}+P\frac{{\partial }^{2}W}{\partial W^2}+K_{0}W-\frac{1}{2}\frac{{\partial }^{2}W}{\partial X^2}{{\displaystyle {\int }_0^L\left(\frac{\partial W}{\partial X}\right)}}^{2}dX=0$$

(4)

If assuming $W(X,t)=v(t)\varphi (X)$ in which $\varphi (X)$ is the first Eigenmode of the beam and using the Galerkin method, then the following governing nonlinear vibration equation of motion for an axially loaded Euler–Bernoulli beam can be obtained

$$\frac{d^{2}v(t)}{d{t}^2}+({\epsilon }_1+P{\epsilon }_2+K_0)v(t)+{\epsilon }_3{v}^3(t)=0$$

(5)

The initial conditions for the center of the beam are

$$v(0)=\Delta ,dv(0)/dt=0$$

(6)

The values of ${\epsilon }_1,{\epsilon }_2$ , and ${\epsilon }_3$ can be obtained as follows

$${\epsilon }_1=\left({\displaystyle {\int }_0^1\left(\frac{{\partial }^4\varphi (X)}{\partial X^4}\right)}\varphi (X)\text{d}X\right)/{\displaystyle {\int }_0^1{\varphi }^2(X)\text{d}X}$$

(7a)

$${\epsilon }_2=\left({\displaystyle {\int }_0^1\left(\frac{{\partial }^2\varphi (X)}{\partial X^2}\right)}\varphi (X)\text{d}X\right)/{\displaystyle {\int }_0^1{\varphi }^2(X)\text{d}X}$$

(7b)

$${\epsilon }_3=\left(\left(-\frac{1}{2}\right){\displaystyle {\int }_0^1\left(\frac{{\partial }^2\varphi (X)}{\partial X^2}{\displaystyle {\int }_0^1{\left(\frac{{\partial }^2\varphi (X)}{\partial X^2}\right)}^2}dX\right)}\varphi (X)\text{d}X\right)/{\displaystyle {\int }_0^1{\varphi }^2(X)\text{d}X}$$

(7c)

The basic idea of IAFF

Consider a generalized nonlinear oscillator in the following form ³³

$$\ddot{v}\mathrm{+}F\left(v\right)=0,v(0)=\Delta ,\dot{v}(0)=0$$

(8)

We use the following two trial functions

$$v_1\left(t\right)=\Delta \cos \left({\omega }_{1}t\right)$$

(9)

and

$$v_2\left(t\right)=\Delta \cos \left({\omega }_{2}t\right)$$

(10)

The residuals are

$$R_1\left(\omega t\right)=-a{\omega }_1{\hspace{0.17em}}^2\cos \left({\omega }_{\text{1 }}t\right)+F\left(a\cos \left({\omega }_{1}t\right)\right)$$

(11)

and

$$R_2\left(\omega t\right)=-a{\omega }_2^2\cos \left({\omega }_{2}t\right)+F\left(a\cos \left({\omega }_{2}t\right)\right)$$

(12)

The original frequency–amplitude formulation reads ³⁴

$${\omega }^2=\frac{{\omega }_1^2{R}_2-{\omega }_2^2{R}_1}{R_2-R_1}$$

(13)

He used the following formulation, and Geng and Cai improved the formulation by choosing another location point ³⁴

$${\omega }^2=\frac{{\omega }_1^2{R}_2\left({\omega }_{2}t=0\right)-{\omega }_2{\hspace{0.17em}}^2{R}_1\left({\omega }_{1}t=0\right)}{R_2-R_1}$$

(14)

This is the improved form given by Geng and Cai. ³⁴

$${\omega }^2=\frac{{\omega }_1^2{R}_2\left({\omega }_{2}t=\pi /3\right)-{\omega }_2^2{R}_1\left({\omega }_{1}t=\pi /3\right)}{R_2-R_1}$$

(15)

The point is: $\cos \left({\omega }_{1}t\right)=\cos \left({\omega }_{2}t\right)=k$

Substituting the obtained $\omega$ into $v\left(t\right)=\Delta \cos \left(\omega t\right)$ , we can obtain the constant k in ${\omega }^2$ equation in order to have the frequency without an irrelevant parameter.

Implementation of IAFF

The following trial functions are considered to solve equation (5)

$$v_1\left(t\right)=\Delta \cos t$$

(16)

and

$$v_2\left(t\right)=\Delta \cos \left(2t\right)$$

(17)

Residual equations are as follows

$$R_1\left(t\right)=\Delta \cos \left(t\right)\left(-1+{\epsilon }_1+p{\epsilon }_2+K_0+{\epsilon }_3{\Delta }^{2}co{s}^2\left(t\right)\right)$$

(18)

and

$$R_2\left(t\right)=\Delta \cos \left(2t\right)\left(-4+{\epsilon }_1+p{\epsilon }_2+K_0+{\epsilon }_3{\Delta }^{2}co{s}^2\left(2t\right)\right)$$

(19)

Considering $\cos \left({\omega }_{1}t\right)=\cos \left({\omega }_{2}t\right)=k,$ we have

$${\omega }^2=\frac{{\omega }_1^2{R}_2-{\omega }_2^2{R}_1}{R_2-R_1}={\epsilon }_1+p{\epsilon }_2+K_0+{\epsilon }_3{\Delta }^2{k}^2$$

(20)

We can rewrite $v\left(t\right)=\Delta \cos \left(\omega t\right)$ in the form

$$v\left(t\right)=\Delta \cos (\sqrt{{\epsilon }_1+p{\epsilon }_2+K_0+{\epsilon }_3{\Delta }^2{k}^2}t)$$

(21)

In view of the approximate solution, the main equation can be rewritten in the form

$$\frac{d^{2}v}{d{t}^2}+\left({\epsilon }_1+P{\epsilon }_2+K_0+{\epsilon }_3{\Delta }^2{k}^2\right)v\left(t\right)={\epsilon }_3{\Delta }^2{k}^{2}v\left(t\right)-{\epsilon }_3{v}^3\left(t\right)$$

(22)

If by any chance equation (21) is the exact solution, then the right side of equation (22) vanishes completely. Considering this approach which is just an approximation one, we set

$${\displaystyle {\int }_0^T\left[{\epsilon }_3{\Delta }^2{k}^{2}v\left(t\right)-{\epsilon }_3{v}^3\left(t\right)\right]\cos \left(\omega t\right)\mathrm{d}t=0,T=\frac{2\pi }{\omega }}$$

(23)

Considering the term $v\left(t\right)=\Delta \cos \left(\omega t\right)$ and substituting the term to equation (23) and solving the integral term, we have

$$k^2=3/4$$

(24)

So, substituting equation (24) into equation (20), the frequency is obtained as follows

$${\omega }_{\mathit{Nonlinear}}=\sqrt{{\epsilon }_1+p{\epsilon }_2+K_0+\frac{3}{4}{\epsilon }_3{\Delta }^2}$$

(25)

$${\omega }_{\mathit{Linear}}=\sqrt{{\epsilon }_1+p{\epsilon }_2+K_0}$$

(26)

We can obtain the following approximate solution

$$v\left(t\right)=\Delta \cos (\sqrt{\left({\epsilon }_1+p{\epsilon }_2+K_0\right)+\frac{3}{4}{\epsilon }_3{\Delta }^2}t)$$

(27)

The ratio of nonlinear frequency ( ${\omega }_{NL}$ ) to linear frequency ${\omega }_L$ is the same as the result was obtained in equation (21)

$$\frac{{\omega }_{NL}}{{\omega }_L}=\frac{1}{2}\frac{\sqrt{4\left({\epsilon }_1+p{\epsilon }_2+K_0\right)+3{\epsilon }_3{\Delta }^2}}{\sqrt{{\epsilon }_1+p{\epsilon }_2+K_0}}$$

(28)

Results and discussion

To verify the IAFF results, a comparison that has been made with the exact frequency of the problem is presented. The exact frequency ${\omega }_{\mathit{Exact}}$ for axially loaded simply supported Euler–Bernoulli beam governed by equation (5) can be derived, as shown in equation (29), as follows

$${\omega }_{\mathit{Exact}}=2\pi /{\displaystyle {\int }_0^{\pi /2}\frac{4\sqrt{2}\Delta \mathit{sin}(t)}{\sqrt{{\Delta }^{2}si{n}^2(t)\left({\Delta }^2{\epsilon }_{3}co{s}^2(t)+2{\epsilon }_1++2p{\epsilon }_2+2K_0+{\Delta }^2{\epsilon }_3\right)}}\text{d}t}$$

(29)

Table 1 represents the comparison of IAFF frequencies and exact frequencies for different normalized soil stiffness and axial loads. The maximum error is less than 0.009%. One of the most advantages of the proposed approximate solution over the numerical integration scheme of the exact solution is providing an understanding of the affecting parameters on the response of the problem.

Table 1.

Comparison of nonlinear frequencies of the IAFF solution with the exact solution for various parameters ( ${\epsilon }_1={\epsilon }_2={\epsilon }_3=1$ ).

$\Delta$	p	$K_0$	${\omega }_{\mathit{IAFF}}$	${\omega }_{\mathit{Exact}}$	Error %
0.5	5	10	4.023369	4.023358	0.00028
1	10	20	5.634714	5.634648	0.00116
1.5	20	50	8.525696	8.525601	0.00112
2	30	100	11.57584	11.57572	0.00104
2.5	50	200	15.99023	15.99012	0.0007
5	80	500	24.48979	24.4893	0.00203
10	100	1000	34.29286	34.28995	0.00847

To better understand the problem and the accuracy of the proposed approach, the results are compared with the numerical solution obtained by Azrar et al. ¹⁹ An intermediate parameter $\beta$ is introduced. This parameter depends on the values of ${\epsilon }_1$ , ${\epsilon }_2$ , ${\epsilon }_3$ , and $p$

$$\beta =\frac{{\epsilon }_3}{\left({\epsilon }_1+p{\epsilon }_2+K_0\right)}$$

(30)

So equation (28) becomes

$$\frac{{\omega }_{NL}}{{\omega }_L}=\sqrt{1+\frac{3}{4}\beta {\Delta }^2}$$

(31)

In dynamical analyses, the first mode of the vibration gives the most behavior of the system. The first mode of the vibration for the simply supported beam is assumed as $\varphi \left(X\right)=\sin (\pi X)$ . By consideration of the trail function, the results of the IAFF are compared with the results obtained by Azrar et al. ¹⁹ for different values of the amplitude and $\beta$ in Table 2. As shown in Table 2, the IAFF results have an agreement with the results of Azrar et al., ¹⁹ and the maximum error is less than 2%. The difference between the IAFF results and Azrar et al. ¹⁹ for the large amplitude of the vibration is related to mid-plane stretching that was ignored by Azrar et al. ¹⁹

Table 2.

Comparison of nonlinear to linear frequency ratio ( ${\omega }_{NL}/{\omega }_L$ ) for simply supported beams.

$\Delta$	β	Present study (IAFF)	Pade approximate P{4,2}19	Pade approximate P{6,4}19	Error (%) IAFF and Pade approximate P{6,4}
1	3	1.80277	1.78468	1.78442	1.028345
1.5	3	2.46221	2.42618	2.42541	1.517269
2	3	3.16227	3.10845	3.10712	1.774956
2.5	3	3.88104	3.80991	3.80802	1.917532

IAFF: Improved Amplitude–Frequency Formulation.

Figure 2 shows the comparison of the displacement time history for two different cases. The motion of the problem is periodical and is a function of the amplitude and soil stiffness. The effects of the soil stiffness from the soft to the stiff conditions on the nonlinear to the frequency ratio of the beam are presented in Figure 3. It is obvious from Figure 3 that the soil condition plays an important role in the nonlinear frequency of the system. This ratio is increased from 1 for very stiff soils to 1.5 for soft soil conditions under the same loading condition. The same trend is indicated for different loading conditions in Figure 4. The nonlinear to linear frequency ratio is decreased for large axial loading values. Figure 5 illustrates the sensitivity of frequency with respect to the axial load and soil stiffness, simultaneously for the normalized values. The frequency of the beam decreases by increasing the axial load and soil stiffness. It has been demonstrated that the IAFF can lead to a highly accurate analytical solution for the whole domain with only one iteration.

Figure 2.

Comparison of IAFF solution of $w(t)$ based on time with the exact solution for simply supported beam (a) $\Delta =0.5,p=20,K_0=200$ and (b) $\Delta =1,p=10,K_0=100$ .

Figure 3.

Influence of $K_0$ on the nonlinear to linear frequency ratio based on $\Delta$ for $p=5$ .

Figure 4.

Influence of axial load on the nonlinear to linear frequency ratio based on $\Delta$ for $K_0=300$ .

Figure 5.

Sensitivity analysis of frequency IAFF solution for various parameters of $\omega ,P\mathrm{and}{K}_0$ .

The effect of $K_0$ on the phase-plan diagram for the case $\Delta =0.5,{\epsilon }_1=1,{\epsilon }_2=1,{\epsilon }_3=1,p=5$ is shown in Figure 6. A sensitive analysis is done to show the effects of ${\epsilon }_1$ on the ratio of the $\frac{{\omega }_{NL}}{{\omega }_l}$ for $\Delta =1,{\epsilon }_2=1,{\epsilon }_3=1,\text{}{K}_0=50$ in Figure 7.

Figure 6.

Effect of $K_0$ on the phase-plan diagram for the case $\Delta =0.5,{\epsilon }_1=1,{\epsilon }_2=1,{\epsilon }_3=1,p=5$ .

Figure 7.

Effect of ${\epsilon }_{{\hspace{0.17em}}_1}$ on the ratio of the $\frac{{\omega }_{NL}}{{\omega }_L}$ for $\Delta =1,{\epsilon }_2=1,{\epsilon }_3=1,\text{}{K}_0=50$ .

Conclusions

A mathematical representation of the axially loaded railway track system has been proposed and solved analytically. A new approximate analytical solution called IAFF has been implemented to achieve an accurate analytical solution for the whole domain. The Winkler spring approach has been used to represent the soil condition to extract the nonlinear and linear frequency of the system. The results of IAFF compared with the exact solution showed the accuracy of the proposed approach. The sensitive analysis has been performed to consider the effects of the axially loading condition and soil stiffness simultaneously. The high accuracy of the first iteration of the proposed approach is the most significant advantage of this method, which is valid for the whole domain.