Cylindrical shear waves in soil around underground pipelines

A method of numerical solution of one-dimensional problem of cylindrical shear wave propagation in an elastic and elastic-plastic soil is developed in the article using the finite difference method. The numerical results obtained are presented in the form of graphs. From the results obtained, the attenuation of the parameters (shear stress, shear strain, and angular velocity) of cylindrical wave propagation with distance in an elastic and elastic-plastic soil was determined. The attenuation of waves with distance is justified by the dissipation of strain energy on the expanding cylindrical soil layer. In the case of load exceeding the elastic limit, plastic strains occur in soil near the point of load application. The boundaries of elastic-plastic strain of soil are determined.


Introduction
The problems of extracting piles, underground pipelines or drill strings from soil by setting torsional loads (moments) arise in the construction, drilling of oil and gas fields. Such problems, in some cases, are solved by reducing the forces of interaction and cohesion between the soil and the objects (piles, pipelines, drill strings) by the rotational motion of the latter. In this case, shear stresses and strains appear at the interaction boundary, and shear waves propagate in soil. The propagation pattern of such waves can affect other objects embedded in soil. Therefore, the study of the shear wave propagation in soil can be considered one of the urgent problems in soil dynamics. Wave propagation in various media was studied in [1][2][3][4][5]. Plane, cylindrical and spherical waves were studied in accordance with the geometry of wave propagation (the shape of the wave front). One-dimensional cylindrical waves were theoretically considered in [2][3][6][7] using the method of characteristics for elastic and elastoplastic media [6 -7]. In [8][9][10][11][12], the wave propagation in a two-dimensional formulation in elastic and elastoplastic media was studied. The reasons for the shear wave initiation and its propagation pattern in the classical case, i.e. in elastic media are reflected in [13][14]. The absorption of longitudinal waves occurs when the plastic or viscous properties of the medium are considered [2][3][15][16][17][18]. Taking into account the viscous properties, the values of the wave amplitude relative to the distance decrease, the maximum values of stress and strain are achieved at different time points due to the rheological properties of soil. When studying cylindrical wave propagation in soils [2,6], under a set shock load, the characteristic lines are straight, and the front is a shock front, and the wave parameters have a discontinuity. Therefore, when such conditions are set on the boundary, the solution to the problems by the method of characteristics is the most suited solution. If an arbitrary load is set at the boundaries, then it is more difficult to construct solutions on the characteristic lines. This article is devoted to the study of a one-dimensional elastic and elastoplastic cylindrical shear wave propagation in soils. The studies the closest to this theme were conducted in [6,18]. In [5], a two-dimensional elastic problem was considered; based on the propagation pattern of cylindrical waves and reflection from the day surface, the possibility of accidents occurring in underground pipelines was investigated [18]. In [6], waves in an elastoplastic medium were investigated by the method of characteristics, and in [18] some errors in the formulation and solution of the problem were noted [19]. The purpose of this study is to develop a solution technique that allows applying arbitrary loads on the boundary and complex strain properties, as well as determining the parameters of propagating cylindrical waves in elastic and elastoplastic soil.

Statement of the problem
Let us assume that there is a long and rigidly fixed underground pipeline with radius = 0 in an infinite soil medium. Let the pipeline begin rotational motion about its axis from the time reference point, and the pipeline strain is ignored due to the differences in the rigidity characteristics of the pipeline and soil. Then, cylindrical shear waves begin to propagate in the soil medium, the parameters of these waves are axisymmetric relative to the axis of the cylindrical body -a pipeline, and they depend only on the radial coordinate and time, i.e. the problem is one-dimensional. Thus, we consider a one-dimensional problem in cylindrical coordinates. The equation of motion of the soil medium in the absence of mass forces in the Euler representation has the form: where is the radial coordinate, = ( , ) is the velocity of soil particles along the angular coordinate, = ( , ) is the shear stress of soil along the radial coordinate. To solve equation (1), where ( , ) and ( , ) are unknown, we introduce the Cauchy relation in the following form To obtain a closed system of equations, it is necessary to add an equation of state to (1) and (2). A specific ratio is taken as an equation of state depending on the consideration of the strain properties. The general form of the equation of state of soil is: Thus, the closed system of equations consists of equations (1) -(2) and a specific form of equation (3) with respect to the unknowns ( , ), ( , ) and ( , ). To solve the problem of the cylindrical shear wave propagation, we take the following initial conditions: for = 0, > 0 : and boundary conditions: for = 0 , ≥ 0: The problem posed is numerically solved by the method of finite differences. The numerical solution is aimed for the future, bearing in mind the solution of similar problems for soil medium with complex properties (structural change, moisture content, viscoplasticity, etc.).

Method for solving the problem posed
Let us compose an algorithm for the numerical solution of the problem. For that, an infinite soil is assumed to be bounded by radius = . In this case, considering velocity of transverse wave propagation, the solution to the problem under study is obtained before the shear waves reach the boundary = , i.e., at time point ≤ ( − 0 )⁄ . The radial segment − 0 is divided into N infinitesimal annular parts, i.e. cylindrical cells. The main focus will be on the soil behavior around the underground pipeline, therefore, the finer cells are located near the pipe; with the distance from it the size of the cells increases. At time point = +1/2 , the particle velocities at the nodal points of the cells = > 0 are denoted by ( , +1/2 ) = ( ) +1/2 . The tangential stress and shear strain determined at the centers of the cells = +1/2 > 0 at the time point = are denoted by ( +1/2 , ) = ( ) +1/2 and ( +1/2 , ) = ( ) +1/2 . Thus, the input of the parameters provides the second order of accuracy of the finite-difference relations [20]. Let us assume that the values of all parameters of the problem are known up to a certain time point = . The same parameters are determined at the next time steps. Using the finite-difference scheme [20], the particle velocities for the time point = +1/2 are determined from equation (1) at the nodal points of the cylindrical cell: where ∆ -is the time step; } .
If the particle velocities are set at the boundary = 0 , ≥ 0, then for = 0 , i.e. for = 0 , instead of relation (6), we set the following velocity at the boundary: in the case of setting the shear stress at the boundary, to calculate the relation (6) for k=0, the following parameter values are used: If the angular displacement of the underground pipeline is specified at the boundary, then the particle velocity at the boundary is determined by the following formula After determining the velocities of the particles, the tangential displacement at time = +1 can be determined by the following formula: To determine the shear strain, first, we find the strain rates by approximating Eq. (2) by the finitedifference relation, where +1/2 = 1 2 ( +1 + ). Then, the shear strain is determined by the following formula Formulas (11) and (12) give the values of shear strain rates and shear strain determined in the center of the cylindrical cell at times = +1/2 and = +1 , respectively. After the shear strain and strain  (3): A pseudo-viscous term +1/2 +1/2 was added to relation (13), to perform a through calculation at the wave front, i.e., to arrive to continuous solutions and to smooth the oscillations in the numerical solution [20]: where 0 = const ≈ 2.
The time step is chosen in accordance with the stability conditions of the difference scheme [20]: In the calculation, the time step can be increased; in such cases the step increase is limited to no more than 10%, that is, if ∆ +3/2 > 1.1∆ +1/2 , then ∆ +3/2 = 1.1∆ +1/2 . Once the step ∆ +3/2 is determined, the time step ∆ +1 can be determined as ∆ +1 = 1 2 (∆ +3/2 + ∆ +1/2 ). Thus, using certain values of velocity, strain and stress at the previous time points, we calculate the values of these parameters for the next time point. By performing such actions in a sequential manner, it is possible to determine the parameters of a wave propagating in the medium (soil) until the required time or until the time point ≤ ( − 0 )⁄ . A program for solving the problem was built and implemented on a computer.

Numerical results and their analysis
Let us present the numerical results in the form of graphs. The solution was obtained for the following initial data: initial soil density ρ 0 =2000 kg•m -3 ; the velocities of longitudinal and transverse wave propagation =2000 ms -1 and =1000 ms -1 , respectively; radius of the cylindrical cavity 0 =1 m.

Elastic cylindrical waves
Let the soil medium be elastic, then the specific form of the constitutive relation (3) is = * + ( − * ) or ̇=̇ , where * , * are the initial (reference) values of shear stress and shear strain. In this case, the finite-difference relation (13) takes the form: Let the particle velocities (5) be given at the boundary in the form 0 ( ) = sin( ), where =0.2 ms -1 , is the frequency of the impact. Figures 1, 2 show changes in particle velocity and shear stress in fixed cylindrical soil layers at different frequencies of the impact. Curves 1-5 correspond to sections r=1.1, 2, 3, 4 and 6 m. Graphs (a) and (b) in figures 1-2 were obtained for = 2 ( − ) and =50 rad•s -1 , T = 0.02 s is the time of the load application, H is the Heaviside function. As seen from figures 1 and 2, in the case of setting the tangential velocity at the boundary = 0 , the amplitude of the particle velocity and shear stress decreases with distance. This decrease (damping) occurs at the first arrival of the wave, then damping in time is not observed in fixed sections. At a positive value of the particle velocity, i.e. under motion (rotation) in one direction, the tangential stress increases, and under reverse motion (at a negative value of the particle velocity), the stresses decrease (the stress state is unloaded). An increase in the frequency of the set particle velocity at the boundary shortens the period of "loading-unloading" oscillations; this indicates a decrease in the maximum values of shear stresses. At that, the frequency and periodicity of the propagating wave do not change.      The presented patterns of changes in particle velocities and shear stresses are similar to figures 1 and 3. Here, wave attenuation with distance from the initial cylindrical section is observed. At the first arrival of the wave, the particle velocity increases abruptly to its maximum value, and then changes over time similar to the previous result. In general, under the elastic strain of soil, cylindrical wave attenuates with distance. Energy dissipation occurs due to its redistribution on the expanding cylindrical layer, which causes wave attenuation.

Elastic-plastic waves
Let the soil under consideration obey the elastic-plastic law of strain (3): loading in the elastic zone and unloading occur according to the following law ̇=̇, and loading in the plastic zone occurs according to ̇= 1̇ for τ ≥ τ S , with τ S varying to the time of the beginning of unloading from the plastic zone for each point. Here 1 -is the shear modulus of plastic strain (coefficient of proportionality between shear stress and shear strain). The finite difference relation according to (13) has the form Let us check the plasticity conditions. If the condition ( * ) +1/2 +1 < , is satisfied, soil particles are strained in the elastic zone and the values of shear stresses are (τ ) k+1/2 n+1 = ( * ) +1/2 +1 . Otherwise, (τ ) k+1/2 n+1 = τ S , = ( ) +1/2 + 1 (̇) +1/2 +1/2 ∆ +1/2 and soil particles are plastically strained.
Here it is necessary to clarify the application of the equation of state in terms of increments: this equation corresponds to strain from the initial state and at some point after plastic strain, unloadingloading occurs; elastic strain should be expressed by the following equation where τ * , * are the reference values of stresses and strains for a new cycle of elastic deformation. If we take ̇= 1̇ , then an account for the reference values is eliminated: for the next point in time, the previous values of stresses and strains are reference values. Consider an option of setting the shear stress at the boundary = 0 according to harmonic law τ 0 ( ) = τ max sin(ω ). The initial data are the same, and τ max =10 МPа, τ =1 МPа, 1 = 10 ⁄ . Figures 7 and 8 show the change in shear stress and shear strain over time in fixed sections r=2, 3, 4, and 6 m, corresponding to curves 1-4, at values of ω: ω = 5π (а) and ω =50π rad•s -1 (b).  Figure 8. Change in shear strain over time. As seen from figure 7, the amplitude of shear stresses in an elastoplastic medium is much (more than 3 times in the considered case) less than the amplitude of stresses in an elastic soil (curves 1 and 2, Figure 7). At ω = 5π and ω =10π rad•s -1 in the considered time interval, particles located further than 3 meters from the boundary are elastically strained, i.e. the plastic wave does not reach these points (curves 3 and 4 in Figure 7). There is also a decrease in the wave parameters (shear stress, shear strain, particle velocity, etc.) with distance from the initial section (Figures 7 and 8). Due to this wave attenuation, plastic strain of soil occurs in a bounded area. This is confirmed by the stress-strain diagrams shown in figure 9 for the considered options. The run of curves 4 (corresponding to sections r=6 m) completely coincides with the curves obtained for stresses and shear strains for elastic soil (the dependence of the shear stress on shear strain remains linear, Figure 9). Analysis of the graphs in figures 7 and 8 showed that some of the breaks in the curves in figures 7 and 8 at the initial points of time are associated with the break point in the stress-strain diagram during the transition to plastic strain or to unloading.

Discussion
From the analysis of the results obtained, it follows that the maximum shear stresses are reached near the application of the load. For all calculation options considered, the distribution of the maximum absolute values of shear stress, as well as the values of elastic limit along the radius, are shown in figure 10 (curves 1-3 correspond to 1 = 10 ⁄ , 50 ⁄ , 100 ⁄ ). Figure 10 clearly demonstrates the shear stress attenuation with distance from the point of applied load and the concentration of maximum values around the initial cylindrical section, where the external load is applied. Here it is also possible to determine whether the soil strain is plastic or elastic -prelimit strain: for the observed cases, the region of elastoplastic strain is no further than 3 m from the point of load application; outside this region, the tangential stresses do not exceed the elastic limits, that is, the elastic limit τ does not change from the original value. The reliability of the methods and programs were verified by comparing the numerical solution with the analytical solution for the case of the Prandtl scheme for the plastic strain of soil under constant tangential load (5) τ 0 ( ) = τ max = 40 KPa. Such a comparison of stresses at the wave front is shown in figure 11. Curve 1 corresponds to the analytical solution obtained in [6] by solving the problem for a shock shear wave (having a discontinuity at the wave front) by the method of characteristics, and curve 2 corresponds to a numerical solution using the developed technique. It should be noted that the comparison of the stress values at the wave front is the most unfavorable aspect for a numerical solution by the finite difference method (the stress discontinuity at the front is blurred in several cells); even in this case, it shows the suitability of the developed solution technique.  Figure 10. Distribution of maximum shear stresses and elastic limit along the radius. Figure 11. Comparison of the numerical result with the analytical solution [6]. As seen from figures 10 and 11, the attenuation intensity of a cylindrical shear wave with distance is approximately inversely proportional to the square root of the distance from the point of load application: attenuation occurs according to the law −γ , where γ ≈ 1/2.

Conclusions
Numerical results were obtained for a one-dimensional problem of cylindrical shear wave propagation in elastic and elastoplastic soil. A numerical method for solving the problem is developed using the finite difference method with a central difference scheme. The developed algorithm is implemented on a computer. From the results obtained, the attenuation of the parameters of cylindrical wave propagation (shear stress, shear strain, and velocity) with distance in an elastic and elastoplastic soil medium was revealed. The attenuation intensity of a cylindrical shear wave with distance is approximately inversely proportional to the square root of the radial coordinate. Attenuation of waves with distance is due to the fact that strain energy dissipates on the expanding cylindrical soil layer. In the case of a load exceeding the elastic limit of soil, the particles undergo plastic strains near the point of load application: in the considered options of the calculation, the boundary of elastic-plastic strain was 3 times the size of the cylindrical radius.