Study on Stability of Elastic Compression Bending Bar in Viscoelastic Medium

Abstract

In the southeastern coastal regions of China, thick layers of marine soft soil are widely distributed, exhibiting characteristics such as high compressibility, high porosity, low strength, high sensitivity, and easy thixotropy, and these viscoelastic behaviors of foundation soil have significant implications for elastic compression bending bar, as evidenced by issues such as post-construction settlement of roadbeds and long-term operation deviation of bridge pile foundations. In this study, a mechanical model of an elastic bar embedded in an elastic and viscoelastic medium, fixed at the base and free at the top, is established based on the Winkler foundation assumption. The deflection function of a bar subjected to both axial force and locally distributed horizontal load is derived using the Rayleigh-Ritz method. Utilizing the elastic-viscoelastic correspondence principle, the viscoelastic medium surrounding the bar is modeled as an elastic medium in which the ground reaction coefficient varies within phase space formulation. This study provides a robust theoretical foundation for soft soil foundation engineering projects and fills a significant gap in the literature by offering a comprehensive framework for understanding displacement in elastic bars within viscoelastic media. Drawing upon the derivation of the deformation function for elastic rods within a viscoelastic medium, the findings of this research hold significant applicability across a range of domains. These include, but are not limited to, the expansion of roadways in regions characterized by coastal soft soil, as well as the monitoring of deformation and lifespan in bridge pile foundations.

1. Introduction

The macroscopic study of soil rheology involves treating soil as a continuous and homogeneous medium and utilizing mathematical and mechanical methods to derive mathematical expressions that describe soil rheological phenomena based on macroscopic data. While this approach lacks the depth of microscopic studies in exploring the underlying mechanisms of soil rheology, it provides an intuitive depiction of the external characteristics of soil rheology. In particular, it enables the determination of rheological properties for soil samples from specific regions, thereby satisfying practical engineering requirements. As such, macroscopic study has long been a focus of soil rheology research. Common methods employed in macroscopic research include element theory and empirical theory. Element theory, also known as model theory, considers the rheological properties of soil to result from the combined action of elasticity, viscosity, and plasticity. These properties are simulated using spring, dashpot, and friction elements, respectively. Element theory is simple and intuitive, providing a comprehensive reflection of various rheological properties of media and is widely applied. Empirical theory, on the other hand, is based on specific soil qualities and conditions. It involves fitting a constitutive model to soil test results and has found some application in engineering.

In engineering practice, diseases such as post-construction settlement of roadbeds and long-term operation deviation of bridge pile foundations reflect the existence of rheological effects of foundation soil. In southeastern coastal areas of China, there are widely distributed thick layers of marine soft soil, which have characteristics such as high compressibility, high porosity, low strength, high sensitivity, and easy thixotropy. When carrying out roadbed, embankment filling and other earthwork operations in soft soil areas, soil consolidation settlement and long-term rheological effects are larger than those in non-soft soil areas, which will have a greater long-term impact on pile foundations of buildings or bridges in nearby areas. Therefore, when studying the impact of roadbed filling near existing bridges in soft soil areas on bridge piles, it is necessary to consider the long-term impact of soft soil rheology on bridge pile internal force state and geometric shape changes.

Notable examples of macroscopic research include the work of Taylor and Merchant [1], who considered the rheological properties of soil in their analysis of soil consolidation and employed the Kelvin model to describe the deformation of the soil skeleton. Tan [2,3] treated clay as a Maxwell body and, building on this theoretical foundation, studied the secondary time effect in the consolidation process of soft soil, proposing the Chen model. Gibson [4] used the Merchant model to fit soil rheological properties and conducted in-depth research on one-dimensional rheological consolidation problems.

In the field of engineering structures, a bar subjected to both axial force and lateral load is referred to as a compression-bending bar. Up to now, extensive investigations have been carried out into the influence and design of lateral braces on columns and beams, including the bracing stiffness, bracing strength requirements, and the corresponding buckling loads of members [5,6,7,8,9,10,11]. Winter [12] theoretically and experimentally developed the principles of minimum stiffness and minimum strength requirements for equally spaced braces by introducing a simplified rigid link model with fictitious hinges at the bracing points, which became the foundation of most of the practical bracing design formulations for braced columns. Timoshenko and Gere [13] derived the buckling loads and obtained the same values of full bracing stiffness as those by Winter for perfect straight columns, using a theoretical derivation based on the segment equilibrium and deformation compatibility. A typical example of such a bar is a pile foundation driven into soil, which often bears additional loads transmitted from the impact of pile loading construction in adjacent areas while simultaneously bearing axial loads. Several researchers have conducted forward-looking studies in this field and have made notable progress [14,15,16,17]. Dai [18] derived equations for calculating the internal force and deformation of horizontally loaded piles using the elastic foundation beam method in combination with the finite difference principle, matlock method (m method), two-parameter method, and p-y curve method. This research demonstrates the maturity of the finite difference method in its application to horizontally loaded piles.

Reese and Welch [19] established differential equations for the force of foundation piles under inclined eccentric load as early as 1975. Sastry et al. [20,21,22] later investigated the ultimate bearing capacity and force characteristics of piles under inclined loads in homogeneous sandy soils, homogeneous clay soils, and stratified soils through indoor model tests. Zhao [23] conducted an in-depth study of the ultimate bearing capacity of foundation piles under inclined load based on the principle of joint action of pile and soil and elasticity theory, using numerical methods and experimental means. By experiments and numerical simulations, some researchers have investigated the flexural behavior of the pile under combined axial and lateral loads. ZHUKOV and BALOV [24] studied pile behavior by in-situ tests, without considering the influence of the pile length above ground. ZHANG et al. [25] deduced the elastoplastic solutions for partially embedded piles, in which the coefficient of subgrade reaction was supposed to be a constant. By assuming that the coefficient of subgrade reaction increases linearly with depth, HAN and FROST [26] obtained the variational solutions for the transversely isotropic pile. The subgrade reaction method is widely used because of its simplicity and reasonable accuracy. However, the coefficient of subgrade reaction is not a fundamental soil property but a model parameter, and the soil resistance is modeled with independent springs without accounting for the continuity of soil [27]. POULOS and DAVIS [28] put forward the elasticity approach for single piles under lateral loads, in which the soil was assumed to be an ideal, homogeneous, isotropic, semi-infinite elastic material, and the differential solutions were obtained.

The research proposed a simplified analysis method for power series calculation based on the m method. However, this method has a large number of calculation parameters and its formulae are lengthy and complex, limiting its application to some extent. Due to the complexity of pile-soil interaction and the nonlinearity of pile perimeter resistance and displacement, single-parameter methods such as the m or c method often struggle to produce results that accurately reflect both pile displacement and internal force calculations. As a result, Wu [29] proposed a comprehensive stiffness principle and two-parameter method for calculating thrust piles, which established a power series solution. However, this solution does not account for the longitudinal and transverse bending that occurs in bridge piles or high bearing piles under horizontal and vertical (or inclined) loads, known as the gravity second-order effect. Obtaining an analytical solution for power series that accounts for longitudinal and transverse bending is challenging. In summary, there are relatively few research results available on elastic pile foundations in viscoelastic media.

In other area outside of soil mechanics, it is imperative to consider the post-buckling behavior when the supporting spring system exhibits nonlinearity or softening effects. In such scenarios, subcritical and unstable post-buckling behavior with localized failure modes may occur, rendering the failure mode highly sensitive to imperfections. Consequently, linear buckling analysis may provide limited guidance in predicting the ultimate load. An example of this phenomenon is the buckling of railways subjected to thermal strain [30,31,32]. And in other area outside of civil engineering, the viscoelastic-elastic complex system is still widely received by researchers [33,34,35,36,37].

In this paper, the comprehensive stiffness principle and two-parameter method are extended to the calculation of longitudinally and transversely bent piles, analogous to longitudinally and transversely bent beams in mechanics of materials [38,39]. The deflection differential equation of a composite loaded pile is established using this method and the corresponding finite difference equation is derived. The finite difference iterative format for pile deformation and internal force is then obtained based on the pile top load, pile end restraint conditions, and the continuity of deformation and internal force at the interface between the free and embedded solid sections of the pile. This study presents a time-history analysis of the displacement behavior of elastic compression bending bar embedded in viscoelastic media. The primary objective of this research is to enhance the long-term stability of elastic bars subjected to lateral disturbances in practical engineering applications. Additionally, this study addresses a significant gap in the existing literature by providing a comprehensive theoretical framework for understanding the temporal variation of displacement in elastic bars situated within viscoelastic media. The findings of this research have important implications for the design and construction of soft soil foundation engineering projects, offering a robust theoretical foundation for future work in this area.

2. Methods

2.1. Deflection Function for an Elastic Compression-Bending Bar within an Elastic Medium

2.1.1. Basic Model of Viscoelastic Mechanics

In the southeastern coastal region of China, where the soft soil layer is relatively thick, rock-socketed piles are commonly employed in the design of bridge pile foundations to ensure the stability of the superstructure. Pile lengths often exceed 50 m and can even reach 100 m in some areas. As such, the stability of single piles under these conditions can be approached as a slender bar stability problem. Unlike a compression bar model with a fixed bottom and free top, the pile body is enveloped by a thick layer of soft soil and cannot be simplified as a compression bar stability research model without lateral support. Similar to the support conditions of an elastic foundation beam, the pile side can be transformed into an elastic foundation support based on the Winkler foundation assumption. In actual working conditions, pile foundations bear not only axial loads transmitted from the bridge superstructure but also additional horizontal stresses caused by lateral deviation of piles, resulting in a state of compression-bending.

The constitutive equation of linear viscoelastic materials can be bifurcated into two primary categories: differential and integral types. The differential type is fundamentally premised on the utilization of ideal elements to construct a variety of models. These models subsequently facilitate the derivation of material constitutive equations, which are based on the inherent properties and combinational forms of the elements. This approach offers a level of convenience in resolving certain problems. As per the foundational theory of the differential type, it can be postulated that the one-dimensional differential constitutive equation of any given linear viscoelastic material can be articulated as follows:

$$P\sigma =Q\epsilon$$

(1)

$$P={\displaystyle \sum _{k=0}^n{p}_k\frac{d^k}{d{t}^k}}, Q={\displaystyle \sum _{k=0}^n{q}_k\frac{d^k}{d{t}^k}}$$

(2)

where: the differential operators P and Q are shown in Equation (2), where p_k and q_k are constants related to the elastic and viscous properties (k is only a counting symbol, not the subgrade reaction coefficient mentioned below).

The constitutive equation of viscoelastic materials is derived from element models, which are formulated through the integration of elastic (spring) and viscous (dashpot) elements in a multitude of configurations. The spring elements adhere to the principles of Hooke’s law.

$$\sigma =E\epsilon , \tau =G\gamma$$

(3)

where σ and τ are the normal stress and shear stress, respectively; ε and γ are the normal strain and shear strain, respectively; E and G are the tensile-compressive elastic modulus and shear elastic modulus.

The elements of a viscoelastic model, known as dashpots, adhere to Newton’s law of viscosity, which is a fundamental principle in fluid mechanics as follow:

$$\sigma =\eta \stackrel{·}{\epsilon }, \tau ={\eta }_1\stackrel{·}{\gamma }$$

(4)

where η and η₁ is the coefficient of viscosity, $\stackrel{·}{\epsilon }$ and $\stackrel{·}{\gamma }$ are the normal strain and shear strain of viscosity. Figure 1 shows the mechanical model of Spring and Dashpot element.

The Kelvin model is composed of a spring element and a dashpot element in parallel, as is shown in Figure 2. And the Van Der Poel model is formed by a spring element and a Kelvin model in series as is shown in Figure 3.

The constitutive equation of Van Der Poel is

$$\sigma +p_1\stackrel{·}{\sigma }=q_0\epsilon +q_1\stackrel{·}{\epsilon }$$

(5)

where:

$$p_1=\frac{{\eta }_1}{E_1+E_2}, q_0=\frac{E_1{E}_2}{E_1+E_2}, q_1=\frac{E_2{\eta }_1}{E_1+E_2}$$

(6)

Given a constant stress σ₀, the creep equation of the Van Der Poel model is as follows:

$$\epsilon \left(t\right)=\frac{{\sigma }_0}{E_2}+\frac{{\sigma }_0}{E_1}\left(1-e^{-\frac{E_1}{\eta }t}\right)$$

(7)

2.1.2. Basic Physical Consumption of Elastic-Viscoelastic Media

The Winkler Foundation Model, also referred to as the Local Elastic Foundation Model, was a seminal concept introduced by the Czechoslovakian engineer Winkler in 1867. This model primarily postulates that the settlement of any given point on the foundation surface is directly proportional to the pressure exerted on the unit area at that point. In accordance with this model, the more rigid rock layer at the base of the foundation is considered a rigid layer, while the upper soil layer, which exhibits lesser stiffness, is discretized into independent springs on the rigid layer, as depicted in Figure 4.

Viscoelastic mechanics encompasses two fundamental aspects: firstly, it provides a description of the viscoelastic properties of materials and expresses their constitutive relations; secondly, it involves the formulation and resolution of viscoelastic boundary value problems. As per the established principles of viscoelastic mechanics, the stress-strain relationship of elastic materials can be construed as a specific instance of viscoelasticity [40,41]. Consequently, a viscoelastic foundation can adhere to Winkler’s assumption that views the foundation soil as discrete support points. By selecting an appropriate viscoelastic element model to replace the soil spring, we can apply the solution of an elastic beam on an elastic foundation to the principle of elastic-viscoelastic correspondence. Following a process of Laplace transform and inverse transformation, we can derive the solution for an elastic beam on a viscoelastic foundation.

2.1.3. Mechanical Model of an Elastic Compression-Bending Bar in an Elastic Medium

Figure 5 depicts the mechanical model of an elastic compression-bending bar in an elastic medium. A constant load with a width of L_p and magnitude of p acts on the surface of the elastic medium at a distance of d_p from one side of the bar to simulate subgrade filling on one side of the pile foundation; q(x) represents the real passive load on the bar body caused by pile loading and soil self-weight; axial force P acts on the end of the bar, inducing a small displacement Δx in the vertical direction at the end of the bar. The bottom of the bar is fixed and the top is free; the length of bar is L, and k(x) represents the ground reaction coefficient at different depths; b₁ denotes the calculated width of the elastic bar.

Drawing upon the relevant theories of elasticity and soil mechanics, it can be shown that, within an elastic medium, the horizontal additional stress σ_x at a given point in the medium induced by a constant load p of a certain width can be determined based on the Boussinesq solution and Flamant solution [42,43]. Subsequently, the real passive load q(x) on the bar body can be obtained demonstrated in Equation (8) which is detailed introduced in Appendix A.

$$q\left(x\right)=\frac{b_{1}p}{\pi }\left(\frac{{\left(L-x\right)}^2\arctan \frac{L_p}{L-x}}{{\left(L-x\right)}^2+\left(L_p+d_p\right)d_p}-\frac{\left(L_p+d_p\right)\left(L-x\right)}{{\left(L-x\right)}^2+{\left(L_p+d_p\right)}^2}+\frac{\left(L-x\right)d_p}{{\left(L-x\right)}^2+d_p{}^2}\right)$$

(8)

where: q(x) is passive load on elastic bar caused by constant load and soil self-weight

b₁ is calculated width of the elastic bar

L is length of bar

L_p is width of constant load

p is magnitude of constant load

d_p is distance from the bar to constant load

The objective of this study is to investigate the displacement response of piles under the influence of ground load in the adjacent area. A crucial step in this analysis is to ascertain the additional pile load induced by the ground load. To this end, ABAQUS is employed as finite element model (FEM) solution software, pile material and geometric parameters are applied with roadbed filling area layout parameters, etc. that are specified at the outset of this section. Detailed parameters are shown in

Table 1. Parameters in ABAQUS model.

Parameter	Magnitude	Unit
Width of bar/b1	1	m
Length of bar/L	50	m
Elastic Modulus of bar/E	3 × 104	MPa
Density of bar/ρ	2385	kg/m3
Poisson’s ratio/λ	0.2	/
Distance from the bar to constant load/dp	1	m
Magnitude of constant force/q	50	kPa
Length of constant force/Lp	20	m
Reaction force coefficient/k	4500	kN/m4

and entire model are shown as Figure 6, which has a grey solid part for soil and a constant force adding in the yellow part.

As is shown in Figure 7, the numerical solution is in good agreement with the analytical solution. Therefore, the passive load of the pile body can be calculated by Equation (8).

2.1.4. A Comprehensive Solution Formula for the Deflection of an Elastic Slender Bar within an Elastic Medium

Based on the law of conservation of energy, the total potential energy of both the bar and the surrounding soil Π_p can be calculated when subjected to an axial force P and a horizontal additional load q(x) as shown in Equation (9) which is detailed introduced in Appendix B:

$${\Pi }_p=\frac{EI}{2}{\displaystyle {\int }_0^L{\left(w^{″}\right)}^{2}dx}+\frac{b_1}{2}{\displaystyle {\int }_0^{L}k{w}^{2}dx}-P\Delta x-b_1{\displaystyle {\int }_0^{L}qwdx}$$

(9)

where: Π_p is total potential energy of both the bar and the surrounding soil

E is elastic modulus of elastic bar

I is moment of inertia of the cross section of elastic bar

k is reaction force coefficient of elastic medium

w is deflection of elastic bar as w(x)

Δx is a small displacement caused by constant load p

P is an axial force as shown in Figure 5

As per the mechanical model depicted in Figure 5, the boundary conditions for the compressed bar are defined as follows: x = 0, w = 0, w′ = 0; x = L, w″ = 0. Given these boundary conditions, a test function w(x) for the deflection curve of the micro-bent state bar can be selected in the form of a triangular series by the Laplace transform as:

$$w(x)={\displaystyle \sum _{n=1}^i{w}_n\left(1-\cos \frac{\left(2n-1\right)\pi x}{2L}\right)}$$

(10)

Substituting Equation (10) into Equation (9) yields:

$$\begin{array}{ll}{\Pi }_p& =\frac{EI}{2}{\displaystyle {\int }_0^L{\left({{\displaystyle \sum _{n=1}^i{w}_n\left(\frac{\left(2n-1\right)\pi }{2L}\right)}}^2\cos \frac{\left(2n-1\right)\pi x}{2L}\right)}^{2}dx}+\frac{b_1}{2}{\displaystyle {\int }_0^{L}k{\left({\displaystyle \sum _{n=1}^i{w}_n}\left(1-\cos \frac{\left(2n-1\right)\pi x}{2L}\right)\right)}^{2}dx}\\ & -\frac{P}{2}{\displaystyle {\int }_0^L{\left({\displaystyle \sum _{n=1}^i{w}_n}\frac{\left(2n-1\right)\pi }{2L}\sin \frac{\left(2n-1\right)\pi x}{2L}\right)}^{2}dx}-b_1{\displaystyle {\int }_0^{L}q{\displaystyle \sum _{n=1}^i{w}_n}\left(1-\cos \frac{\left(2n-1\right)\pi x}{2L}\right)dx}\end{array}$$

(11)

In accordance with the principle of minimum potential energy, for a system in a stable equilibrium state, Π_p must be minimized, which means the first-order derivative of Π_p is 0. This can be expressed as follows:

$$\begin{array}{ll}\frac{\partial {\Pi }_p}{\partial w_r}=& {\displaystyle \sum _{n=1}^i{w}_n{\displaystyle {\int }_0^L\left(\mathit{EI}\frac{{\left(2n-1\right)}^2{\left(2r-1\right)}^2{\pi }^4}{16L^4}\cos \frac{\left(2n-1\right)\pi x}{2L}\cos \frac{\left(2r-1\right)\pi x}{2L}\right.}}\\ & +b_{1}k\left(1-\cos \frac{\left(2n-1\right)\pi x}{2L}\right)\left(1-\cos \frac{\left(2r-1\right)\pi x}{2L}\right)\\ & \left.-P\frac{\left(2n-1\right)\left(2r-1\right){\pi }^2}{4L^2}\sin \frac{\left(2n-1\right)\pi x}{2L}\sin \frac{\left(2r-1\right)\pi x}{2L}\right)dx\\ & -b_1{\displaystyle {\int }_0^{L}q\left(1-\cos \frac{\left(2r-1\right)\pi x}{2L}\right)dx}=0,r=1,2,3,\dots ,i\end{array}$$

(12)

Take:

$$\begin{array}{ll}{W}_{nr}=& {\displaystyle {\int }_0^L\left(\mathit{EI}\frac{{\left(2n-1\right)}^2{\left(2r-1\right)}^2{\pi }^4}{16L^4}\cos \frac{\left(2n-1\right)\pi x}{2L}\cos \frac{\left(2r-1\right)\pi x}{2L}\right.}\\ & +b_{1}k\left(1-\cos \frac{\left(2n-1\right)\pi x}{2L}\right)\left(1-\cos \frac{\left(2r-1\right)\pi x}{2L}\right)\\ & \left.-P\frac{\left(2n-1\right)\left(2r-1\right){\pi }^2}{4L^2}\sin \frac{\left(2n-1\right)\pi x}{2L}\sin \frac{\left(2r-1\right)\pi x}{2L}\right)dx\end{array}$$

(13)

$$Q_r=b_1{\displaystyle {\int }_0^{L}q\left(1-\cos \frac{\left(2r-1\right)\pi x}{2L}\right)dx}$$

(14)

Then Equation (11) can be transformed into Equation (14):

$${\displaystyle \sum _{r=1}^i{\displaystyle \sum _{n=1}^i\left(w_n{W}_{nr}-Q_r\right)}}=\left[\begin{array}{cccc}{W}_{11}& W_{12}& \cdots & W_{1i}\\ W_{21}& W_{12}& \cdots & W_{i2}\\ ⋮& ⋮& \cdots & ⋮\\ W_{i1}& W_{i2}& \cdots & W_{ii}\end{array}\right]\left[\begin{array}{c}{w}_1\\ w_2\\ ⋮\\ w_i\end{array}\right]-\left[\begin{array}{c}{Q}_1\\ Q_2\\ ⋮\\ Q_i\end{array}\right]=0$$

(15)

Equation (15) delineates the phase space formulation representation of an elastic bar within an elastic medium by the Laplace transform, which is prepared for viscoelastic medium illustrated in Equation (23). The time-dependency and spatial-dependency of the viscoelastic medium are transposed into a phase space via the Laplace transform. This process effectively decouples the time-dependency from the deformation process inherent in the viscoelastic medium. Consequently, the deformation issue concerning the elastic bar within the viscoelastic medium, as referenced in this study, is reduced to a problem devoid of time-dependency. Upon completion of the solution, the Laplace inverse transformation is employed. This reconfigures the deformation function of the bar from the resolved phase space formulation back to a general solution exhibiting both time-dependency and spatial-dependency. The outcome is a universal function capable of predicting terminal deformations of the bar.

By solving Equation (15), the coefficients w_n of the deflection curve test function for the bar body can be obtained. This yields an approximate solution w(x) for the true displacement curve of the bar under varying axial loads P. For instance, when i = 2 and k = m(L − x), solving Equation (15) results in:

$$w_1=\frac{288L^3{\pi }^2\left(81{\pi }^4\left(9EI{\pi }^2-4L^{2}P\right)Q_1+8L^{5}m{b}_1\left(\begin{array}{l}\left(27{\pi }^2-28\right)Q_1\\ +2\left(62-9{\pi }^2\right)Q_2\end{array}\right)\right)}{\begin{array}{l}729{\pi }^8\left(EI{\pi }^2-4L^{2}P\right)\left(9EI{\pi }^2-4L^{2}P\right)\\ +16L^{5}m{b}_1\left(\begin{array}{l}41328L^{2}P{\pi }^4-180\left(511EI+27L^{2}P\right){\pi }^6\\ +9963EI{\pi }^8+4L^{5}m{b}_1\left(405{\pi }^4-3096{\pi }^2-8320\right)\end{array}\right)\end{array}}$$

(16)

$$w_2=\frac{288L^3{\pi }^2\left(9{\pi }^4\left(EI{\pi }^2-4L^{2}P\right)Q_2+8L^{5}m{b}_1\left(\begin{array}{l}2\left(62-9{\pi }^2\right)Q_1\\ +9\left(3{\pi }^2-28\right)Q_2\end{array}\right)\right)}{\begin{array}{l}729{\pi }^8\left(EI{\pi }^2-4L^{2}P\right)\left(9EI{\pi }^2-4L^{2}P\right)\\ +16L^{5}m{b}_1\left(\begin{array}{l}41328L^{2}P{\pi }^4-180\left(511EI+27L^{2}P\right){\pi }^6\\ +9963EI{\pi }^8+4L^{5}m{b}_1\left(405{\pi }^4-3096{\pi }^2-8320\right)\end{array}\right)\end{array}}$$

(17)

Thus, it can be inferred that within the context of a second-order triangular series, the deflection function for the bar body is:

$$w(x)=w_1\left(1-\cos \frac{\pi x}{2L}\right)+w_2\left(1-\cos \frac{3\pi x}{2L}\right)$$

(18)

2.2. Deflection Function for an Elastic Compression-Bending Bar in a Viscoelastic Medium

2.2.1. Mechanical Model of an Elastic Compression-Bending Bar in a Viscoelastic Medium

The interaction between elastic compression bending bar and the surrounding soil is complex and can be either static or dynamic. To simplify calculations, this interaction is often modeled using an elastic foundation beam. However, this approach has limitations as soil is a discrete, multiphase system with nonlinear and irrecoverable stress-strain characteristics under external forces. As such, modeling soil solely with springs may not accurately capture its behavior. Viscoelastic mechanics can better characterize the creep properties of soft soils by revealing their rheological properties. This typically involves establishing the viscoelastic instantaneous relationship of the material and determining boundary conditions and their solutions. By transforming the parameters characterizing the foundation soil’s material properties in the elastic foundation model into time-dependent viscoelastic parameters, the springs in the elastic foundation can be replaced with a suitable viscoelastic model while maintaining the boundary conditions of the elastic foundation model. This allows for the solution of the beam on the elastic foundation to be transformed into a time-dependent viscoelastic solution.

In this part, adopts the method of the Winkler foundation hypothesis, selects a suitable viscoelastic constitutive model, and simplifies the soil around the slender bar shown in Figure 5 into discrete viscoelastic element supports. Taking the Van Der Poel model as an example, as shown in Figure 8.

2.2.2. A Comprehensive Solution Formula for the Deflection of an Elastic Slender Bar within a Viscoelastic Medium

In accordance with the principle of elastic-viscoelastic correspondence, the viscoelastic foundation elastic beam problem may be considered as an elastic foundation elastic beam problem with a reaction force coefficient of $\overline{k}\left(s,x\right)$. By analogy, utilizing the m-method to determine the value of the foundation reaction force coefficient, the foundation reaction force coefficient within the phase space formulation is established to be:

$$\overline{k}(s,x)=A{\overline{k}}_0(s)\left(L-x\right)$$

(19)

where: $\overline{k}$ is reaction force coefficient of viscoelastic medium in kelvin model

A is constant

s is displacement of elastic bar

Within Equation (19), A represents a constant that governs the rate of change for the foundation reaction force coefficient $\overline{k}\left(s,x\right)$ from the ground surface to the depth of the soil. By substituting Equation (19) into Equation (13), the result is obtained:

$$\begin{array}{ll}{\overline{W}}_{nr}=& {\displaystyle {\int }_0^L\left(\mathit{EI}\frac{{\left(2n-1\right)}^2{\left(2r-1\right)}^2{\pi }^4}{16L^4}\cos \frac{\left(2n-1\right)\pi x}{2L}\cos \frac{\left(2r-1\right)\pi x}{2L}\right.}\\ & +A{\overline{k}}_0\left(s\right)\left(L-x\right)b_1\left(1-\cos \frac{\left(2n-1\right)\pi x}{2L}\right)\left(1-\cos \frac{\left(2r-1\right)\pi x}{2L}\right)\\ & \left.-P\frac{\left(2n-1\right)\left(2r-1\right){\pi }^2}{4L^2}\sin \frac{\left(2n-1\right)\pi x}{2L}\sin \frac{\left(2r-1\right)\pi x}{2L}\right)dx,n=1,2,\dots ;r=1,2,\dots \end{array}$$

(20)

In the event of an instantaneous increase in the side load of the bar, namely q(t,x) = q(x)H(t), upon substitution into Equation (14), the result is obtained:

$$Q_r\left(t\right)=H\left(t\right)b_1{\displaystyle {\int }_0^{L}q\left(1-\cos \frac{\left(2r-1\right)\pi x}{2L}\right)dx=H\left(t\right)}{Q}_r$$

(21)

where: W, Q and H are generalized function

Upon conducting Laplace transform on Equation (21), the result is obtained:

$${\overline{Q}}_r=Q_r/s$$

(22)

Analogous to Equation (15), the subsequent matrix expression may be derived:

$${\displaystyle \sum _{r=1}^i{\displaystyle \sum _{n=1}^i\left({\overline{w}}_n\left(s\right){\overline{W}}_{nr}-{\overline{Q}}_r\right)}}=\left[\begin{array}{cccc}{\overline{W}}_{11}& {\overline{W}}_{12}& \cdots & {\overline{W}}_{1i}\\ {\overline{W}}_{21}& {\overline{W}}_{12}& \cdots & {\overline{W}}_{i2}\\ ⋮& ⋮& \cdots & ⋮\\ {\overline{W}}_{i1}& {\overline{W}}_{i2}& \cdots & {\overline{W}}_{ii}\end{array}\right]\left[\begin{array}{c}{\overline{w}}_1\left(s\right)\\ {\overline{w}}_2\left(s\right)\\ ⋮\\ {\overline{w}}_i\left(s\right)\end{array}\right]-\left[\begin{array}{c}{\overline{Q}}_1\\ {\overline{Q}}_2\\ ⋮\\ {\overline{Q}}_i\end{array}\right]=0$$

(23)

where: $\overline{W}$ and $\overline{Q}$ are generalized function under reaction force coefficient of viscoelastic medium

By solving Equation (23), the coefficients ${\overline{w}}_n\left(s\right)$ can be determined. Subsequently, the solution w(s,x) for the deflection function of the elastic bar within the phase space formulation can be obtained. By performing an inverse Laplace transform on this solution, the deflection function w(t,x) of the elastic bar within a viscoelastic medium can be derived.

2.2.3. Deflection Solution of an Elastic Slender Bar in a Kelvin Medium

The constitutive equation for a Kelvin body is given by $\sigma =q_0\epsilon +q_1\dot{\epsilon }$, where q₀ represents the Young’s modulus (E) and q₁ represents the viscosity coefficient (η). By applying a Laplace transform to this constitutive equation and subsequently rearranging it, the reference value for the ground reaction coefficient ${\overline{k}}_0\left(s\right)$ can be determined:

$${\overline{k}}_0\left(s\right)=\frac{\overline{\sigma }}{\overline{\epsilon }}=E+\eta s$$

(24)

where: ${\overline{k}}_0$ is reference value for the ground reaction coefficient

$\overline{\sigma }$ is viscoelastic stress

$\overline{\epsilon }$ is viscoelastic strain

Consequently, the ground reaction coefficient $\overline{k}\left(s, x\right)$ within the phase space formulation, as derived from the Kelvin model, is given by:

$$\overline{k}(s,x)=A\left(E+\eta s\right)\left(L-x\right)$$

(25)

Continuing with the example where i = 2, Equation (23) can be solved to determine the expression for the coefficient ${\overline{w}}_n\left(s\right)$, where E_k denotes the elastic modulus of the spring element within the Kelvin model.

$${\overline{w}}_1\left(s\right)=\frac{288L^3{\pi }^2\left(81{\pi }^4\left(9EI{\pi }^2-4L^{2}P\right)Q_1+8A{b}_1{L}^5\left(\begin{array}{l}\left(27{\pi }^2-28\right)Q_1\\ +2\left(62-9{\pi }^2\right)Q_2\end{array}\right)\left(E_k+\eta s\right)\right)}{s\left(\begin{array}{l}729{\pi }^8\left(EI{\pi }^2-4L^{2}P\right)\left(9EI{\pi }^2-4L^{2}P\right)\\ +144A{b}_1{L}^5{\pi }^4\left(4592L^{2}P-20\left(511EI+27L^{2}P\right){\pi }^2+1107EI{\pi }^4\right)\left(E_k+\eta s\right)\\ +64A^2{b}_1{}^2{L}^{10}\left(405{\pi }^4-3096{\pi }^2-8320\right){\left(E_k+\eta s\right)}^2\end{array}\right)}$$

(26)

$${\overline{w}}_2\left(s\right)=\frac{288L^3{\pi }^2\left(36L^{2}P{\pi }^4{Q}_2-9EI{\pi }^6{Q}_2+8A{b}_1{L}^5\left(\begin{array}{l}2\left(9{\pi }^2-62\right)Q_1\\ +9\left(28-3{\pi }^2\right)Q_2\end{array}\right)\left(E_k+\eta s\right)\right)}{s\left(\begin{array}{l}256A^2{b}_1{}^2{L}^{10}{\left(62-9{\pi }^2\right)}^2{\left(E_k+\eta s\right)}^2\\ -9\left(\begin{array}{l}EI{\pi }^6-4L^{2}P{\pi }^4+\\ 8A{b}_1{L}^5\left(3{\pi }^2-28\right)\left(E_k+\eta s\right)\end{array}\right)\left(\begin{array}{l}81\left(9EI{\pi }^6-4L^{2}P{\pi }^4\right)\\ +8A{b}_1{L}^5\left(27{\pi }^2-28\right)\left(E_k+\eta s\right)\end{array}\right)\end{array}\right)}$$

(27)

By applying an inverse Laplace transform to Equations (26) and (27), the two coefficients w₁(t) and w₂(t) for the deflection function w(t,x) can be determined as:

$$w_1\left(t\right)=\frac{144{\mathrm{e}}^{-N_{7}t}{L}^3{\pi }^2}{\left(N_3+E_k\left(N_4+E_k{N}_5\right)\right)N_6}\left(\begin{array}{l}\left({\mathrm{e}}^{N_{8}t}-1\right)\left(\begin{array}{l}{N}_2\left(2N_3+E_k{N}_4\right)\\ -N_1\left(N_4+2E_k{N}_5\right)\end{array}\right)\eta \\ +\left(N_1{N}_6+E_k{N}_2{N}_6\right)\left(2{\mathrm{e}}^{N_{7}t}-{\mathrm{e}}^{N_{8}t}-1\right)\end{array}\right)$$

(28)

$$w_2\left(t\right)=R_1\left(R_{11}+\frac{{\mathrm{e}}^{R_{9}t}\left(3R_{12}\left(R_2+E_k{R}_3\right)+R_{13}\eta \right)+{\mathrm{e}}^{R_{10}t}\left(3R_{12}\left(R_2+E_k{R}_3\right)-R_{13}\eta \right)}{6R_{12}\left(9R_5{R}_7+E_k\left(9R_6{R}_7-E_k{R}_4+9R_8\left(R_5+E_k{R}_6\right)\right)\right)}\right)$$

(29)

where:

$$\begin{array}{l}{N}_1=81{\pi }^4\left(9EI{\pi }^2-4L^{2}P\right)Q_1; N_2=8A{b}_1{L}^5\left(\left(27{\pi }^2-28\right)Q_1+2\left(62-9{\pi }^2\right)Q_2\right)\\ N_3=729{\pi }^8\left(EI{\pi }^2-4L^{2}P\right)\left(9EI{\pi }^2-4L^{2}P\right);\\ N_4=144A{b}_1{L}^5{\pi }^4\left(4592L^{2}P-20\left(511EI+27L^{2}P\right){\pi }^2+1107EI{\pi }^4\right);\\ N_5=64A^2{b}_1{}^2{L}^{10}\left(405{\pi }^4-3096{\pi }^2-8320\right);\\ N_6=\sqrt{\left(N_4{}^2-4N_3{N}_5\right){\eta }^2}; N_7=\frac{N_4\eta +2E_k{N}_5\eta +N_6}{2N_5{\eta }^2}; N_8=\frac{N_6}{N_5{\eta }^2};\\ R_1=288L^3{\pi }^2; R_2=9{\pi }^4\left(4L^{2}P-EI{\pi }^2\right)Q_2;\\ R_3=8A{b}_1{L}^5\left(2\left(9{\pi }^2-62\right)Q_1+9\left(28-3{\pi }^2\right)Q_2\right);\\ R_4=256A^2{b}_1{}^2{L}^{10}{\left(62-9{\pi }^2\right)}^2; R_5=EI{\pi }^6-4L^{2}P{\pi }^4; R_6=8A{b}_1{L}^5\left(3{\pi }^2-28\right)\\ R_7=81\left(9EI{\pi }^6-4L^{2}P{\pi }^4\right);\\ R_8=8A{b}_1{L}^5\left(27{\pi }^2-28\right); R_9=\frac{18E_k{R}_6{R}_8+9R_5{R}_8+9R_6{R}_7-2E_k{R}_4}{2\left(R_4-9R_6{R}_8\right)\eta }-\frac{3R_{12}}{2\left(R_4-9R_6{R}_8\right){\eta }^2};\\ R_{10}=\frac{18E_k{R}_6{R}_8+9R_5{R}_8+9R_6{R}_7-2E_k{R}_4}{2\left(R_4-9R_6{R}_8\right)\eta }+\frac{3R_{12}}{2\left(R_4-9R_6{R}_8\right){\eta }^2};\\ R_{11}=\frac{R_2+E_k{R}_3}{E_k{}^2{R}_4-9R_5{R}_7-9E_k{R}_6{R}_7-9E_k{R}_5{R}_8-9E_k{}^2{R}_6{R}_8};\\ R_{12}=\sqrt{\left(4R_4{R}_5{R}_7+9R_6{}^2{R}_7{}^2-18R_5{R}_6{R}_7{R}_8+9R_5{}^2{R}_8{}^2\right){\eta }^2}\\ R_{13}=9\left(R_2{R}_6{R}_7+R_2{R}_5{R}_8-2R_3{R}_5{R}_7\right)-9E_k{R}_3\left(R_6{R}_7+R_5{R}_8\right)-2E_k{R}_2\left(R_4-9R_6{R}_8\right)\end{array}$$

As a result, the deflection function w(t,x) for an elastic bar surrounded by a Kelvin body within a second-order triangular series has been determined. If the surrounding medium of the bar is either a Van Der Poel or Burgers body, the deflection function w(t,x) for the elastic bar can also be calculated at varying levels of solution accuracy using the aforementioned method.

3. Results and Discussion

According to the rheological test data obtained from soft soil samples in the engineering site selection area, the Burgers model was employed as the viscoelastic constitutive model for the bar-side medium in theoretical calculations. The constitutive model parameters were determined to be E_M = 116.2270 MPa, E_K = 7.0203 MPa, η_M = 511.5671 MPa·d, and η_K = 8.6031 MPa·d. The soil sample utilized in this rheological test was extracted from a depth of 9–11 m below ground level, with a median value of 10 m being selected. Consequently, in the phase space formulation, ${\overline{k}}_0\left(s\right)$ represents the elastic ground reaction coefficient at a depth of 10 m below ground level. As per the definition of the ground reaction coefficient by the matlock method, the ground reaction coefficient increases linearly from 0 at ground level to the underground depth. Based on this definition, it was determined that the value of A in Equation (25) is 0.1. All solution is solved by MATLAB based on extension of the Burgers model’s solution in the Section 2 above, by transforming K in Equation (25) into a basic burgers ground reaction coefficient which is treated under Laplace transform as $\overline{k}\left(s,x\right)=A\frac{E_1{\eta }_2\left(E_3+{\eta }_{3}s\right)s}{{\eta }_2\left(E_3+{\eta }_{3}s\right)s+E_1\left(E_3+\left({\eta }_2+{\eta }_3\right)s\right)}\left(L-x\right)$. Detail is shown in Appendix C.

To enhance the accuracy of the calculation results, a solution accuracy of a seventh-order triangular series was employed. At this level of accuracy, the form of the coefficient ${\overline{w}}_n\left(s\right)$ solved according to Equation (23) is highly complex and obtaining an analytical solution for w_n(t) through inverse Laplace transform is not straightforward. Abate J and Whitt W summarized existing research on Laplace numerical inverse transform in literature [44] and proposed a unified theoretical framework for Laplace inverse transformation numerical inversion. In this study, a MATLAB program plug-in developed by Tucker McClure based on this theoretical framework was utilized to write a calculation program and complete subsequent case calculations.

To investigate the compressive stability of longitudinal and transverse bending elastic bars in viscoelastic media, the value of bar end displacement w(L) under long-term axial force and lateral loading was taken as a key parameter. This study examines the compressive stability of elastic compression bending bars in viscoelastic media.

3.1. Time Dependence of the Effect of Axial Force on the Lateral Deflection of the Bar Top

When an elastic bar is situated within a viscoelastic medium, its deflection over time can be determined using the deflection function Equation (10) of the longitudinal and transverse bending elastic bar. In this section, we investigate the compressive stability of the longitudinal and transverse bending elastic bar within the viscoelastic medium using the value of the bar end displacement w(L) under long-term lateral pile load as a key parameter.

Taking a lateral pile load strength of 50 kPa and a bar length of 50 m as examples, and with a background of a typical bridge design life of 100 years in engineering, the total calculation time was set to 100 years. The change in transverse deflection, w(L), of the bar top was investigated when different magnitudes of axial force, P, were applied to the bar top, as shown in Figure 9a. Analysis of Figure 9a reveals that within the 100-year calculation period, when P is small (e.g., 50 kN, 100 kN, 500 kN, 1000 kN, 2500 kN), w(L) maintains an approximate linear growth trend within the long-term load action time after the initial stage.

When P is 5000 kN or 7500 kN, w(L) exhibits a linear growth stage in the early period within the calculation period of 100 years, followed by a trend of nonlinear acceleration growth with time in the later period. When P is 10,000 kN, w(L) does not reach 100 years and tends towards +∞ at a certain time point during this period. This indicates that when different axial forces P are applied to the bar top and if the load action time is sufficiently large, they will all cause the bar top displacement w(L) to tend towards +∞.

When P is small, only its linear growth stage can be observed within the calculation period; when P increases to a certain level, w(L) will exhibit nonlinear acceleration growth within the calculation period but w(L) → +∞ has not yet occurred; when P is large, w(L) → +∞ can be observed within the calculation period. By analogy with the definition of critical force P_cr that causes simple compression bars to become unstable, we define the time that causes bars to become unstable as critical instability time t_cr.

The above results show that if w(L) → +∞ is used to define bar instability, different axial forces P correspond to corresponding critical instability times t_cr; alternatively, specifications can be referred to and w(L) = 10 cm can be used as the basis for determining critical instability time t_cr.

In order to clarify the characteristics of bar displacement throughout the process and in conjunction with Burgers body’s creep law, it is necessary to study the deformation morphology of long-term loaded bars. Figure 9b shows the w(L)-t curve of bars subjected to axial force and horizontal additional load in their initial stage. Analysis of this figure reveals that instantaneous deformation occurs in bars at the moment of load action and that larger axial forces result in larger instantaneous deformation amounts.

As time progresses, w(L) experiences rate of increase with time first increases and then decreases before finally entering a stable uniform growth stage. Larger axial forces result in larger growth amplitudes at each time point during nonlinear growth stage of w(L) and also result in larger rates during its uniform growth stage. On the graph it shows that when axial force is large, w(L)-t curve envelopes w(L)-t curve when axial force is small from top to bottom.

In this section, considered a scenario as where the bar-side pile loading intensity is 50 kPa and the bar length is 50 m. The total calculation time is set to 100 years. We investigate the temporal evolution of lateral deflection w(L) at the top of the bar under varying magnitudes of axial force P as depicted in Figure 9.

Figure 9a,b illustrate that the bar undergoes instantaneous deformation upon the application of load. As time progresses, the rate of increase of w(L) with respect to time initially increases before decreasing, ultimately reaching a stable and uniform growth phase. This displacement process closely resembles the three-stage creep characteristics of the Burgers model.

Over the course of the 100-year calculation period, when P is small, w(L) exhibits an approximately linear growth trend during the long-term load application period following the initial stage. When P is larger, such as 5000 kN or 7500 kN, w(L) undergoes a linear growth phase in the early stage and exhibits a nonlinear acceleration trend with time in the later stage. When P continues to increase to 10,000 kN, w(L) does not reach 100 years and tends towards +∞ at some point during that period. This indicates that when different axial forces P are applied to the bar top and the load application time is sufficiently large, it will cause the bar top displacement w(L) to tend towards +∞.

When P is small, only its linear growth stage can be observed within the calculation period. When P increases to a certain level, within the calculation period, w(L) will exhibit nonlinear acceleration growth but has not yet shown the phenomenon of w(L) → +∞. When P is large, w(L) → +∞ can be observed within the calculation period. Analogous to the definition of the critical force P_cr that causes instability in a compressed bar, the time that causes bar instability is defined as the critical instability time t_cr.

These results demonstrate that if w(L) → +∞ is used to define bar instability, different axial forces P correspond to their respective critical instability times t_cr. However, given that w(L) → +∞ is inconsistent with engineering practice and that w(L) → +∞ results cannot be obtained under common values of P for vertical ultimate bearing capacity of pile foundations within a 100-year calculation period, it is recommended that w(L) = 10 cm, as stipulated by specification [45], be utilized as the basis for determining the critical instability time t_cr.

Figure 10 illustrates that if w(L) = 10 cm, as stipulated by the specification, is employed as the control condition for defining t_cr, the magnitude of t_cr decreases as the axial force increases. This result suggests that when the bar length and bar-side pile loading intensity are held constant, an increase in axial force has a detrimental effect on the stability of the bar and reduces the time required for the bar displacement to exceed its limit.

3.2. Time Dependence of the Effect of Bar Length on the Lateral Deflection of the Bar Top

By maintaining the lateral pile load strength of the bar at a constant value of 50 kPa and using axial force values of P = 2500 kN, 5000 kN, 7500 kN, and 10,000 kN, the transverse deflection w(L) of the bar top over time was calculated for bar lengths ranging from 50 m to 140 m. Figure 11 illustrated the results obtained over a calculation period of 100 years. Analysis of the data reveals that when the axial force assumes smaller values (e.g., 2500 kN and 5000 kN), the correlation between bar length and the time course of bar displacement is minimal within the calculation period. Conversely, when the axial force assumes larger values (e.g., 7500 kN and 10,000 kN), an increase in bar length results in an initial slight decrease in the time course of bar instability, followed by an increase within the calculation period.

By maintaining the loading intensity of the bar-side pile at 50 kPa and considering axial forces, denoted as P, of 2500 kN, 5000 kN, 7500 kN, and 10,000 kN for computing the temporal evolution of the lateral deflection, denoted as w(L), at the top of the bar for bar lengths ranging from 50 m to 140 m over a period of 100 years. Based on the obtained results, when w(L) reaches 10 cm, the relationship between bar length and the critical instability time, denoted as t_c as shown in Figure 12.

As illustrated in Figure 12, the t_cr-L curve remains relatively horizontal within the small bar length interval when the axial force is 2500 kN, 5000 kN, and 7500 kN. With an increase in L, t_cr experiences a slight decrease before beginning to increase. For larger axial forces such as 10,000 kN and 12,500 kN, t_cr remains relatively constant before continuously increasing with L. Given that pile foundation lengths between 50 m and 80 m are prevalent in deep soft foundation areas, it can be inferred that within this range of bar lengths, the instability time of bars is minimally impacted by bar length under equal pile loading intensity. However, for bar lengths exceeding 80 m, the relationship between bar length L and bar instability time t_cr under varying axial load levels should be carefully considered.

Figure 13 depicts the relationship between axial force P and critical instability time t_cr for different values of bar length L. It is evident that as axial force P increases, the critical instability time t_cr decreases more rapidly for smaller bar lengths L. Additionally, for bar lengths L less than 80 m, there is a high degree of overlap in the t_cr-L curve. This observation aligns with the conclusion that for bars with a length less than 80 m, bar length has minimal impact on its instability time.

3.3. Time Dependence of the Effect of Bar-Side Pile Loading Intensity on the Lateral Deflection of the Bar Top

For a bar length L of 50 m, the change in lateral deflection w(L) of the bar top over time was calculated for axial forces of 2500 kN, 5000 kN, 7500 kN, and 10,000 kN respectively and pile loading intensities on the side of the bar of 10 kPa, 20 kPa, 30 kPa, 40 kPa, and 50 kPa. The calculation period spanned 100 years. Based on these results, the relationship between pile loading intensity and critical instability time t_cr when w(L) = 10 cm was determined as in Figure 13.

In order to accurately reflect the real-world deviation of pile foundations, the critical time t_cr for bar instability is defined using the control condition w(L) = 10. For instance, when P = 10,000 kN and the pile load intensity is 10 kPa, 20 kPa, 30 kPa, 40 kPa, and 50 kPa, the corresponding values of t_cr-P when w(L) reaches 10 cm are 11.4411 years, 6.81096 years, 4.81096 years, 3.67397 years, and 2.96164 years respectively. These results demonstrate that lower pile load intensities result in longer times for bars to exceed code requirements for deviation under the combined influence of axial force and horizontal additional load. As pile load intensity increases, this time decreases at a diminishing rate. Furthermore, t_cr-P curves corresponding to lower pile load intensities exhibit larger average slopes within the same interval compared to those corresponding to higher pile load intensities. This indicates that as axial force increases, the time required for bars to exceed code requirements for deviation decreases more rapidly at lower pile load intensities.

As illustrated in Figure 14, the t_cr-P curve is inversely related to pile loading intensity, with the magnitude of the downward shift decreasing continuously. For instance, when P = 10,000 kN and pile loading intensity is 10 kPa, 20 kPa, 30 kPa, 40 kPa, and 50 kPa respectively, the time t_cr at which w(L) reaches 10 cm is 11.4411 years, 6.81096 years, 4.81096 years, 3.67397 years, and 2.96164 years. This observation suggests that for smaller pile loading intensities, the time required for a bar to produce an offset exceeding specification requirement under the combined action of axial force and additional horizontal load is longer. As pile loading intensity increases, this time decreases continuously with a decreasing rate of decrease. Furthermore, for t_cr-P curves with smaller pile loading intensities, their average slope within the same interval is larger than that for larger pile loading intensities. This indicates that for smaller pile loading intensities, an increase in axial force results in a faster decrease in time for bars to produce an offset exceeding specification requirement.

4. Conclusions

This study constructs a mechanical model of an elastic bar constrained by elastic and viscoelastic medium, fixed at the base and free at the top under Winkler foundation hypothesis, The Rayleigh-Ritz method is employed to derive the deflection function of the bar when subjected to both axial force and local horizontal distributed load. The primary conclusions are as follows:

• A mechanical model of an elastic slender compression-bending bar, constrained by elastic-viscoelastic medium wrapping and fixed at the base while free at the top, was established. Utilizing the energy method criterion and the Rayleigh-Ritz method, an approximate solution for the deflection function of the bar body was derived when subjected to axial force P at the top and horizontal additional load caused by pile loading at the side.
• Employing the elastic-viscoelastic correspondence principle, the approximate solution for the deflection function of the bar body in an elastic medium was transformed into an approximate solution for the deflection function of an elastic compression-bending bar in a viscoelastic medium. Using a Kelvin body as an example for the side medium of the bar, a deflection function under a second-order triangular series was derived. MATLAB is applied in this study to compute the corresponding theoretical solutions under the burgers model.
• The influence of the magnitude of axial force on the change of w(L) over time under equal bar length and equal pile load strength is studied. The calculation results show that when there is constant load on the side of the bar, axial force acts on the top of the bar, and the larger the axial force, the faster the rate of change of w(L) with time, and the earlier it develops towards instability. The influence of pile load strength on the change of w(L) over time under equal bar length condition is studied. The results show that the change of bar length has little effect on the stability of elastic bar in viscoelastic medium, and the research should focus on the influence of axial force and pile load strength on the surface of side medium on bar stability.

This research presents the derivation of the deflection equation of an elastic rod in a viscoelastic medium and discusses the potential implications for the durability of a long rod under this assumption. However, the working condition assumed in this study has not been validated or compared with the actual engineering scenarios. In future research, more realistic engineering situations can be examined by incorporating boundary conditions, and contrasted with the actual engineering cases, to demonstrate the applicability and accuracy of this method in practical work.