Introduction

Much of the world’s oil and gas resources are transported through pipelines. In cold regions, frozen soil can cause freezing damage to buried pipelines, resulting in oil and gas leaks that have economic and environmental impacts1. Freezing damage to pipelines is mainly caused by differential frost heaving of the soil mass. The mechanical behaviors of buried pipelines subjected to frost heaving force are a research hotspot and a key concern of engineers. Some parts of a pipeline may be located in regions with frost heave, while others are not. This type of differential frost heave scenario can cause considerable stresses to pipelines, potentially breaking them2,3,4. Hence, it is important to investigate the mechanical behaviors of buried pipelines subjected to differential frost heaving. Pipelines play a vital role in energy transportation, and the development of a continuous monitoring and evaluation model for pipeline safety is a critical area of focus for future research5,6,7,8.

Many studies have investigated pipeline–soil interactions under differential frost heaving. In comparison with theoretical methods, experimental analysis represents the most direct approach9,10. In 1999, Nixon and Burgess11 studied the effects of frost heaving and melt settlement deformation on petroleum pipelines in Norman Wells. They found that pipelines in the heave regions developed 0.3% bending strain and the bending deformation of pipelines might be caused by high axial stress and differential frost heaving caused by temperature differences. Additionally, an analysis of pipeline strain indicated that about 0.2 m of frost heaving movement was enough to cause uplift buckling of the critical uplift length (22–25 m). The data findings once again validate that frost heave has a considerable impact on the safety of pipelines. Huang, et al.12 and Tanaka, et al.13carried out a field test of frost heaving characteristics and thermodynamic behaviors of large-diameter natural gas pipelines passing through a transition belt between unfrozen soils and permafrost in Fairbanks, Alaska, this monitoring persisted for approximately five years. The maximum pipeline displacement occurred in the transition region, which was attributed to frozen pipeline-induced heat gaps between unfrozen soils and permafrost. The on-site monitoring data reveals the thermodynamic behavior characteristics of pipelines under frost heave circumstances. The large-scaled field test of pipeline–soil interactions provided valuable data that are used in the present study. However, the experiments had relatively poor repeatability and operability and were time- and labor-consuming. From an economic perspective, this research endeavors to discover a straightforward and effective approach to obtain the deformation and stress outcomes of pipelines. By comparing pipelines in Caen, Rajani and Morgenstern14 designed models of two steel tubes with different radii and wall thicknesses, and lengths that were 32–40 times the diameter. The steel tubes were embedded into a test chamber and loads were applied to simulate the influences of frost heaving on deformation of the tubes. The force and deformation results of the pipeline at different times were obtained. Although the results of model experiments have certain reference value, there is no comparison with actual field observation results, which limits the accuracy. Xu, et al.15 carried out a laboratory test on a model of the Mohe-Daqing Section of the Sino-Russia crude oil pipelines. They monitored the displacement and axial strain of the pipeline as well as the stress between the pipeline and soils by controlling the environmental and oil temperatures. The results showed that soil freezing can cause pipeline uplifts and increase their axial strain, as well as the stress between the pipelines and soils. The stress–strain gauges deployed in this experiment are limited; thus, the results obtained are not the stress–strain outcomes of the entire pipeline section over time and temperature. Laboratory model tests are often greatly influenced by the test sites, devices and manual operations used, which limits the usefulness of their results.

Nixon, et al.16 pointed out that differential frost heaving can cause pipeline deformation when a frozen pipeline passes from frozen to unfrozen soils. They suggested regarding the soil mass as an elastic or nonlinear viscous continuity and applying nonlinear boundary conditions to express the frost heaving behavior of the soil. They also hypothesized that pipelines are completely passive structural components, this is inconsistent with the reality. Rajani and Morgenstern17 pointed out that Nixon did not investigate pipeline–soil interaction problems. Selvadurai and Shinde18 modeled buried pipelines as flexible round beams and considered the creeping properties of frozen soil regions and the elastic properties of unfrozen soils. They assumed that the frost heaving process and pipeline–soil interactions could be considered independently. Obviously, pipelines and soil frost heaving form a mutually restrictive coupled relationship. Therefore, it is unrealistic to view pipelines as completely passive components, or to consider pipelines and soil independently. Selvadurai, et al.19,20proposed a mechanical behavior analysis model for buried pipelines in discontinuous frozen soil regions. In this model, pipelines are viewed as beam units with axial, shearing, and rigid properties. Frozen soils are simulated with a viscoelastic constitutive model, while unfrozen soils use an elastic constitutive model. Wu, et al.21 constructed a three-dimensional model of interactions between petroleum pipelines and frost heaving mounds using ANSYS (ANSYS, Inc., 1997) finite element software. In this model, the pipelines are viewed as linear components, while the soil mass is viewed as a nonlinear medium. Pipeline yield is estimated using the Von Mises yield criterion. Similarly, Wu, et al.22 built a nonlinear elastic–plastic model for a pipeline–soil system using ANSYS (ANSYS, Inc., 2005), in which both components are viewed as elastic–plastic materials. The Mises plastic yielding criterion is employed to estimate the plastic yield of the pipelines, while the Drucker-Prager model is applied to simulate the soil mass. In these models, soils are assumed to swell at sub-zero temperatures. Although such a hypothesis provides a convenient way to study pipeline–soil interactions, it is unrealistic to simply express soil frost heaving as linearly related with temperature. This fails to reflect the genuine relationship between soil frost heave and temperature.

Buried pipelines are deformed by the action of overlying soils. Such deformation is restricted by the surrounding soil mass. Hence, soils not only generate loads on pipelines, but also act as a medium that maintains or increases pipeline strength and rigidity. Some researchers have summarized these scenarios as typical soil-structure interaction problems, which are usually solved by the Winkler theory. In other words, to analyze the stresses and strains on pipelines, they are viewed as elastic foundation beams, while the soils are viewed as a series of discontinuous springs described by elastic–plastic load–displacement equations. The Winkler foundation beam has been widely applied to theoretical models of the mechanical behaviors of buried pipelines for years. Kanie, et al.23 pointed out that the mechanical behaviors of pipelines can often be separated from the frost heaving of the soil mass when studying pipeline–soil interactions. The mechanical behaviors of pipelines are often described in terms of beam theory on an elastic foundation, while soil frost heaving is often calculated by segregation potential theory or the Clapeyron equation. To simplify the calculations, Kanie suggested estimating soil frost heaving by the Takashi equation. The deformation and bending stress of a pipeline in Calgary were predicted well by combining the Takashi equation and an elastic foundation beam. This approach is to initially predict the frost heave of the soil, estimate the corresponding frost heave force, and subsequently introduce it into the elastic foundation beam system to assess the deformation and stress of the pipeline. Nixon24, Mikkola and Hartikainen25, and Zheng, et al.26,27 predicted deformation in buried pipelines caused by frost heaving of the underlying soil by using a discrete ice lens model, thermodynamic model, and Takashi frost heaving empirical model. In these theoretical studies, soil masses and pipelines are often studied independently. However, pipelines and soil have a mutual interaction that should be considered together. Results inconsistent with the stress characteristics of some buried pipelines can be gained through independent analyses of pipelines and soil masses. With consideration of existing shortcomings in the analysis and calculation theory of the buried pipeline structures, Liu28 reviewed associated field experiments and laboratory simulation test data. A tensile stress-free Winkler model, Vlazov model, and combined model were used to analyze buried pipeline–soil interactions, and a method was proposed to determine the relevant model parameters. By analyzing the stress characteristics and deformation performance of the buried pipelines, a new method for calculating the internal stress and deformation of the round profile of pipelines was developed based on interaction analysis theory with consideration of the actual working state. However, this theory is directed at the analysis of the force and deformation of the pipeline’s cross-section, and what is of greater concern is the axial force and deformation of the pipeline. Gong, et al.29 employed the Winkler beam theory to calculate pipeline deformation and stress under overlying loads. She also analyzed the influences of the overloading position, buried depth, and rigidity of pipelines, as well as soil properties, on the internal stresses and deformation of pipelines using a case study. This study examines the interaction between pipes and soil under ordinary soil conditions, while the study of pipe-soil interaction under frozen soil conditions is more complicated.

The summary of research methods for evaluating or estimating the interaction behavior between buried pipelines and frost heave/thawing is shown in Table 1. So far, normal designs have been imperfect and require some theoretical analysis and guidance30,31,32. In theoretical analyses, the pipeline–soil interaction is generally studied using the Winkler foundation beam. In thawed soils, there are many studies on pipeline–soil interactions. In frozen regions, pipelines can still be modeled using the Winkler model (if the basic assumption is met), although with some differences: (1) frozen soils require a much different foundation coefficient to thawed soils33,34. (2) Frost heaving forces on pipelines generate new loads (additional frost heaving loads). Under such loads, buried pipelines interact with the surrounding soil mass. Due to non uniformity in frost heaving, the mechanical analysis of pipeline–soil interactions are complex due to the dependence between frost heaving-induced deformation of the pipeline and the frost heaving force. Hence, conventional foundation beam theory requires improvement.

Table 1 Summary of the research methods for assessing or estimating the interaction behaviours between a buried pipeline and frost heave/thaw.
Full size table

This study focuses on improving models based on the conventional Winkler foundation beam. A discrete forecasting model for analyzing the mechanical behavior of buried pipelines is proposed by introducing the concept of frost heaving spring. This model is used to analyze the mechanical behaviors of buried pipelines at different moments under the influences of differential frost heaving. Although frost heaving spring cannot reflect the frost heaving process, the pipeline deformation and stress conditions at each moment can be analyzed by using the frost heaving model. The Takashi empirical equation can indirectly estimate the relationship between the frost heaving amount and force according to development of freezing depth at each moment, thus obtaining the nonlinear properties of frost heaving spring. Tests of the famous Caen and Alaska pipelines were carried out and provided useful calculation results that were consistent with observed data.

Physical model for stress analysis of buried pipelines under frost heaving conditions

Unlike pipeline–soil interactions in ordinary soils, frost heaving conditions bring the additional influence of frost heaving force. How to apply it to the pipeline soil interaction system is an important issue. In the Winkler foundation beam model, Winkler simulated a foundation as a series of independent springs on a rigid base. When the foundation surface suffers a pressure of p, there is local settlement (y) at that point. It can be seen from Fig. 1 that the ordinary soil is replaced by ordinary foundation springs and frost-heaving-sensitive soil is replaced by frost heaving springs. The frost heaving force is the driving force. Imagine that the frost heaving force is provided by a kind of pre-compressed special frost heaving spring, which induces deformation of the foundation beams. The freezing process in the frost heaving section can be viewed as consistent. The stress paths at different points have mild influences on the relationship between the frost heaving force and frost heaving amount. Instead, this relationship is mainly decided by the final stress (maximum stress). Here, such a spring is defined as a physical model called a frost heaving spring. Its properties are manifested by the relationship between the frost heaving force and amount (rate). The frost heaving spring has different mechanism from an elastic foundation spring. In the beginning, a frost heaving spring is in a state of pre-compression. Due to pipeline–soil interactions, pipelines generate corresponding displacements. Displacement of a pipeline section in a frost heaving region causes pipeline sections in non-frost heaving regions to develop corresponding displacements. In this case, the soil mass in the non-frost heaving region is squeezed by the pipeline, thus developing corresponding elastic responses. Based on the above processes, the left spring is an elastic foundation spring that responds to the foundation, whereas the right spring is an imitatively bouncing spring that is released by the frost heaving force. Substantially, the mechanisms of these two springs are completely different and the springs have different properties. This is different from conventional foundation beam theory. The Winkler foundation beam theory is improved to analyze the mechanics of buried pipelines under frost heaving conditions. The frost heaving spring meets the following basic conditions:

  1. (1)

    A frost heaving section is equal to the action state of a pre-compressed frost heaving spring.

  2. (2)

    The pre-compression force of a frost heaving spring is equal to the maximum frost heaving force of the soil mass at the analysis nodes.

  3. (3)

    Upon elastic deformation of a frost heaving spring, its force on the pipeline is always equal to the frost heaving force. The elastic deformation of the spring is equal to the frost heaving deformation under the corresponding constraint of frost heaving force.

  4. (4)

    The frost heaving spring can be compressed but it cannot be stretched and it cannot be compressed again after rebounding.

Fig. 1
figure 1

Schematic diagram of physical model of buried pipelines under differential frost heave conditions.

Full size image

Construction and solving of the mathematical model

  1. (1)

    Properties of frost heaving springs and construction of the mathematical model

The properties of a physical spring are reflected by the relationship between its stresses and deformation. The properties of a frost heaving spring are reflected by the relationship between the frost heaving force and frost heaving amount (rate). This relationship has been extensively studied. Japanese scholars pointed out a higher frost heaving rate results in a higher frost heaving force. Hence, the maximum frost heaving force can be calculated under a given frost heaving rate. Some research also demonstrates that this relationship presents no linear changes. Different expressions can be obtained by fitting different models to test results35.

Liu36,37 derived the relationship between the frost heaving rate (\(\eta\)) of foundation bottom and a normal frost heaving force (\(F\)) based on tests as \(F = a\eta^{ - b}\), where \(a\) and \(b\) are constants of the soil. In a test, Huang4 proposed the concept of the constraint ratio with respect to the relationship between frost heaving force and deformation. The constraint ratio is \(\theta = \frac{{\eta_{0} - \eta_{i} }}{{\eta_{0} }} \times 100\), where \(\theta\) is the constraint ratio (%), \(\eta_{0}\) is the free frost heaving ratio under no load and no constraint (%), and \(\eta_{i}\) is the frost heaving rate under a constrained condition (%). In the experimental analysis, Huang eliminated the influences of time factors and gained a fitting formula of \(F - \theta\) under different paths. Based on foundation frost heaving tests of their soil model, the relationship between frost heaving force and frost heaving amount was disclosed. This relationship was fitted by an exponential function under three different paths.

$$F = - 12.27 \times \left( {1 - e^{0.044\theta } } \right),R^{2} = 0.9977$$
(1)
$$F = - 26.24 \times \left( {1 - e^{0.037\theta } } \right),R^{2} = 0.9811$$
(2)
$$F = - 40.50 \times \left( {1 - e^{0.032\theta } } \right),R^{2} = 0.9930$$
(3)

In this model, a frost heaving spring has nonlinear properties and may be fitted by many formulas. The key is to reflect the relationship between the actual frost heaving amount and frost heaving force. For the unification of forms and convenience of calculation, the properties of nonlinear frost heaving springs mentioned in the following text can all be calculated according to the simplest relationship between frost heaving force and frost heaving amount: \(F = P_{0} + Ae^{ - a\omega \left( x \right)}\). In this relationship, \(A\) and \(a\) are coefficients obtained from testing that are related to soil frost heaving (\(P_{0}\) is an adjustment coefficient that makes the fitted curve more closely approximate the measured values), and b (m) refers to the calculated width of the pipelines. The detailed working conditions are shown in Fig. 2.

Fig. 2
figure 2

Working condition diagram.

Full size image

The pipeline is located in the frost heave sensitive zone and non-frost heave zone (coordinate system establishment is shown in Fig. 2).

When 0 < x < L(frost heave zone) ,

$$E_{p} I_{p} \frac{{d^{4} \omega \left( x \right)}}{{dx^{4} }} = b\left( {P_{0} + Ae^{ - a\omega \left( x \right)} } \right)$$
(4)

When x < 0(non frost heave zone) ,

$$E_{p} I_{p} \frac{{d^{4} \omega \left( x \right)}}{{dx^{4} }} = - k_{2} \omega \left( x \right)$$
(5)
  1. (2)

    Solving of the pipeline mechanical model.

The analytical solution is as follows:

In order to facilitate the solution, they are defined as the characteristic coefficients \(\beta\) as:

$$\beta = \sqrt[4]{{\frac{{k_{2} }}{{4E_{p} I_{p} }}}}$$
(6)

The units of \(\beta\) are \(\frac{1}{{\text{m}}}\), then define \(\left[ L \right]_{\beta } = \frac{1}{\beta }\) as the characteristic lengths, the units of \(\left[ L \right]_{\beta }\) are \({\text{m}}\).

When x < 0, the general solution is:

$$y = e^{\beta x} \left( {A\cos \beta x + B\sin \beta x} \right) + e^{ - \beta x} \left( {C\cos \beta x + D\sin \beta x} \right)$$
(7)

Now there are 4 coefficients \(A,B,C,D\), which are obtained according to boundary conditions and continuity conditions.

Since the parameters of the frost heaving region are expressed by an exponential function, it is a boundary problem described by a fourth-order ordinary differential equation. The solving process is too complicated to obtain an analytic solution. Instead, it can only provide numerical solutions through calculation. Since the boundary calculation of the fourth-order ordinary differential equation is relatively complicated, one mathematical software with the bvp5c function was used. The bvp5c function calculates numerical solutions for two-point boundary problems. The basic principle is to convert the boundary problems into initial-value problems through the target dichotomy, and then calculate the initial values by the Runge–Kutta method. Specifically, it has to find numerical solutions for the following n first-order ordinary differential equations:

$$\frac{{dy_{j} }}{dx} = f_{j} \left( {x,y_{1} ,y_{2} ...y_{n} ,q} \right),j = 1,2,3...n$$
(8)

where x has the range of \(a \le x \le b\). The initial conditions in the interval are: \(y_{j} \left( a \right) = a_{j}\) and \(y_{j} \left( b \right) = b_{j}\), where \(j = 1,2,3...n\), \(a_{j} ,b_{j}\) are constants.

In finding numerical solutions with the mathematical software, boundary conditions must be defined. As shown in Fig. 3, the bvp5c function is employed to solve the boundary problems of this fourth-order ordinary differential equation.

Fig. 3
figure 3

Bvp5c operation function solving process.

Full size image

The boundary conditions and continuity conditions listed are as follows (a total of 8 conditions):

At that time \(x = - \infty\),\(y = 0\),\(\theta = 0\);

From the continuity condition of points, we can obtain:

When the two solutions are equal at point \(x = 0\), and the 1st to 3rd derivatives of the solutions are equal at point \(x = 0\).

When \(x = L\) (the point with the highest frost heave in the middle), the corner \(\theta = 0\) and the shear force \(Q = 0\).

Model discretization

The description of the frost heaving spring provided above only considers the stress conditions of pipelines at the default end of frost heaving. The relationship between the frost heaving force and frost heaving amount from the beginning to the defaulted end of frost heaving is viewed in terms of the properties of the frost heaving spring. Although the frost heaving spring cannot reflect the frost heaving process, it provides information on pipeline deformation and stress conditions at any moment when used in the frost heaving model. In this study, the mentioned spring state refers to the different states of the frost heaving springs of a particular row under consistent freezing conditions. The free frost heaving amount is ω0. This implies that the freezing process calculation is for the node where the free frost heaving amount is ω0 (known as the calculation point, such as the year when the maximum frost heaving occurred, the calculation point in the freezing process, or the end-point of the experiment). This corresponds to the concept of there being a maximum frost heaving force at each calculation node. In other words, if there is a complete constraint against the soil mass (frost heave = 0) in this freezing process, the soil mass will generate the maximum frost heaving force to the constraint. This corresponds to the frost heaving spring; that is, the pre-compressed constraint state. During calculation of the mechanical behaviors at different nodes (e.g. 1 month or 2 months after freezing, or when the frost heaving amounts are 1/4ω, 1/2ω, and ω, where ω is the frost heaving amount), the corresponding maximum frost heaving forces vary and the frost heaving springs are different. However, a row of frost heaving springs in the frost heaving section have the same values in the calculation of each node due to the consistent frost heaving performance and consistent freezing conditions. The frost heaving springs represent the relationship between the frost heaving amount and frost heaving force at the calculation node. The frost heaving section uses a row of such springs. It corresponds to another node and has another row of frost heaving springs with the same performance. The results of each node are connected to represent the discretization process.

There are \(n\) stages in Fig. 4. Specifically, Stage 1 occurs from \(t = 0\) to \(t = t_{1}\), Stage 2 from \(t = 0\) to \(t = t_{2}\), …… and Stage \(n\) from \(t = 0\) to \(t = t_{n}\). It is important to note that these stages are hypothesized as independent processes rather than as a single continuously changing process. The end-point of each stage is viewed as the final moment. The frost heaving spring is analyzed at each stage. In this case, the relationship between the frost heaving force \(y\) and frost heaving amount \(\omega \left( x \right)\) is \(y = F\left( {\omega \left( x \right)} \right)\). In this way, expressions of the properties of the nonlinear frost heaving springs corresponding to each stage can be written as Eqs. (9–11). From the perspective of a physical model, the frost heaving springs in each stage can be viewed as having different frost heaving properties. For example, the pipeline deformation and stress conditions at \(t_{1}\), \(t_{2}\), and \(t_{n}\) can be calculated by using properties of the springs at Stages 1, 2, and n, respectively. The function \(F\left( {\omega \left( x \right)} \right)\) involves parameters related to frost heaving, which can be determined by the curve derived from fitting the frost heaving model. Figure 4a shows the fitting curve of the frost heave force and frost heave amount corresponding to stages 1 to n. Figure 4b shows the pipe soil interaction before frost heave, with (b) to (d) representing stage 1, (b) to (e) representing stage 2, and (b) to (f) representing stage n. Figure 4c shows the physical meaning represented by the springs in the figure.

$$Phase{\mkern 1mu} 1{\text{(}}t = 0\;to\;t = t_{1} ):y = F_{1} \left( {\omega \left( x \right)} \right)$$
(9)
$$Phase{\mkern 1mu} 2{\text{(}}t = 0\;to\;t = t_{2} )y = F_{2} \left( {\omega \left( x \right)} \right):$$
(10)
$$Phase{\mkern 1mu} n{\text{(}}t = 0\;to\;t = t_{n} ):y = F_{n} \left( {\omega \left( x \right)} \right)$$
(11)
Fig. 4
figure 4

Schematic diagram of pipe-soil interaction at the stages1 to n.

Full size image

The assumption that nonlinear frost heave springs with different properties are constantly replaced over time has supplemented and improved the research on the stress state of pipelines at any moment. The properties of nonlinear frost heaving springs at each moment can be estimated through the frost heaving model, thus enabling the analysis of pipeline deformation and stress conditions at each moment. The Takashi equation can solve this problem by approximating the relationship between frost heaving force and amount according to the frost heaving rate and the relationship between the constraint stress along the freezing direction and the freezing rate. Based on comprehensive one-dimensional laboratory tests, Takashi, et al.38 from Japan proposed a practical formula to estimate frost heaving. In the Takashi equation, the freezing rate and effective constrained stress in the freezing direction are major factors for predicting frost heaving. Three constants represent the soil mass properties. The Takashi equation has been applied by the Japan Geotechnical Society (JGS); therefore, there are regulations for using it in practical engineering. Moreover, the Takashi equation has been extensively applied in practical projects in Japan, with good effects. The Takashi equation was used as the basis for frost heaving prediction in the present study. It connects the frost heaving rate with the constrained stress in the freezing direction and freezing rate. As shown in Eq. (12) and highlighted by Takashi, et al.38, the frost heaving rate is correlated with the constrained stress in the freezing direction and freezing rate:

$$\xi = \xi_{0} + \frac{{\sigma_{0} }}{\sigma }\left( {1 + \sqrt {\frac{{U_{0} }}{U}} } \right)$$
(12)

where \(\xi \left( \% \right)\) is frost heaving rate and (defined as the ratio of the volumetric growth of samples to their initial volume), \(\sigma \left( {{\text{kPa}}} \right)\) is the effective constraint stress in the freezing direction, and \(U\left( {{{\text{m}} \mathord{\left/ {\vphantom {{\text{m}} {\text{s}}}} \right. \kern-0pt} {\text{s}}}} \right)\) is the freezing rate. Variables \(\xi_{0}\), \(\sigma_{0}\) and \(U_{0}\) are representative constants of soil properties gained via the JGS’s adjustment experiment (JGS 0171-2003)39, which provided their empirical values. To obtain the aforementioned parameters, the Takashi method is employed for frost heave analysis, which demands at least thermal analysis and the acquisition of corresponding parameters through experiments as stipulated in JGS. The relationship between frost heaving rate and effective constraint stress is shown in Eq. (12), while the frost heaving rate can be expressed as the ratio of frost heaving amount to freezing depth. Based on Eq. (12), the relationship between the frost heaving amount and constraint force (frost heaving force) can be acquired and fitted by the exponential function \(F = P_{0} + Ae^{ - a\omega \left( x \right)}\), thus obtaining the properties of the nonlinear frost heaving springs. The parameters of the Takashi equation can be gained by looking up norms.

Case study

Two cases were used for verification: Case 1 is the Caen pipeline test and Case 2 is the Alaska pipeline test.

Case 1. The Caen pipeline test was a large laboratory test carried out in France and Canada. The test devices included a large refrigerator (length = 18 m, width = 8 m, height = 5 m) with controlled temperature and a groove with a filling depth of 2 m. At the bottom of the groove, the internal heating power and hydraulic power systems were separated from the natural ground conditions. The testing groove was filled with soil to the depth of 1.75 m. This test was conducted to compare pipeline deformation and stresses under the influence of soil masses with different frost heaving sensitivities. Hence, two types of soils were used. The frost heave-sensitive soil in Caen was composed of 13–20% clay, 65–75% silt, and 10–20% sandy soil. The non-sensitive soil was SNEC sandy soil, which was composed of < 10% silt and 80–90% sand. The Caen silt and SNEC sand accounted for half the volume of the test groove. The test pipeline was composed of steel tubes (length = 18 m, diameter = 273 mm, wall thickness = 0.5 mm). This pipeline had an independent cooling system. The test instrumentation was buried in an 18 m-long groove at a depth of 33 cm, and then covered with the same soils. For long-distance pipeline modeling, there was no constraint at the two ends of the pipeline19,20,40.

Case 2. Frost heaving is a major design problem in natural gas pipeline projects along highways in Alaska. This project involved the design, construction, and operation of a large-diameter natural gas pipeline running through Alaska from the north slope, across the Yukon Region and Western Canada. It transports natural gas to downstream markets in the USA. In the continuous frozen soil region and some discontinuous frozen soil regions, the natural gas is cooled to < 0 ℃. The cooled gas flows from Prudhoe Bay on the north slope of Alaska to the Yukon region. If the above mentioned problems are not fully considered in the design and operation of the pipeline system, a frozen pipeline which runs through non-frozen regions in discontinuous permafrost areas may suffer serious frost heaving problems during operation. Huang, et al.12 carried out an outdoor test in Fairbanks, Alaska, from 1999 to 2003 using a large-diameter pipeline for cold gas transport. In the freezing stage in 1999, large deformations and stresses of the pipeline, as well as the temperature distributions underground and in the pipeline, were observed. The pipeline was buried at the permafrost boundary. Cold air (− 10 ℃) was circulated in the pipeline. A frost layer only encased the pipeline in non-permafrost areas. Field grading analysis of 16 soil samples found them to be generally composed of 7.9% sand, 78.1% silt, and 14.0% clay, 0.02% of the particle size is finer than 41%. The detailed parameters of the Caen and Alaska pipelines are listed in Table 212,41.

Table 2 Assumed calculation parameters.
Full size table

The development of freezing depth in the regions of the Caen and Alaska pipelines is shown in Fig. 5a. There were single freezing-depth monitoring points at the south and north ends of the of Alaska pipeline. Clearly, the freezing depths at the southern and northern monitoring points differ slightly. The black line is the calculated freezing depth reported by Kim, et al.42. It differs slightly from the observed values. Both can provide references for the following calculation.

Fig. 5
figure 5

(a) Development curve of freezing depth (b) Fitting curve of frost heave amount and frost heave force of Caen pipeline (c) Fitting curve of frost heave amount and frost heave force of Alaska pipeline.

Full size image

Case 1: Caen test.

According to Eq. (12), the connection between frost heave and overlying load could be acquired. For situation 1 (Fig. 5b), to obtain the relationship between frost heave and overlying load at 50 d, 100 d, 150 d, 250 d, and 350 d, the corresponding freezing depth at each time point were determined based on Fig. 5a. This case constitutes a field experiment executed by previous researchers. It was manifestly impracticable for this study to obtain the corresponding soil parameters anew for theoretical calculation outcomes. Thus, this manuscript made reference to some literatures. It was existed precedents in these literatures that employ Takashi’s empirical equation to analyze this case and have attained auspicious validation results within the theoretical method proposed in the paper. As a consequence, these parameter values were reliable. Based on reference literature and trial calculations, the parameters were set as \(\xi_{0} = 3.5 \times 10^{ - 2}\),\(\sigma_{0} = 3{\text{kPa}}\),\(U_{0} = 6.9 \times 10^{ - 7} {\text{m/s}}\), and \(U = 5.6 \times 10^{ - 7} {\text{m/s}}\)27. By substituting the freezing depth and the aforementioned parameters into Eq. (12), the characteristics of the hypothesized frost heave spring corresponding to different time periods could be derived, namely the correspondence between the frost heave force and the frost heave amount. With the aim of enhancing the accuracy of the fitting curve, the magnitude of the free frost heave was set as \(\omega_{0}\) = 0.085 m, 0.12 m, 0.14 m, and 0.18 m, respectively19,20,40. The relationship between load and frost heaving amount can be expressed as the general formula previously mentioned \(y = p_{0} + A \cdot e^{ - ax}\). The x-axis reflects the frost heaving amount, while the y-axis shows the overlying loads. The values of parameters \(\omega_{0}\), \(p_{0}\), \(A\), and \(a\) of the fitted curve between frost heaving amount and overlying load as shown in Fig. 5b are listed in Table 3.

Table 3 Assumed calculation parameters.
Full size table

Case 2: Alaska test.

The freezing depth in Case 2 are shown in Fig. 5a. With reference to the literature with appropriate adjustment, the parameters were set as \(\xi_{0} = 1.5 \times 10^{ - 2}\), \(\sigma_{0} = 1.2{\text{kPa}}\),\(U_{0} = 6.9 \times 10^{ - 7} {\text{m/s}}\) and \(U = 5.6 \times 10^{ - 7} {\text{m/s}}\)27. Similarly, by substituting the freezing depth and the aforementioned parameters into Eq. (12), the relationship between frost heave and overlying load at 377 d, 512 d, and 754 d were derived. With the aim of enhancing the accuracy of the fitting curve, In this case, \(\omega_{0}\) was set as 0.09 m, 0.12 m, and 0.14 m, respectively12,41. Clearly, the frost heaving amount decreases exponentially with increases in overlying load. Under this circumstance, the relationship between frost heaving amount and overlying load can be expressed as \(y = p_{0} + A \cdot e^{ - ax}\). The x-axis reflects the frost heaving amount, while the y-axis shows the overlying load. The values of the parameters \(\omega_{0}\), \(p_{0}\), \(A\), and \(a\) of the relationship between frost heaving amount and overlying load as shown in Fig. 5c are listed in Table 3. After determining the relationship between frost heave force and frost heave amount at each instant, mathematical model formulas (4) and (5) would be employed to solve and acquire the deformation and stress results of the pipeline.

According to the soil properties in unfrozen soil regions in Case 1, it was determined that \(k_{2} = 1.63{\text{MPa}}\) according to the foundation coefficients and pipeline diameter43. The curves of pipeline deflection, rotation angle, bending moment, shear force, and deformation at 50 d, 100 d, 150 d, 250 d, and 350 d were calculated by the discrete forecast model and are compared with the measured results in Fig. 6a. A comparison between the calculated and measured bending stresses is also shown in Fig. 6a. According to the calculation results, the pipeline deflection, rotation angle, bending moment, shear force, and loading peak all increased with time. In Fig. 6a, the maximum pipeline deflection occurred in the frost heaving region on the left of the pipeline. The deflection curve increased gradually from 50 to 350 d. The maximum deflection of the pipeline at 50 d was 0.079 m and the maximum deflection at 350 m was 0.18 m. In the coordinate system, the maximum rotation angle of the pipeline occurred at 8 m of the x-axis; that is, the middle position of the pipeline. The rotation angle curve increased gradually from 50 to 250 d. The maximum rotation angle of the pipeline at 50 d was 0.015, which increased to 0.03 at 350 d. The bending moment curve of the pipeline also increased gradually from 50 to 350 d. The maximum bending moment of the pipeline at 50 d occurred at 6.5 m of the x-axis, being 35 \({\text{kN}} \cdot {\text{m}}\). The minimum bending moment of the pipeline occurred at 10 m of the x-axis, with a value of − 27 \({\text{kN}} \cdot {\text{m}}\). The maximum bending moment (55 \({\text{kN}} \cdot {\text{m}}\)) of the pipeline at 350 d occurred at 6.5 m of the x-axis, while the minimum bending moment (− 55 \({\text{kN}} \cdot {\text{m}}\)) occurred at 10 m. The shear force of the pipeline peaked in the middle. The maximum force at 50 d was 46 kN and the maximum at 350 d was 77 kN. As time went on, the freezing depth increased gradually. As a result, the frost heaving amount of the soil mass increased gradually, accompanied by gradual increases in pipeline deformation and stresses. Additionally, the variations in calculated deformation from 50 to 350 d were similar to the measured results, peaking in the frost heaving region on the left of the pipeline. The calculated maximum displacement of the pipeline at 50 d was 0.079 m and the measured maximum displacement was 0.09 m. At 350 d, the calculated maximum displacement of the pipeline was 0.18 m and the measured maximum displacement was 0.189 m. The calculated deformations were largely similar to the measured results. The calculated bending stresses from 50 to 250 d fluctuated within a smaller range than the measured results. The calculated maximum bending stress (96 kPa) of the pipeline at 50 d was at the 10 m of x-axis, and the measured result was 95 kPa. The calculated minimum bending stress (− 78 kPa) of the pipeline at 50 d was at the 6.5 m of the x-axis, and the measured result was − 77 kPa. The calculated maximum bending stress (200 kPa) of the pipeline at 350 d was at 10 m of the x-axis, and the measured result was 189 kPa. The calculated minimum bending stress (− 201 kPa) of the pipeline at 350 d was at 6.5 m of the x-axis, and the measured result was − 186 kPa. Generally speaking, the calculated results did not deviate, greatly from the measured ones and were relatively robust. To sum up, the calculated pipeline deformation and bending stress curves show similar trends to those of the measured results, this attests to the rationality and feasibility of the model. The calculated results are in a reasonable error range and are relatively robust, showing similar trends and high agreement with the measured results.

Fig. 6
figure 6

(a) Pipeline mechanics analysis results of deflection, rotation angle, bending moment, shear force, and comparison of pipeline stress and displacement with measured values of Case1; (b) Comparison between predicted and measured displacement changes of three test points over time of Case2; (c) Pipeline mechanics analysis results of deflection, rotation angle, bending moment, shear force, and comparison of pipeline stress and displacement with measured values of Case2.

Full size image

In Fig. 6b, pipeline displacement data were acquired from 28 lifters (HR1–HR28) to monitor pipeline movement in Case 2. Lifter HR-1 was placed at a vertical distance of 8.5 m from the inlet and HR-28 was placed at 104.5 m. Many lifters were placed near the hypothesized transition region to acquire greater accuracy there. A total of 25 lifters were welded directly onto the top external surface of the pipeline between 0 and 58 m of the vertical pipe. The measured displacement of HR-1 was < 0 before 500 d, indicating that the pipeline initially subsided to some extent. The displacement at 600 d was increased to 0.02 m and then remained largely constant. The variation in the calculated displacement increased gradually to about 0.008 m and showed a slight difference to the measured results. However, the displacement gap was small, which might be attributed to errors in displacement data acquisition. The measured displacement of HR-25 at 58.805 m to the inlet increased gradually, reaching the maximum (0.144 m) at about 1050 d. The calculation results were slightly higher than the measurement results before 335 d, but they were slightly lower than the measurement results after. The calculated maximum displacement (0.14 m) was observed at about 1172 d. The measured displacement of HR-22 at 46.615 m to the inlet increased gradually, reaching the maximum (0.185 m) at about 1050 d. The calculation results were slightly higher than the measurement results before 418 d, but they were slightly lower than the measurement results after. The calculated maximum displacement (0.178 m) was observed at about 1172 d. Upon comparing the calculated results with the measured values, it can be observed that there is a minor fluctuation in the displayed curve of the measured values, while the calculated curve is relatively smooth. The fluctuation in the measured values might be attributed to data collection errors. Nevertheless, the trend of the displacement change curve over time calculated from the three observation points is close to the measured values, and the numerical error of the displacement calculation is not significant, which substantiates the rationality of the calculation.

According to the soil properties in unfrozen soil regions in Case 2, it was determined that k2 = 2.6 MPa according to the foundation coefficients and pipeline diameter43. Figure 6c compares the curves of pipeline deflection, rotation angle, bending moment, shear force, and deformation at 377 d, 517 d, and 754 d, as calculated by the discrete forecast model and measured results. The calculated bending stress is also shown in Fig. 6c. The maximum rotation angle of the pipeline occurred at 0 m of the x-axis; that is, the border between the frost heaving and non-frost heaving zones. The curve of rotation angle increased gradually from 377 to 754 d. The maximum rotation angle of the pipeline at 377 d was 0.007, and it increased to 0.01 at 754 d. The curve of the pipeline bending moment also increased gradually from 377 to 754 d. The maximum bending moment (558 kN·m) of the pipeline at 377 d occurred at 3 m of the x-axis, while the minimum (− 453 kN·m) occurred at -4.46 m of the x-axis. At 754 d, the maximum bending moment (726 kN·m) of the pipeline was at 3 m of the x-axis, while the minimum (-675 kN·m) was at − 4.46 m of the x-axis. At the border between the frost heaving and non-frost heaving zones, the maximum shear force of the pipeline at 377 d was 345 \({\text{kN}}\), which increased to 467 kN at 754 d. The calculated and measured maximum displacements of the Alaskan pipelines were 0.078 m and 0.0727 m, respectively, at 377 d, 0.11 m and 0.1029 m at 517 d, and 0.124 m and 0.1225 m at 754 d. The errors between the calculated and measured results were generally < 10%, indicating that the two datasets are similar. The calculated maximum and minimum bending stresses of the pipeline were 79 kPa and − 97 kPa, respectively, at 377 d, 108 kPa and − 122 kPa at 517 d, and 11 kPa and -127 kPa at 754 d, respectively.

Based on the results of the two case studies, the modeled and measured pipeline deformation curves are similar, as are the results for bending stress. This substantiates that the proposed theoretical model is cogent and dependable. In the aforementioned two cases, the displacement and stress outcomes of the pipeline have an error of approximately 10%. Although they are not precisely identical to the measured values, the measured data may also entail certain errors due to factors such as testing apparatuses and human collection. The aim of this study is to furnish a reference for pipeline design and disease prevention, and such an error range is acceptable for engineering design.

In Case 2, the pipeline outer diameter (\(b = 0.5{\text{m}}\), \(b = 0.914{\text{m}}\) and \(b = 1.5{\text{m}}\)) and wall thickness (\(h = 0.005{\text{m}}\), \(h = 0.009{\text{m}}\) and \(h = 0.015{\text{m}}\)) of the pipeline were chosen for calculation. Since the calculations of pipeline deflection are slightly influenced by diameter and wall thickness, a diagram of the variations in calculated deflection, bending moment, shear force, and bending stress peak for pipelines with different diameters and wall thicknesses is presented by Fig. 7. Clearly, the pipeline’s maximum rotation angle is negatively related with its diameter and positively related with its wall thickness. The positive bending moment peak of the pipeline shows positive correlations with both outer diameter and wall thickness. The absolute value of the negative bending moment peak also presents positive correlations with both outer diameter and wall thickness. The shear force of the pipeline increases with both outer diameter and wall thickness. The bending stress of the pipeline decreases with increases in outer diameter (slightly) and wall thickness. In summary, larger diameters and wall thicknesses lower the bending stress of the pipeline, making it more damage-resistant. The theoretical model in this study can be employed for pre-analyzing and pre-calculating the stress and deformation of the pipeline at diverse times during the pipeline design stage. The weakest link of the pipeline can be distinctly identified, and subsequently, the weakest link can be tackled in advance, such as foundation treatment, enhanced insulation measures, or modification of pipeline size. As analyzed formerly, pipelines with larger diameters and wall thicknesses are less susceptible to damage. However, choosing appropriate pipeline dimensions will be a trade-off between material usage (cost) and safety.

Fig. 7
figure 7

Calculation results of peak points of rotation angle, bending moment, shear force, and bending stress.

Full size image

Conclusions

In this study, a discrete forecast model of the mechanical behaviors of pipelines buried under differential frost heaving conditions was proposed. Two case studies were made to apply the model and the calculation results were discussed. Some major conclusions can be drawn:

  1. (1)

    The concept of frost heaving spring is proposed based on the Winkler foundation beam. The frost heaving behaviors of soil masses are connected with the stress and deformation of pipelines. The properties of the frost heaving spring depend on the relationship between frost heaving force and frost heaving amount. A mathematical model for analyzing the mechanical behavior of buried pipelines was constructed for two working conditions and solved using the mathematical software.

  2. (2)

    Pipeline deformation and stress conditions at each moment were analyzed by introducing the frost heaving model. The end-point of each stage was used as the final moment to analyze the frost heaving springs of the current stage. In other words, the frost heaving springs in each stage are viewed as springs with different frost heaving properties. The approach of varying the properties of nonlinear frost heaving springs is advancing research on pipeline stress states at a given moment. Later, the Takashi empirical equation was introduced to the frost heaving model. The model can indirectly estimate the properties of nonlinear frost heaving springs. The frost heaving springs in each stage are analyzed by viewing the end-point as the final moment. This supplements the pipeline deformation and stress condition information at each moment, thus achieving dynamic prediction of the mechanical behaviors of buried pipelines. The famous Caen and Alaska pipelines were employed as case studies. The results show that the freezing depth increases gradually, thus gradually increasing the soil frost heaving amount and intensifying pipeline deformation and stress. The calculated variations in pipeline deformation and bending stress show similar trends to the measured curves. The calculation results are within a reasonable error range and agree well with the measured results. The discretized mechanical model of the buried pipeline is effective.

  3. (3)

    According to a sensitivity analysis of the modeling results with different pipeline sizes, a larger diameter and wall thickness lower the bending stress of the pipeline, thus making it more damage-resistant. However, large diameters and wall thicknesses require more materials, making it important to balance strength with safety.