Hermite Finite Element Method for Vibration Problem of Euler-Bernoulli Beam on Viscoelastic Pasternak Foundation

Pengfei Ji, Zhe Yin^*

School of Mathematics and Statistics, Shandong Normal University, Jinan, China.

DOI: 10.4236/eng.2024.1610025                

Abstract

Viscoelastic foundation plays a very important role in civil engineering. It can effectively disperse the structural load into the foundation soil and avoid the damage caused by the concentrated load. The model of Euler-Bernoulli beam on viscoelastic Pasternak foundation can be used to analyze the deformation and response of buildings under complex geological conditions. In this paper, we use Hermite finite element method to get the numerical approximation scheme for the vibration equation of viscoelastic Pasternak foundation beam. Convergence and error estimation are rigourously established. We prove that the fully discrete scheme has convergence order $O\left({\tau }^2+h^4\right)$ , where $\tau$ is time step size and $h$ is space step size. Finally, we give four numerical examples to verify the validity of theoretical analysis.

Keywords

Viscoelastic Pasternak Foundation, Beam Vibration Equation, Hermite Finite Element Method, Error Estimation, Numerical Simulation

1. Introduction

Beams are one of the most common components in mechanical equipment and construction. With the rapid development of science and technology, a variety of large bridges have come out one after another. Bridges will produce bending deformation and violent vibration under various types of loads, hiding huge safety risks. Therefore, the damping problem must be considered in the design and construction of bridges. The viscoelastic foundation can effectively disperse the structural load to the foundation soil and reduce the deformation, which ensures the stability and safety of the whole structure. Studying the characteristics of viscoelastic foundation beams can adapt to the specific geological conditions and make the foundation design more accurate and reasonable.

The beam vibration model on viscoelastic foundation is widely used in engineering practice. For example, in the field of aerospace and automotive manufacturing, viscoelastic beam structures can both store and dissipate energy, which is of great significance for noise and vibration control. In fact, many common structures in civil engineering, such as long bridges, tall buildings, tunnels and tracks, can be regarded as the viscoelastic foundation beam model. In order to study the static deflection and dynamic response of beams on different viscoelastic foundations, many researchers have used various models to do a lot of research. Senalp et al. [1] investigated the dynamic response of a finite-length Euler-Bernoulli beam on linear and nonlinear viscoelastic foundation which is subjected to the moving concentrated force. Babilio [2] studied the nonlinear dynamics of beams which rest on a linear viscoelastic foundation. Hörmann et al. [3] focused on investigating the initial-boundary value problem for an Euler-Bernoulli beam model characterized by a discontinuous bending stiffness, which is positioned on a viscoelastic foundation. Beskou and Muho [4] investigated the dynamic response of a simply supported elastic beam on viscoelastic Winkler foundation under a point load with variable speed. Elhuni and Basu [5] present a new method for dynamic analysis of Euler-Bernoulli beams resting on multi-layered viscoelastic foundation. Lu et al. [6] obtained the differential equation of the system for lateral vibration of a viscoelastic rotating beam with fractional derivative by using Hamilton principle. Snehasagar et al. [7] studied the effects produced by the viscoelastic modelling of pavement on the dynamics of vehicle-pavement coupled system. Praharaj and Datta [8] studied the transient response of plates on fractionally damped Kelvin-Voigt viscoelastic foundation model.

However, the equation which can describe the beam vibration problem is a fourth-order partial differential equation. It is difficult to obtain the analytical solution of this kind of equation. In view of the important application value of beam vibration model, it also attracts many researchers to study its numerical method. Frýba et al. [9] made stochastic finite-element analysis for the beam on a random foundation with uncertain damping under moving force. Lou [10] proposed several new finite element methods to calculate the section force of Euler-Bernoulli beam on continuous viscoelastic foundation under concentrated moving load. Taeprasartsit [11] developed a finite element model of large amplitude free vibrations of thin functionally graded beams with immovably supported ends. Sánchez et al. [12] demonstrated the numerical simulation of non-symmetric laminated beams using the three-dimensional finite element program. Wang et al. [13] used mixed finite volume element method to solve the time-fractional damping beam vibration problem. Sun et al. [14] used the Hermite finite element method to solve the vibration problem of a beam with time fractional damping term.

Euler-Bernoulli beams are widely used in construction and engineering because of their simple and refined theoretical expressions. In the following, we take a typical Euler-Bernoulli beam as an example to briefly introduce the derivation of the vibration equation of viscoelastic Pasternak foundation beam.

First of all, the axis of the beam without deformation is taken as the x-axis on the premise of ignoring shear deformation and rotation of the section around the neutral axis. Assuming that the beam has a symmetric plane, the direction perpendicular to the x-axis in the symmetric plane is taken as the y-axis. Considering the force on the element with thickness $\text{d}x$ , we can obtain the dynamic equation of the element along the y direction

$\rho A\left(x\right)\text{d}x\frac{{\partial }^{2}u}{\partial t^2}=F_s-\left(F_s+\frac{\partial F_s}{\partial x}\text{d}x\right)+f\left(x,t\right)\text{d}x,$ (1)

where $\rho A\left(x\right)$ is the mass per unit length, $F_s$ is the shear force, and $f\left(x,t\right)$ is the load on the beam. When the influence of shear deformation and section rotation is not considered, the micro-element body meets the condition of moment balance. Taking any point on the section as the center of moment, we get

$\left(M+\frac{\partial M}{\partial x}\text{d}x\right)-M-F_s\text{d}x-f\left(x,t\right)\text{d}x\frac{\text{d}x}{2}=0,$ (2)

where $M$ is the bending moment. Ignoring high order small quantities, we can derive

$F_s=\frac{\partial M}{\partial x}.$ (3)

According to the analysis of material mechanics and the relationship between bending moment and deflection, we have

$M\left(x,t\right)=EI\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2},$ (4)

where $E$ is the Young’s modulus, $I$ is the area moment of the beam with respect to the neutral axis. Substituting formula (3) and (4) into formula (1), we can obtain the bending vibration equation of the beam.

$\frac{{\partial }^2}{\partial x^2}\left[EI\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial x^2}\right]+\rho A\left(x\right)\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}=f\left(x,t\right).$ (5)

If the beam is of equal section, it becomes

$EI\frac{{\partial }^{4}u\left(x,t\right)}{\partial x^4}+\rho A\frac{{\partial }^{2}u\left(x,t\right)}{\partial t^2}=f\left(x,t\right).$ (6)

Since Winkler first studied the characteristics of foundations, viscoelastic foundation beams have been widely concerned in science and engineering. He built a model of the foundation, which simulates the elastic factor by assuming infinitely close continuous linear springs. The image of Winkler model is shown in Figure 1.

The vibration equation of the viscoelastic Winkler foundation beam can be obtained as follows:

Figure 1. The model of viscoelastic Winkler foundation.

$EI\frac{{\partial }^{4}u}{\partial x^4}+\rho A\frac{{\partial }^{2}u}{\partial t^2}+\mu \frac{\partial u}{\partial t}+ku=f\left(x,t\right),x\in \left(0,L\right),t\in \left(0,T\right],$ (7)

where $\rho$ is the density of beam material, $A$ is the cross sectional area, $\mu$ is the viscosity coefficient, $k$ is the elastic coefficient, $L$ is the length of the beam, and $T$ is the length of the time interval.

However, the Winkler model doesn’t do well in approximating the mechanical behavior of real foundation, which mainly fails to take into account the continuity or cohesion of the foundation. And it ignores the influence of foundation acting on the side of the beam. Pasternak proposed a new model that reconsiders the shear effect on the basis of Winkler model. Wang and Stephens [15] studied the natural frequencies of transverse vibration under a variety of classical boundary conditions by means of variable separation method, which concluded that the frequency of Pasternak-type foundation beams affected by shear is greater than that of Winkler-type foundation beams. Yu et al. [16] gave the dynamic response of an infinite beam on Pasternak foundation under arbitrary dynamic load in the form of an analytical solution. Cai et al. [17] proposed a fractional Pasternak-type foundation model to characterize the time-dependent properties. Miao et al. [18] obtained the analytical solution for the dynamic response of an infinite beam supported on Pasternak foundation under inclined travelling loads. All these results show that damping factors and shear factors have obvious effects on the transverse response of beams. Therefore, we introduced the Pasternak foundation model which is based on the comprehensive consideration of elasticity, damping and shear. The image of Pasternak model is shown in Figure 2.

Figure 2. The model of viscoelastic Pasternak foundation.

The vibration equation of the viscoelastic Pasternak foundation beam can be obtained as follows:

$EI\frac{{\partial }^{4}u}{\partial x^4}+\rho A\frac{{\partial }^{2}u}{\partial t^2}+\mu \frac{\partial u}{\partial t}+ku-G_p\frac{{\partial }^{2}u}{\partial x^2}=f\left(x,t\right),x\in \left(0,L\right),t\in \left(0,T\right],$ (8)

where $G_p$ is the shear coefficient.

In this paper, the numerical method of equation (8) is discussed. And it is organized as follows: In section 2, we use Hermite finite element method to deal with the vibration equation of viscoelastic Pasternak foundation beam. Semi-discrete scheme and fully discrete scheme are proposed, and the convergence order is proved respectively. In section 3, we present four numerical examples to verify the conclusions of the previous theoretical analysis. In section 4, the main conclusions of this paper are summarized.

2. Hermite Finite Element Scheme

Consider the following vibration equation of viscoelastic Pasternak foundation beam

$EI\frac{{\partial }^{4}u}{\partial x^4}+\rho A\frac{{\partial }^{2}u}{\partial t^2}+\mu \frac{\partial u}{\partial t}+ku-G_p\frac{{\partial }^{2}u}{\partial x^2}=f\left(x,t\right),x\in \left(0,L\right),t\in \left(0,T\right].$ (9)

Initial value conditions are given:

$u\left(x,0\right)=\phi \left(x\right),\frac{\partial u}{\partial t}\left(x,0\right)=\psi \left(x\right),x\in \left[0,L\right].$ (10)

The boundary value conditions about beam structure with fixed supports at both ends can be written as

$u\left(0,t\right)=u\left(L,t\right)=0,\frac{\partial u}{\partial x}\left(0,t\right)=\frac{\partial u}{\partial x}\left(L,t\right)=0,t\in \left[0,T\right],$ (11)

where $u\left(x,t\right)$ is the deflection, $\phi \left(x\right)$ , $\psi \left(x\right)$ and $f\left(x,t\right)$ are known smooth functions. The physical meanings of the remaining parameters are given in Section 1.

When the finite element method is used to solve these equations, the time variable is usually taken as the parameter. We let $I=\left[0,L\right]$ and introduce the Sobolev space $H_0^2\left(I\right)=\left\{u|u\in H^2\left(I\right),u\left(0,t\right)=u\left(L,t\right)=0,u_x\left(0,t\right)=u_x\left(L,t\right)=0\right\}$ . Multiplying both sides of Equation (9) by $v\left(x\right)\in H_0^2\left(I\right)$ and integrating, we have

$\left(EI\frac{{\partial }^{4}u}{\partial x^4},v\right)+\left(\rho A\frac{{\partial }^{2}u}{\partial t^2},v\right)+\left(\mu \frac{\partial u}{\partial t},v\right)+\left(ku,v\right)-\left(G_p\frac{{\partial }^{2}u}{\partial x^2},v\right)=\left(f,v\right),$ (12)

where $\left(\cdot ,\cdot \right)$ is the $L^2\left(\text{Ω}\right)$ inner product. Using Green’s formula and the initial-boundary value conditions, we can obtain the variational form: $\forall v\in H_0^2\left(I\right)$ , find $u\in H_0^2\left(I\right)$ to satisfy

$\{\begin{array}{l}\left(EI{u}_{xx},v_{xx}\right)+\left(\rho A{u}_{tt},v\right)+\left(\mu u_t,v\right)+\left(ku,v\right)+\left(G_p{u}_x,v_x\right)=\left(f,v\right).\\ u\left(x,0\right)=\phi \left(x\right),u_t\left(x,0\right)=\psi \left(x\right).\end{array}$ (13)

2.1. Semi-Discrete Scheme and Error Analysis

The finite element semi-discrete scheme is derived below. Let $I_h:0=x_0<x_1<\cdots <x_M=L$ be a uniform subdivision of interval $\left[0,L\right]$ , $h=\frac{L}{M}$ , $x_j=jh$ , $j=0,1,2,\cdots ,M$ . The finite element space is

$V_h=\left\{v\in H_0^2\left(I\right):v_h\in P_3\right\},$ (14)

where $P_3$ is the piecewise cubic Hermite polynomial on $I_h$ . The Hermite interpolation method can satisfy the global continuity of the first derivative of u, which corresponds to two basis functions at each node. One is

${\phi }_0^{\left(0\right)}\left(x\right)={\left(1-\frac{x-x_0}{h_1}\right)}^2\left(2\frac{x-x_0}{h_1}+1\right),\text{if}{x}_0\le x<x_1,$

${\phi }_i^{\left(0\right)}\left(x\right)=\{\begin{array}{ll}{\left(1-\frac{x_i-x}{h_i}\right)}^2\left(2\frac{x_i-x}{h_i}+1\right),\hfill & \text{if}{x}_{i-1}\le x<x_i,\hfill \\ {\left(1-\frac{x-x_i}{h_{i+1}}\right)}^2\left(2\frac{x-x_i}{h_{i+1}}+1\right),\hfill & \text{if}{x}_i\le x<x_{i+1},\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$

${\phi }_M^{\left(0\right)}\left(x\right)={\left(1-\frac{x_M-x}{h_M}\right)}^2\left(2\frac{x_M-x}{h_M}+1\right),\text{if}{x}_{M-1}\le x<x_M.$

The other is

${\phi }_0^{\left(1\right)}\left(x\right)=\left(x-x_0\right){\left(\frac{x-x_0}{h_1}-1\right)}^2,\text{if}{x}_0\le x<x_1,$

${\phi }_i^{\left(1\right)}\left(x\right)=\{\begin{array}{ll}\left(x-x_i\right){\left(\frac{x_i-x}{h_i}-1\right)}^2,\hfill & \text{if}{x}_{i-1}\le x<x_i,\hfill \\ \left(x-x_i\right){\left(\frac{x-x_i}{h_{i+1}}-1\right)}^2,\hfill & \text{if}{x}_i\le x<x_{i+1},\hfill \\ 0,\hfill & \text{otherwise},\hfill \end{array}$

${\phi }_M^{\left(1\right)}\left(x\right)=\left(x-x_M\right){\left(\frac{x_M-x}{h_M}-1\right)}^2,\text{if}{x}_{M-1}\le x<x_M,$

where $i=1,2,3,\cdots ,M-1$ . The finite element space $v_h$ is constructed on the basis functions of all ${\phi }^{\left(0\right)}\left(x\right)$ and ${\phi }^{\left(1\right)}\left(x\right)$ .

We propose a finite element semi-discrete scheme: $\forall v_h\in V_h$ , find $u_h\in V_h$ to satisfy

$\left(EI{u}_{h,xx},v_{h,xx}\right)+\left(\rho A{u}_{h,tt},v_h\right)+\left(\mu u_{h,t},v_h\right)+\left(k{u}_h,v_h\right)+\left(G_p{u}_{h,x},v_{h,x}\right)=\left(f,v_h\right).$ (15)

We define bilinear form as follows:

$a\left(u,v\right)={\displaystyle {\int }_0^L\left(EI\frac{{\partial }^{2}u}{\partial x^2}\frac{{\partial }^{2}v}{\partial x^2}+kuv+G_p\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}\right)\text{d}x}.$ (16)

We let $u_i=u_h\left(x_i,t\right)$ , ${u^{\prime }}_i=\frac{\partial u_h\left(x_i,t\right)}{\partial x}$ , and $u_h$ can be expressed as:

$u_h={\displaystyle \sum }_{i=0}^M\left[u_i{\phi }_i^{\left(0\right)}\left(x\right)+{u^{\prime }}_i{\phi }_i^{\left(1\right)}\left(x\right)\right].$ (17)

Substituting it into the finite element semi-discrete scheme (15) and letting $v_h={\phi }_j^{\left(l\right)}\left(x\right)$ , $\left(l=0,1\right)$ , we have

$\begin{array}{l}{\displaystyle \sum }_{i=0}^M\{u_{i}a\left({\phi }_i^{\left(0\right)},{\phi }_j^{\left(l\right)}\right)+{u^{\prime }}_{i}a\left({\phi }_i^{\left(1\right)},{\phi }_j^{\left(l\right)}\right)+\frac{{\partial }^2{u}_i}{\partial t^2}\left(\rho A{\phi }_i^{\left(0\right)},{\phi }_j^{\left(l\right)}\right)+\frac{{\partial }^2{u^{\prime }}_i}{\partial t^2}\left(\rho A{\phi }_i^{\left(1\right)},{\phi }_j^{\left(l\right)}\right)\\ +\frac{\partial u_i}{\partial t}\left(\mu {\phi }_i^{\left(0\right)},{\phi }_j^{\left(l\right)}\right)+\frac{\partial {u^{\prime }}_i}{\partial t}\left(\mu {\phi }_i^{\left(1\right)},{\phi }_j^{\left(l\right)}\right)\}=\left(f,{\phi }_j^{\left(l\right)}\right),j=0,1,2,\cdots ,M\end{array}$ (18)

According to the boundary conditions that $u_0=0,u_M=0$ , we can define

$U={\left({u^{\prime }}_0,u_1,{u^{\prime }}_1,u_2,{u^{\prime }}_2,\cdots ,u_{M-1},{u^{\prime }}_{M-1},{u^{\prime }}_M\right)}^{\text{T}}.$ (19)

$F={\left[\left(f,{\phi }_0^{\left(1\right)}\right),\left(f,{\phi }_1^{\left(0\right)}\right),\left(f,{\phi }_1^{\left(1\right)}\right),\cdots ,\left(f,{\phi }_{M-1}^{\left(0\right)}\right),\left(f,{\phi }_{M-1}^{\left(1\right)}\right),\left(f,{\phi }_M^{\left(1\right)}\right)\right]}^{\text{T}}.$ (20)

Then Equation (18) can be rewritten in the following form

$A\frac{{\partial }^{2}U}{\partial t^2}+B\frac{\partial U}{\partial t}+CU=F,$ (21)

where $A,B,C$ are all seven diagonal matrices. Then we can use numerical methods to solve the linear equations and get numerical solutions.

To estimate the error, we define $Du=\frac{\partial u}{\partial x}$ and introduce the biharmonic projection $R_h:H_0^2\left(I\right)\to V_h$ which satisfies

$EI\left(D^2\left(u-R_{h}u\right),D^2{v}_h\right)+k\left(\left(u-R_{h}u\right),v_h\right)+G_p\left(D\left(u-R_{h}u\right),D{v}_h\right)=0.$ (22)

The projection has the following estimation [19].

Lemma 2.1 $k$ ($1\le k\le 3$ ) is the degree of piecewise polynomial, and for any $u\in H_0^2\cap H^{k+1}$ , we have the following conclusion:

$\Vert u-R_{h}u\Vert +h{\Vert u-R_{h}u\Vert }_1+h^2{\Vert u-R_{h}u\Vert }_2\le C{h}^{k+1}{\Vert u\Vert }_{k+1}.$ (23)

Here, we use the standard notation $W^{m,q}\left(\text{Ω}\right)$ for Sobolev space on $\text{Ω}$ with norm ${\Vert \cdot \Vert }_{m,q}$ . When $q=2$ , we denote $H^m\left(\text{Ω}\right)=W^{m,2}\left(\text{Ω}\right)$ , and in the format of norm it is ${\Vert \cdot \Vert }_m={\Vert \cdot \Vert }_{m,2}$ . For $m=0$ , we denote ${\Vert \cdot \Vert }_0=\Vert \cdot \Vert$ . $C$ is a general constant. We are now ready to prove the error estimation.

Theorem 2.1 Assuming that $u$ and $u_h$ are solutions to equation (13) and (15) respectively, , we have

$\Vert u-u_h\Vert \le C{h}^4\left[{\Vert u\Vert }_4+{\left({\displaystyle {\int }_0^t\left({\Vert u_{tt}\Vert }_4^2+{\Vert u_t\Vert }_4^2\right)\text{d}s}\right)}^{\frac{1}{2}}\right],$ (24)

where $C$ is a constant, independent of $h$ .

Proof. From Lemma 2.1, we can get

$\Vert u-R_{h}u\Vert \le C{h}^4{\Vert u\Vert }_4.$ (25)

Letting $v=v_h$ in the variational form (13) and subtracting (15) from (13), we have

$\begin{array}{l}\left(EI\left(u_{xx}-u_{h,xx}\right),v_{h,xx}\right)+\left(\rho A\left(u_{tt}-u_{h,tt}\right),v_h\right)\\ +\left(\mu \left(u_t-u_{h,t}\right),v_h\right)+\left(k\left(u-u_h\right),v_h\right)+\left(G_p\left(u_x-u_{h,x}\right),v_{h,x}\right)=0.\end{array}$ (26)

We introduce the following notation

$u-u_h=\left(u-R_{h}u\right)+\left(R_{h}u-u_h\right)=\eta +\theta .$ (27)

Using the projection $R_h$ in (22), we can rewrite (26) as

$\begin{array}{l}\left(EI{\theta }_{xx},v_{h,xx}\right)+\left(\rho A{\theta }_{tt},v_h\right)+\left(\mu {\theta }_t,v_h\right)+\left(k\theta ,v_h\right)+\left(G_p{\theta }_x,v_{h,x}\right)\\ =-\left(\rho A{\eta }_{tt},v_h\right)-\left(\mu {\eta }_t,v_h\right).\end{array}$ (28)

Letting $v_h={\theta }_t$ , we get

$\begin{array}{l}\left(EI{\theta }_{xx},{\theta }_{t,xx}\right)+\left(\rho A{\theta }_{tt},{\theta }_t\right)+\left(\mu {\theta }_t,{\theta }_t\right)+\left(k\theta ,{\theta }_t\right)+\left(G_p{\theta }_x,{\theta }_{t,x}\right)\\ =-\left(\rho A{\eta }_{tt},{\theta }_t\right)-\left(\mu {\eta }_t,{\theta }_t\right).\end{array}$ (29)

Using Cauchy-Schwarz inequality and $\epsilon$ -inequality, we obtain

(30)

According to the finite element method and the theory of biharmonic projection, we can suppose that $u_h^0=R_{h}u\left(t_0\right)$ . Integrating the above formula from 0 to $t$ with respect to time, we get

$EI{\Vert {\theta }_{xx}\Vert }^2+\rho A{\Vert {\theta }_t\Vert }^2+k{\Vert \theta \Vert }^2+G_p{\Vert {\theta }_x\Vert }^2\le {\displaystyle {\int }_0^t\frac{{\rho }^2{A}^2}{\mu }{\Vert {\eta }_{tt}\Vert }^2\text{d}s}+{\displaystyle {\int }_0^t\mu {\Vert {\eta }_t\Vert }^2\text{d}s}.$ (31)

According to Lemma 2.1, the following estimates can be obtained

$\Vert {\eta }_{tt}\Vert =\Vert u_{tt}-R_h{u}_{tt}\Vert \le C{h}^4{\Vert u_{tt}\Vert }_4.$ (32)

$\Vert {\eta }_t\Vert =\Vert u_t-R_h{u}_t\Vert \le C{h}^4{\Vert u_t\Vert }_4.$ (33)

Substituting them into (31) and after simplification, we have

$\Vert \theta \Vert \le C{h}^4{\left({\displaystyle {\int }_0^t\left({\Vert u_{tt}\Vert }_4^2+{\Vert u_t\Vert }_4^2\right)\text{d}s}\right)}^{\frac{1}{2}}.$ (34)

Using triangle inequality, we can get

$\Vert u-u_h\Vert \le \Vert \eta \Vert +\left|\theta \right|\le C{h}^4\left[{\Vert u\Vert }_4+{\left({\displaystyle {\int }_0^t\left({\Vert u_{tt}\Vert }_4^2+{\Vert u_t\Vert }_4^2\right)\text{d}s}\right)}^{\frac{1}{2}}\right].$ (35)

We complete the proof of this theorem.

2.2. Fully Discrete Scheme and Error Analysis

In this section, we propose a fully discrete Hermite finite element approximation scheme for the vibration equation of Pasternak-type viscoelastic foundation beam. The error convergence order is obtained and proved.

First, we consider that $0=t_0<t_1<\cdots <t_N=T$ is a uniform partition over the interval $\left[0,T\right]$ , $\tau =\frac{T}{N}$ , $t_n=n\tau$ , $n=0,1,2,\cdots ,N$ . Then, partial derivatives with

respect to time are approximated using the central difference scheme. So we can get the fully discrete scheme: $\forall v_h\in V_h$ , find $u_h^n\in V_h$ to satisfy

$\begin{array}{l}\left(EI\frac{u_{h,xx}^{n+1}+u_{h,xx}^{n-1}}{2},v_{h,xx}\right)+\left(\rho A\frac{u_h^{n+1}-2u_h^n+u_h^{n-1}}{{\tau }^2},v_h\right)+\left(\mu \frac{u_h^{n+1}-u_h^{n-1}}{2\tau },v_h\right)\\ +\left(k\frac{u_h^{n+1}+u_h^{n-1}}{2},v_h\right)+\left(G_p\frac{u_{h,x}^{n+1}+u_{h,x}^{n-1}}{2},v_{h,x}\right)=\left(f^n,v_h\right).\end{array}$ (36)

Now we are ready to present and prove the main convergence theorem.

Theorem 2.2 Assuming that $u^n$ and $u_h^n$ are solutions to equation (13) and (36) respectively, $u\in H_0^2\cap H^4$ , we have

$\Vert u^n-u_h^n\Vert \le C\left({\tau }^2+h^4\right),$ (37)

where $C$ is a constant, independent of $\tau$ and $h$ .

Proof. Choosing $v=v_h$ in (13) and subtracting (36) from (13), we get the error equation

$\begin{array}{l}\left(EI\left(u_{xx}^n-\frac{u_{h,xx}^{n+1}+u_{h,xx}^{n-1}}{2}\right),v_{h,xx}\right)+\left(\rho A\left(u_{tt}^n-\frac{u_h^{n+1}-2u_h^n+u_h^{n-1}}{{\tau }^2}\right),v_h\right)\\ +\left(\mu \left(u_t^n-\frac{u_h^{n+1}-u_h^{n-1}}{2\tau }\right),v_h\right)+\left(k\left(u^n-\frac{u_h^{n+1}+u_h^{n-1}}{2}\right),v_h\right)\\ +\left(G_p\left(u_x^n-\frac{u_{h,x}^{n+1}+u_{h,x}^{n-1}}{2}\right),v_{h,x}\right)=0.\end{array}$ (38)

According to the notation of (27), we can rewrite the error Equation (38) in the following equivalent form.

$\begin{array}{l}\left(EI\frac{{\theta }_{xx}^{n+1}+{\theta }_{xx}^{n-1}}{2},v_{h,xx}\right)+\left(\rho A\frac{{\theta }^{n+1}-2{\theta }^n+{\theta }^{n-1}}{{\tau }^2},v_h\right)+\left(\mu \frac{{\theta }^{n+1}-{\theta }^{n-1}}{2\tau },v_h\right)\\ +\left(k\frac{{\theta }^{n+1}+{\theta }^{n-1}}{2},v_h\right)+\left(G_p\frac{{\theta }_x^{n+1}+{\theta }_x^{n-1}}{2},v_{h,x}\right)\\ =\left(EI\left(\frac{u_{xx}^{n+1}+u_{xx}^{n-1}}{2}-u_{xx}^n\right),v_{h,xx}\right)-\left(EI\frac{{\eta }_{xx}^{n+1}+{\eta }_{xx}^{n-1}}{2},v_{h,xx}\right)\\ +\left(\rho A\left(\frac{u^{n+1}-2u^n+u^{n-1}}{{\tau }^2}-u_{tt}^n\right),v_h\right)-\left(\rho A\frac{{\eta }^{n+1}-2{\eta }^n+{\eta }^{n-1}}{{\tau }^2},v_h\right)\\ +\left(\mu \left(\frac{u^{n+1}-u^{n-1}}{2\tau }-u_t^n\right),v_h\right)-\left(\mu \frac{{\eta }^{n+1}-{\eta }^{n-1}}{2\tau },v_h\right)\\ +\left(k\left(\frac{u^{n+1}+u^{n-1}}{2}-u^n\right),v_h\right)-\left(k\frac{{\eta }^{n+1}+{\eta }^{n-1}}{2},v_h\right)\\ +\left(G_p\left(\frac{u_x^{n+1}+u_x^{n-1}}{2}-u_x^n\right),v_{h,x}\right)-\left(G_p\frac{{\eta }_x^{n+1}+{\eta }_x^{n-1}}{2},v_{h,x}\right).\end{array}$ (39)

Using the projection $R_h$ in (22), we can obtain

$\left(EI\frac{{\eta }_{xx}^{n+1}+{\eta }_{xx}^{n-1}}{2},v_{h,xx}\right)+\left(k\frac{{\eta }^{n+1}+{\eta }^{n-1}}{2},v_h\right)+\left(G_p\frac{{\eta }_x^{n+1}+{\eta }_x^{n-1}}{2},v_{h,x}\right)=0.$ (40)

Substitute (40) into (39) and let $v_h=\frac{{\theta }^{n+1}+{\theta }^{n-1}}{2}$ . Using Cauchy-Schwarz inequality and $\epsilon$ -inequality, we obtain

$\begin{array}{l}\frac{\rho A}{2}\left({\Vert \frac{{\theta }^{n+1}-{\theta }^n}{\tau }\Vert }^2-{\Vert \frac{{\theta }^n-{\theta }^{n-1}}{\tau }\Vert }^2\right)+\frac{\mu }{4\tau }\left({\Vert {\theta }^{n+1}\Vert }^2-{\Vert {\theta }^{n-1}\Vert }^2\right)\\ \le \frac{2{\rho }^2{A}^2}{k}{\Vert \frac{{\eta }^{n+1}-2{\eta }^n+{\eta }^{n-1}}{{\tau }^2}\Vert }^2+\frac{2{\rho }^2{A}^2}{k}{\Vert \frac{u^{n+1}-2u^n+u^{n-1}}{{\tau }^2}-u_{tt}^n\Vert }^2\\ +\frac{2{\mu }^2}{k}{\Vert \frac{{\eta }^{n+1}-{\eta }^{n-1}}{2\tau }\Vert }^2+\frac{2{\mu }^2}{k}{\Vert \frac{u^{n+1}-u^{n-1}}{2\tau }-u_t^n\Vert }^2+\frac{k}{2}{\Vert \frac{u^{n+1}+u^{n-1}}{2}-u^n\Vert }^2\\ +\frac{EI}{4}{\Vert \frac{u_{xx}^{n+1}+u_{xx}^{n-1}}{2}-u_{xx}^n\Vert }^2+\frac{G_p}{4}{\Vert \frac{u_x^{n+1}+u_x^{n-1}}{2}-u_x^n\Vert }^2.\end{array}$ (41)

According to Taylor expansion formula and Lemma 2.1, we have the following proved estimations.

${\Vert \frac{{\eta }^{n+1}-2{\eta }^n+{\eta }^{n-1}}{{\tau }^2}\Vert }^2\le \frac{C{h}^8}{\tau }{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_{tt}\Vert }_4^2\text{d}s}$ (42)

${\Vert \frac{{\eta }^{n+1}-{\eta }^{n-1}}{2\tau }\Vert }^2\le \frac{C{h}^8}{\tau }{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_t\Vert }_4^2\text{d}s}$ (43)

${\Vert \frac{u^{n+1}-2u^n+u^{n-1}}{{\tau }^2}-u_{tt}^n\Vert }^2\le \frac{{\tau }^3}{126}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \frac{{\partial }^{4}u}{\partial t^4}\Vert }^2\text{d}s}$ (44)

${\Vert \frac{u^{n+1}-u^{n-1}}{2\tau }-u_t^n\Vert }^2\le \frac{{\tau }^3}{80}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \frac{{\partial }^{3}u}{\partial t^3}\Vert }^2\text{d}s}$ (45)

${\Vert \frac{u^{n+1}+u^{n-1}}{2}-u^n\Vert }^2\le \frac{{\tau }^3}{6}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_{tt}\Vert }^2\text{d}s}$ (46)

${\Vert \frac{u_{xx}^{n+1}+u_{xx}^{n-1}}{2}-u_{xx}^n\Vert }^2\le \frac{{\tau }^3}{6}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \frac{{\partial }^{4}u}{\partial x^2\partial t^2}\Vert }^2\text{d}s}$ (47)

${\Vert \frac{u_x^{n+1}+u_x^{n-1}}{2}-u_x^n\Vert }^2\le \frac{{\tau }^3}{6}{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \frac{{\partial }^{3}u}{\partial x\partial t^2}\Vert }^2\text{d}s}$ (48)

Substituting them into (41), we get

$\begin{array}{l}\frac{\rho A}{2}\left({\Vert \frac{{\theta }^{n+1}-{\theta }^n}{\tau }\Vert }^2-{\Vert \frac{{\theta }^n-{\theta }^{n-1}}{\tau }\Vert }^2\right)+\frac{\mu }{4\tau }\left({\Vert {\theta }^{n+1}\Vert }^2-{\Vert {\theta }^{n-1}\Vert }^2\right)\\ \le \frac{C{h}^8}{\tau }{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_{tt}\Vert }_4^2\text{d}s}+C{\tau }^3{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \frac{{\partial }^{4}u}{\partial t^4}\Vert }^2\text{d}s}+\frac{C{h}^8}{\tau }{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_t\Vert }_4^2\text{d}s}\\ +C{\tau }^3{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \frac{{\partial }^{3}u}{\partial t^3}\Vert }^2\text{d}s}+C{\tau }^3{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert u_{tt}\Vert }^2\text{d}s}+C{\tau }^3{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \frac{{\partial }^{4}u}{\partial x^2\partial t^2}\Vert }^2\text{d}s}\\ +C{\tau }^3{\displaystyle {\int }_{t_{n-1}}^{t_{n+1}}{\Vert \frac{{\partial }^{3}u}{\partial x\partial t^2}\Vert }^2\text{d}s}.\end{array}$ (49)

Multiplying both sides by $\frac{4\tau }{\mu }$ and summing the above formula from 1 to $n$ , we obtain

$\begin{array}{l}\frac{2\rho A\tau }{\mu }\left({\Vert \frac{{\theta }^{n+1}-{\theta }^n}{\tau }\Vert }^2-{\Vert \frac{{\theta }^1-{\theta }^0}{\tau }\Vert }^2\right)+\left({\Vert {\theta }^{n+1}\Vert }^2+{\Vert {\theta }^n\Vert }^2-{\Vert {\theta }^1\Vert }^2-{\Vert {\theta }^0\Vert }^2\right)\\ \le C{h}^8{\displaystyle {\int }_{t_0}^{t_{n+1}}{\Vert u_{tt}\Vert }_4^2\text{d}s}+C{\tau }^4{\displaystyle {\int }_{t_0}^{t_{n+1}}{\Vert \frac{{\partial }^{4}u}{\partial t^4}\Vert }^2\text{d}s}+C{h}^8{\displaystyle {\int }_{t_0}^{t_{n+1}}{\Vert u_t\Vert }_4^2\text{d}s}\\ +C{\tau }^4{\displaystyle {\int }_{t_0}^{t_{n+1}}{\Vert \frac{{\partial }^{3}u}{\partial t^3}\Vert }^2\text{d}s}+C{\tau }^4{\displaystyle {\int }_{t_0}^{t_{n+1}}{\Vert u_{tt}\Vert }^2\text{d}s}+C{\tau }^4{\displaystyle {\int }_{t_0}^{t_{n+1}}{\Vert \frac{{\partial }^{4}u}{\partial x^2\partial t^2}\Vert }^2\text{d}s}\\ +C{\tau }^4{\displaystyle {\int }_{t_0}^{t_{n+1}}{\Vert \frac{{\partial }^{3}u}{\partial x\partial t^2}\Vert }^2\text{d}s}.\end{array}$ (50)

According to the finite element method and the theory of biharmonic projection, we can suppose that $u_h^0=R_{h}u\left(t_0\right)$ and $u_h^1=R_{h}u\left(t_1\right)$ . So ${\theta }^0=0$ and ${\theta }^1=0$ . (50) can be rewritten as

${\Vert {\theta }^n\Vert }^2\le C\left({\tau }^4+h^8\right).$ (51)

So

$\Vert {\theta }^n\Vert \le C\left({\tau }^2+h^4\right).$ (52)

Using triangle inequality and Lemma 2.1, we get

$\Vert u^n-u_h^n\Vert =\Vert {\eta }^n+{\theta }^n\Vert \le \Vert {\eta }^n\Vert +\Vert {\theta }^n\Vert \le C\left({\tau }^2+h^4\right).$ (53)

We complete the proof.

3. Numerical Example

In this section, we present four examples. The first example tests the validity of Hermite finite element method. The last three examples explore the influence of viscoelastic Pasternak foundation parameters on beam vibration.

Example 1 Consider the following initial-boundary value problem of Euler-Bernoulli beam on viscoelastic Pasternak foundation

$\{\begin{array}{ll}EI{u}_{xxxx}+\rho A{u}_{tt}+\mu u_t+ku-G_p{u}_{xx}=f\left(x,t\right),\hfill & \left(x,t\right)\in \left[0,1\right]\times \left(0,1\right],\hfill \\ u\left(x,0\right)=0,u_t\left(x,0\right)=0,\hfill & x\in \left[0,1\right],\hfill \\ u\left(0,t\right)=u\left(1,t\right)=0,u_x\left(0,t\right)=u_x\left(1,t\right)=0,\hfill & t\in \left[0,1\right].\hfill \end{array}$ (54)

We choose $u\left(x,t\right)=t^2\left(1-\cos \left(2\pi x\right)\right)$ as the exact solution, which satisfies the initial-boundary value condition. Table 1 lists the physical parameters of Euler-Bernoulli beams and viscoelastic Pasternak foundation.

For convenience of calculation, we process the parameters to get dimensionless quantity and substitute them into the equation. First of all, we fix time step $\tau =1/10000$ and change space step to obtain the error results under different norms in Table 2. The results show that the space convergence order in the sense of $L^2$ -norm is approximately 4, which is consistent with the theoretical results in Theorem 2.1.

Table 1. The physical parameters in the equation.

Name	Symbol	Value	Unit
Young modulus	$E$	210	GPa
Moment of inertia	$I$	3.055 × 10−5	m4
Density	$\rho$	7850	kg/m3
Cross sectional area	$A$	7.69 × 10−3	m2
Elastic coefficient	$k$	3.5 × 107	N/m2
Viscosity coefficient	$\mu$	1732.5	N∙s/m2
Shear coefficient	$G_p$	5.92 × 108	N

Table 2. $L^2$ -norm, $H^1$ -norm and $H^2$ -norm error and space convergence order.

$h$	$\tau$	${\Vert u-u_h\Vert }_0$	Order	${\Vert u-u_h\Vert }_1$	Order	${\Vert u-u_h\Vert }_2$	Order
$\frac{1}{2^3}$	$\frac{1}{10000}$	2.9201 × 10−4	-	1.1324 × 10−2	-	6.7444 × 10−1	-
$\frac{1}{2^4}$	$\frac{1}{10000}$	2.0402 × 10−5	3.8392	1.5140 × 10−3	2.9030	1.6185 × 10−1	2.0590
$\frac{1}{2^5}$	$\frac{1}{10000}$	1.3452 × 10−6	3.9229	1.9300 × 10−4	2.9716	4.0156 × 10−2	2.0110
$\frac{1}{2^6}$	$\frac{1}{10000}$	8.0162 × 10−8	4.0687	2.4451 × 10−5	2.9807	1.0121 × 10−2	1.9882

And then, we fix space step $h=1/1000$ and change time step to obtain the error results under different norms in Table 3. The results show that the time convergence order in the sense of $L^2$ -norm is approximately 2, which is consistent with the theoretical results in Theorem 2.2.

Table 3. $L^{\infty }$ -norm and $L^2$ -norm error and time convergence order.

$h$	$\tau$	${\Vert u-u_h\Vert }_{\infty }$	Order	${\Vert u-u_h\Vert }_0$	Order
$\frac{1}{1000}$	$\frac{1}{2^3}$	5.4539 × 10−2	-	3.6863 × 10−2	-
$\frac{1}{1000}$	$\frac{1}{2^4}$	1.2855 × 10−2	2.0849	9.2058 × 10−3	2.0016
$\frac{1}{1000}$	$\frac{1}{2^5}$	3.1744 × 10−3	2.0178	2.2935 × 10−3	2.0050
$\frac{1}{1000}$	$\frac{1}{2^6}$	7.9838 × 10−4	1.9913	5.7227 × 10−4	2.0028
$\frac{1}{1000}$	$\frac{1}{2^7}$	1.9973 × 10−4	1.9990	1.4294 × 10−4	2.0013

We can obtain the three-dimensional images of numerical solution and exact solution in Figure 3.

Figure 3. 3D images of numerical solution and exact solution of the deflection $u$ .

Example 2 In order to explore the influence of viscosity coefficient of foundation on beam vibration, We let $f\left(x,t\right)=0$ , that is, the beam is assumed to be in a state of free vibration. A suitable initial value is assigned to the beam vibration equation. We change the value of viscosity coefficient and fix the remaining parameters. The vibration image of the beam midpoint is obtained, as shown in Figure 4. It shows that the viscosity coefficient can affect the amplitude, but it does not change the free vibration frequency. The amplitude of beam vibration decreases with the increase of viscosity coefficient.

Figure 4. The deflection at the midpoint of beam corresponding to different viscosity coefficients $\mu$ .

Example 3 In order to explore the influence of elastic coefficient of foundation on beam vibration, we change the elastic coefficient and fix the remaining parameters. The vibration image of the beam midpoint is obtained, as shown in Figure 5. It shows that the vibration frequency increases with the increase of elastic coefficient.

Figure 5. The deflection at the midpoint of beam corresponding to different elastic coefficients $k$ .

Example 4 The biggest difference between Pasternak foundation model and Winkler foundation model is that the former takes shear force into account. To explore the impact of this distinction, we plot the image of beam vibration on Winkler foundation ($G_p=0$ ) and on Pasternak foundation ($G_p=5.92\times {10}^8$ ), respectively, as shown in Figure 6. We can draw a conclusion that the added consideration of shear force leads to greater stiffness between the structure and foundation, resulting in an increase in the free vibration frequency.

Figure 6. The deflection at the midpoint of beam on Winkler foundation and Pasternak foundation.

4. Conclusions

Considering that the beam vibration equation requires the first derivative of deflection $u$ to be continuous on the whole, we choose Hermite finite element method, which has some advantages. Compared with other numerical methods, Hermite finite element method has a higher order of error convergence, which can reach $O\left(h^4\right)$ in the sense of $L^2$ -norm (Theorem 2.1). In addition, it can ensure that the interpolation function has continuous first derivative. Due to the characteristics of Hermite finite element method, when we calculate the deflection $u$ , we can also directly get the $u^{\prime }$ ((19)-(21)), that is, the angle value. This also reduces the error caused by the intermediate process.

In this article, we propose the Hermite finite element scheme to solve the vibration equation of viscoelastic Pasternak foundation beam. Semi-discrete and fully discrete schemes are given and their error estimation in the sense of $L^2$ norm is proved, which has fourth-order convergence accuracy in space and second-order in time. The results of theoretical analysis are verified in numerical examples, and the effects of foundation coefficients on beam vibration are investigated. It can be seen from the experiment that the shear coefficient can affect the free vibration characteristics of the foundation beam, including vibration frequency and amplitude. In engineering practice, it is necessary to focus on the effects of shear forces to ensure the stability and safety of the structure.

Acknowledgments

The work is supported by Natural Science Foundation of Shandong Province (contract grant number: ZR2021MA063). The authors would like to thank the referees and editors for their valuable comments and suggestions.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Senalp, A.D., Arikoglu, A., Ozkol, I. and Dogan, V.Z. (2010) Dynamic Response of a Finite Length Euler-Bernoulli Beam on Linear and Nonlinear Viscoelastic Foundations to a Concentrated Moving Force. Journal of Mechanical Science and Technology, 24, 1957-1961. https://doi.org/10.1007/s12206-010-0704-x
[2]	Babilio, E. (2014) Dynamics of Functionally Graded Beams on Viscoelastic Foundation. International Journal of Structural Stability and Dynamics, 14, Article ID: 1440014. https://doi.org/10.1142/s0219455414400148
[3]	Hörmann, G., Konjik, S. and Oparnica, L. (2013) Generalized Solutions for the Euler–Bernoulli Model with Zener Viscoelastic Foundations and Distributional Forces. Analysis and Applications, 11, Article ID: 1350017. https://doi.org/10.1142/s0219530513500176
[4]	Beskou, N.D. and Muho, E.V. (2018) Dynamic Response of a Finite Beam Resting on a Winkler Foundation to a Load Moving on Its Surface with Variable Speed. Soil Dynamics and Earthquake Engineering, 109, 222-226. https://doi.org/10.1016/j.soildyn.2018.02.033
[5]	Elhuni, H. and Basu, D. (2019) Dynamic Soil Structure Interaction Model for Beams on Viscoelastic Foundations Subjected to Oscillatory and Moving Loads. Computers and Geotechnics, 115, Article ID: 103157. https://doi.org/10.1016/j.compgeo.2019.103157
[6]	Lu, T., Wang, Z. and Liu, D. (2019) Analysis of Complex Modal Characteristics of Fractional Derivative Viscoelastic Rotating Beams. Shock and Vibration, 2019, Article ID: 5715694. https://doi.org/10.1155/2019/5715694
[7]	Snehasagar, G., Krishnanunni, C.G. and Rao, B.N. (2019) Dynamics of Vehicle-Pavement System Based on a Viscoelastic Euler-Bernoulli Beam Model. International Journal of Pavement Engineering, 21, 1669-1682. https://doi.org/10.1080/10298436.2018.1562189
[8]	Praharaj, R.K. and Datta, N. (2020) On the Transient Response of Plates on Fractionally Damped Viscoelastic Foundation. Computational and Applied Mathematics, 39, Article No. 256. https://doi.org/10.1007/s40314-020-01285-6
[9]	Frýba, L., Nakagiri, S. and Yoshikawa, N. (1993) Stochastic Finite Elements for a Beam on a Random Foundation with Uncertain Damping under a Moving Force. Journal of Sound and Vibration, 163, 31-45. https://doi.org/10.1006/jsvi.1993.1146
[10]	Lou, P. (2007) Finite-Element Formulae for Calculating the Sectional Forces of a Bernoulli-Euler Beam on Continuously Viscoelastic Foundation Subjected to Concentrated Moving Loads. Shock and Vibration, 15, 147-162. https://doi.org/10.1155/2008/309216
[11]	Taeprasartsit, S. (2013) Nonlinear Free Vibration of Thin Functionally Graded Beams Using the Finite Element Method. Journal of Vibration and Control, 21, 29-46. https://doi.org/10.1177/1077546313484506
[12]	Sánchez, E.D., Nallim, L.G., Bellomo, F.J. and Oller, S.H. (2020) Generalized Viscoelastic Model for Laminated Beams Using Hierarchical Finite Elements. Composite Structures, 235, Article ID: 111794. https://doi.org/10.1016/j.compstruct.2019.111794
[13]	Wang, T., Jiang, Z., Zhu, A. and Yin, Z. (2022) A Mixed Finite Volume Element Method for Time-Fractional Damping Beam Vibration Problem. Fractal and Fractional, 6, Article No. 523. https://doi.org/10.3390/fractalfract6090523
[14]	Sun, X., Zhu, A., Yin, Z. and Ji, P. (2023) Hermite Finite Element Method for Time Fractional Order Damping Beam Vibration Problem. Fractal and Fractional, 7, Article No. 739. https://doi.org/10.3390/fractalfract7100739
[15]	Wang, T.M. and Stephens, J.E. (1977) Natural Frequencies of Timoshenko Beams on Pasternak Foundations. Journal of Sound and Vibration, 51, 149-155. https://doi.org/10.1016/s0022-460x(77)80029-1
[16]	Yu, H., Cai, C., Yuan, Y. and Jia, M. (2017) Analytical Solutions for Euler-Bernoulli Beam on Pasternak Foundation Subjected to Arbitrary Dynamic Loads. International Journal for Numerical and Analytical Methods in Geomechanics, 41, 1125-1137. https://doi.org/10.1002/nag.2672
[17]	Cai, W., Chen, W. and Xu, W. (2016) Fractional Modeling of Pasternak-Type Viscoelastic Foundation. Mechanics of Time-Dependent Materials, 21, 119-131. https://doi.org/10.1007/s11043-016-9321-0
[18]	Miao, Y., He, H., Yang, Q. and Shi, Y. (2020) Analytical Solution Considering the Tangential Effect for an Infinite Beam on a Viscoelastic Pasternak Foundation. Applied Mathematical Modelling, 85, 231-243. https://doi.org/10.1016/j.apm.2020.03.031
[19]	Ciarlet, P.G. and Oden, J.T. (1978) The Finite Element Method for Elliptic Problems. Journal of Applied Mechanics, 45, 968-969. https://doi.org/10.1115/1.3424474