PHYSICAL REVIEW E 91, 022204 (2015)
Numerical investigation of the cylinder movement in granular matter Xue Zhang,* Daichao Sheng, George P. Kouretzis, Kristian Krabbenhoft, and Scott W. Sloan ARC Centre of Excellence for Geotechnical Science and Engineering, The University of Newcastle, Newcastle, Australia (Received 12 November 2014; revised manuscript received 6 January 2015; published 13 February 2015) We investigate numerically the mechanisms governing horizontal dragging of a rigid cylinder buried inside granular matter, with particular emphasis on enumerating drag and lift forces that resist cylinder movement. The recently proposed particle finite element method is employed, which combines the robustness of classical continuum mechanics formulations in terms of representing complex aspects of the material constitutive behavior, with the effectiveness of discrete element methods in simulating ultralarge deformation problems. The investigation focuses on the effect of embedment depth, cylinder roughness, granular matter macromechanical properties, and of the magnitude of the cylinder’s horizontal displacement on the amplitude of the resisting forces, which are discussed in light of published experimental data. Interpretation of the results provides insight on how the material flow around the cylinder affects the developing resistance, and a mechanism is proposed to describe the development of a steady-state drag force at large horizontal movements of the cylinder. DOI: 10.1103/PhysRevE.91.022204
PACS number(s): 45.70.Mg, 83.80.Fg, 46.15.−x
I. INTRODUCTION
The movement of an object inside a granular material is of relevance in many areas of engineering and applied science. Typical examples include soil-buried pipeline interaction problems in geotechnical engineering, stirring of granular media in pharmaceutical applications, and the locomotion of living organisms and robots buried in sand. Over the past decades, a number of experimental studies have been conducted to investigate the mechanics governing the movement of a cylinder fully immersed in granular matter. The experiments simulated dragging of a cylinder buried inside granular matter under gravitational stresses [1–6], or investigated the response of a fixed cylinder crossing the flow path of a sliding granular mass [7–10]. Despite the considerable number of studies devoted to this problem, several aspects related to the forces developing on the cylinder during relative movements still remain unclear, which is ascribable to certain limitations of the physical models used. One of particular interest is the influence of the embedment depth on the drag and lift forces developing during the movement of a cylinder on the horizontal plane. Trautmann and O’Rourke [11] and Ding et al. [12] simulated horizontal dragging of a cylinder in sand and glass beads, respectively, and measured drag and lift forces that were proportional to the embedment depth. On the contrary, Guillard et al. [9,10] reported that both drag and lift forces developing on the cylinder during experimental tests mimicking the stirring process of granular media eventually became depth independent, as the depth of embedment increased. Nevertheless Guillard et al. [9,10] refer mainly to the residual (steady-state) drag and lift forces, after “several rotations” of the cylinder, when arching effects could alter the gravitational stress field around the cylinder; the peak forces developing during the first half rotation of the cylinder increase with increasing embedment of the cylinder during stirring too. Another aspect of the problem not investigated deeply in the literature is how the macromechanical material properties,
*
[email protected]
1539-3755/2015/91(2)/022204(9)
associated with the granulometry as well as the volume density of the material (the ratio of the volume of solids over the total volume), contribute to and affect the forces developing to resist cylinder movement. Computational methods have a significant advantage in this respect, as they allow multiple “numerical experiments” to be performed to investigate different setups of the problem parametrically. Different methods were employed in the past to investigate this particular problem, with the bulk of the researchers adopting the discrete element method (DEM). Zhou et al. [13] used the DEM to simulate horizontal dragging of a cylinder through a confined granular bed, and found that the drag force is independent of the considered friction coefficient at the particle-cylinder interface, the roughness of the container, and the velocity of the cylinder (at relatively low speeds); it chiefly depends on the diameter of the cylinder, the volume density, the size of the particles, and the friction coefficient of the particles’ surface. Guillard et al. [10] also employed the DEM to study the mechanism triggering the development of lift forces, and attributed it to the gravitational pressure gradient resulting in nonsymmetric resistance to cylinder movement. Potiguar et al. [14] examined inertial effects on the resistance developing on a horizontally translating cylinder using the DEM too. They found that the drag force increases monotonically with the velocity and the diameter of the cylinder, unlike the lift force which appears to be less sensitive to the increase of the translation velocity. Potiguar et al. [14] attributed this feature to the decrease of granules being in contact with the cylinder at high speeds. Other research groups [2,15–17] working on pipeline engineering problems have followed the continuum modeling approach to simulate the behavior of pipelines buried in sand during permanent ground surface deformations, and used the classical displacement finite element method (FEM) and the finite element limit analysis method [16]. The above studies, which model the granular matter as a continuum and focus on the drag force only, do not discuss the development of the lift force on a pipeline. Apart from that, the range of application of classical continuum models is restricted to relatively small movements of the cylinder; for large displacements the granules will collapse into the cavity formed downstream of the cylinder, leading to severe mesh
022204-1
©2015 American Physical Society
ZHANG, SHENG, KOURETZIS, KRABBENHOFT, AND SLOAN
distortion and convergence problems, as discussed in detail by Yimsiri et al. [17]. Consequently, the ultimate drag force obtained from continuum models is the one at, or slightly after the “critical” displacement of the cylinder, defined herein as the horizontal displacement required for the complete failure surface to develop within the granular matter, when the reaction forces on the cylinder become independent of the elastic deformation properties of the granules. In reality, as the displacement of the cylinder increases beyond its critical value, the reaction forces may continue to increase due to the evolution of the failure surface, a phenomenon that cannot be addressed with classical continuum modeling formulations. In the following we refer to the state where the horizontal displacement has reached its critical value as the critical state, and the associated reaction forces as “critical forces.” In this study we will attempt to address certain limitations of the above-mentioned numerical studies, by employing the conic programming formulation proposed by Krabbenhoft et al. [18,19] and extended by Zhang et al. [20,21], under the framework of the recently developed continuum particle finite element method (PFEM) [20–24]. Compared to the DEM, the PFEM is a continuum approach that can make use of available continuum models, for which the calibration of the material parameters is less complicated, and the simulation is less time consuming. Besides that, another merit of this continuum approach is that simulations are not subjected to the typical limitations of the finite element method associated with mesh distortion and frictional contact issues. The basic idea [20,23] behind the PFEM is that mesh nodes are treated as a cloud of particles that can move freely, and even separate from the domain to which they originally belong. During a typical time step, the PFEM algorithm first identifies the material domain based on the position of the particles, and regenerates the mesh accordingly. Then a standard finite element procedure is employed to solve the governing equations on the new mesh. As a consequence, the PFEM combines the flexibility of mesh-free particle approaches in handling arbitrarily large deformation problems with the solid mathematical foundation of the traditional finite element method. This allows for the robust representation of complex aspects of the macromechanical behavior of the material, via implementation of appropriate constitutive laws. Following a short description of the problem and the adopted formulation, the presentation focuses on the drag and lift forces developing during the dragging of a fully buried cylinder horizontally through granular matter. We discuss the effect of burial depth; macromechanical properties of the granular material; roughness of the cylinder’s external surface; and of the magnitude of the horizontal displacements, on the amplitude of the resisting forces, and the flow mechanisms governing their development. Results from this exercise are compared against published data, and insights are discussed in view of existing experimental and numerical studies. II. PROBLEM DESCRIPTION
Figure 1 portrays the problem of a cylinder of diameter D which is fully buried in a granular bed with a free surface, at a depth H > 0.5D The cylinder is dragged along the horizontal plane by applying a horizontal displacement ux . Assuming
PHYSICAL REVIEW E 91, 022204 (2015)
Cylinder
D
H
ux
y x
FIG. 1. Schematic representation of horizontal dragging of a cylinder in granular matter.
the length of the cylinder is large compared to its diameter, the problem can be treated as plane strain when solving the governing equations that follow. The momentum conservation equations for the granular matter are written as ¨ in V , ∇ T σ + b = ρ u,
(1)
with the natural boundary conditions being NT σ = t, on S,
(2)
where σ are the stresses, b are body forces, t are tractions, ρ is the density of the granular matter, V is the domain under consideration with its boundary represented by S, and N corresponds to the component of the outward normal vector of the natural boundaries. A superposed dot denotes differentiation with respect to time, and ∇ is the linear operator relating strains ε and displacements u as ε = ∇u.
(3)
The cylinder is represented as a weightless rigid body, thus ignoring section deformation effects on the resistance forces. The mechanical behavior of the granular matter is represented by a rigid, perfectly-plastic constitutive model as F (σ ) 0, ˙ σ G(σ ), ε˙ = λ∇ ˙ (σ ) = 0, λF
(4)
λ˙ 0,
where F is the yield function, G is the plastic potential, and λ˙ is the plastic multiplier. For F = G, the flow rule is associated; otherwise it is nonassociated flow which is typical for granular materials. In this paper, both F and G obey the Mohr-Coulomb failure criterion: 2 + (σ F = (σxx − σyy )2 + 4σxy xx + σyy ) sin ϕ − 2c cos ϕ, 2 + (σ G = (σxx − σyy )2 + 4σxy (5) xx + σyy ) sin ψ, where φ is the friction angle, c is the cohesion, and ψ is the dilation angle, with φ = 0 implying that plastic flow takes place without any volume change. The assumed constitutive behavior ignores any elastic, reversible behavior of the granular material, and any arbitrary small displacement imposed on the cylinder will result in the cylinder-granular matter system reaching the critical state, defined earlier. This first-order approach simplifies the estimation of the critical drag and lift forces developing on the cylinder, which are of
022204-2
NUMERICAL INVESTIGATION OF THE CYLINDER . . .
PHYSICAL REVIEW E 91, 022204 (2015)
particular interest here, and its effectiveness has been proved in various similar studies. Albeit we utilize a simple rigid, perfectly-plastic model here, the consideration of elasticity, strain softening, and strain rate dependence is possible using the adopted framework. The contact between the cylinder’s external surface and its surrounding granular matter is simulated as frictional, obeying the classical Coulomb friction law [25]: (6)
where gN is the gap between the granular matter and the cylinder surface, p is the contact pressure which is positive corresponding to compression, q is the tangential stress, and μ is the friction coefficient. The first line of the constraints shown in Eq. (6) is known as the Hertz-Signorini-Moreau condition for frictionless contact, the meaning of which is obvious and can be expressed as follows: When the cylinder’s external surface is in contact with the granular matter one has the conditions that the gap gN = 0 and the normal force p 0, whereas the relations gN > 0 and p = 0 hold if there is a gap between the granular matter and the surface. The second line of (6) is a constraint on the tangential force which is governed by the standard Coulomb friction law [25]. Various numerical schemes can be used to solve the above governing equations. In the present work we treat them as a mathematical programming problem. In other words, the governing equations are discretized and then cast in terms of a standard optimization program which can be solved using a high-performance optimizer. The solution scheme utilized here was originally proposed by Krabbenhoft et al. [19] and Zhang et al. [20] and is documented briefly in the Appendix for the sake of completeness. III. CRITICAL-STATE DRAG AND LIFT FORCES
In this section we investigate the effect of the embedment depth, the macromechanical properties of the granular material, and the roughness of the cylinder on the criticalstate reaction forces due to an arbitrarily small horizontal displacement of the cylinder. Focusing first on the effect of shallow-to-moderate embedment depths (0.6D < H < 20D), we perform a set of benchmark analyses with fixed material properties to replicate experiments performed by Trautmann and O’Rourke [11]. These are outlined in Table I, and correspond to loose Cornell filter sand [11]. The diameter of the cylinder is D = 10.2 cm as per the tests of Trautmann and O’Rourke [11] and in the simulation both the width and the length of the granular bed are large enough
FIG. 2. Finite element meshes for the simulation.
to eliminate boundary effects. The cylinder-sand interface is assumed to be purely frictional with a friction coefficient μ = tan(0.5φ), in line with similar studies [16,17,26]. Figure 2 shows the mesh discretization of the computational domain. As can be seen, a fine mesh is utilized for the area close to the cylinder to ensure a suitable accuracy while a coarser one is used elsewhere to reduce the computational time of the modeling. In the simulation, a quadratic displacement–linear stress triangular element is utilized. In such an element, the displacement interpolation points are located at the corner and midside nodes, while the stress interpolation points are fixed at (λj−1 ,λj ,λj+1 ) = ( 16 , 46 , 16 ), j = 1, 2, 3, with λj being the area coordinates. Further details of the properties of this element can be found in Ref. [19]. Figure 3 presents the dimensionless critical-state drag force ND = FD /γ H DL (L is the length of the cylinder which is equal to L = 1 m in the plane strain model), and the dimensionless critical-state lift force NL = FL /γ H DL on the cylinder during horizontal dragging at different depths. The results indicate that the effect of the embedment depth on the developing forces is not monotonic, and two zones (A and B) of diverse behavior can be identified in Fig. 3. The criticalstate drag force remains practically constant as the cylinder is being dragged close to the surface bed (zone A, 0.6 < H/D < 2.2), whereas it increases linearly with H/D in zone B (2.2 < H/D < 20). Notably, due to the lack of experiments conducted for the case of the cylinder embedment located in -25 -20
20
ND (This study)
15
NL (This study)
-10
10
-5
5
0
0 0
14.8 kN/m3 31◦ 0◦ 0 kPa
25
B ND (Experimental data [11])
-15
TABLE I. Granular material properties considered in the benchmark analyses to investigate the effect of cylinder embedment on the developing drag and lift forces. Physical properties Unit weight, γ Critical friction angle, φ Dilation angle, ψ Cohesion, c
A
5
10
H/D
15
Critical state lift force, NL=FL/γHDL
|q| − μp 0
Critical state drag force, ND=FD/γHDL
gN 0, p 0, p gN = 0,
20
FIG. 3. (Color online) Variation of the dimensionless criticalstate drag force ND and the lift force NL with the embedment depth for the benchmark problem. The sign of the forces follows the coordinate system of Fig. 1. 022204-3
ZHANG, SHENG, KOURETZIS, KRABBENHOFT, AND SLOAN
(a)
(b)
Displacement increment (cm) 0.1 0
0.2
FIG. 4. (Color online) Critical-state flow patterns developing during horizontal translation of a cylinder buried at depths (a) H = 2D and (b) H = 12D.
zone A in Ref. [11], such a zone (which is identified in our studies) was not observed by Trautmann and O’Rourke [11]; nevertheless, the finding for zone B in our studies is compatible with those in Ref. [11], and a quantitative comparison of the computed and measured drag forces (Fig. 3) endorses the validity of our numerical approach. The dimensionless lift force, which was not measured by Trautmann and O’Rourke due to limitations of the experimental setup they developed, appears to peak when the cover of the cylinder crown is a minimum (H/D = 0.6), and decreases linearly in zone A. The dimensionless lift force developed remains constant throughout zone B, which implies that the absolute value of the lift force increases with the embedment depth. This behavior is in line with the lift measurements obtained by Ding et al. [12] while dragging a cylinder horizontally in granular matter at depths H < 6.5D. Guillard et al. [9,10], on the other hand, reported that drag and lift forces reached a plateau at a depth of about H = 10D, and remained constant as the embedment of the cylinder increased. At first reading this may seem contradictory but, as mentioned in the Introduction, Guillard et al. [9,10] experimentally modeled rotation (instead of dragging) of the cylinder on the horizontal plane. In addition they refer mainly to the steady-state rather than the critical-state forces, when stress redistribution due to arching is possible in the stirring problem. Figure 4 illustrates indicative flow patterns of granular matter, which define the failure surface around the cylinder, for different embedment depths. The transition from a surface failure mode [Fig. 4(a)] to a deep, localized failure mode [Fig. 4(b)] is compatible with the different zones identified
PHYSICAL REVIEW E 91, 022204 (2015)
in Fig. 3. Interestingly, the flow pattern in Fig. 4(b) is not symmetric, albeit the geometry of the problem is symmetric and the failure surface does not propagate up to the granular bed. However, one must bear in mind that the gradient of the gravitational stresses, which gives rise to the lift force during horizontal movement, results in nonsymmetric flow of the granular matter around the pipe. Indeed, if instead of a gravitational stress field we simulate cylinder movement under an isotropic stress field, the flow patter becomes symmetric with respect to the horizontal axis, and the lift force vanishes. Next we investigate parametrically the sensitivity of the critical-state reaction forces to the macromechanical properties of the granular matter (i.e., density, friction angle, dilation angle) as well as to the roughness of the cylinder’s external surface. In these simulations we keep only the diameter of the cylinder constant, and vary all the other parameters of the problem; the examined ranges are depicted in Fig. 5. A first, perhaps expected, observation from Fig. 5 is that the dimensionless reaction forces do not depend on the density of the granular matter, something that proves that the dimensionless terms are sufficient to describe the problem. It is worth mentioning here that the parametric investigation of the sensitivity of the reaction forces to the macrodensity of the material was performed while keeping the rest of the macromechanical properties constant. For an assembly of granular matter, the change of density is always associated with the change of the volume density (or void ratio), and thus of the shear strength of the material (for instance, the friction angle). Figure 5(b) illustrates the effect of the surface roughness of the cylinder on the induced drag and lift forces. It appears that the effect of the cylinder roughness is rather trivial, in line with the experimental findings of Trautmann and O’Rourke [11] where cylinders covered with coarse sandpaper, polyethylene plastic, or coated with gear oil all exhibited very similar behavior when dragged laterally in sand. The parametric investigations of Yimsiri et al. [17] (who used the finite element method) and Zhou et al. [13] (who used the discrete element method) which studied the effect of contact friction on the drag force, resulted in similar conclusions. On the contrary, it appears that both the critical-state drag and lift forces clearly depend on the macrofriction angle of the material [Fig. 5(c)]. The transition embedment depth defining the two distinct zones of behavior A and B first observed in Fig. 3 is not altered, and the dimensionless lift force remains constant at embedment depths of about H > 2.2D. Note that the drag force increasing with the macrofriction angle echoes the finding in the DEM simulations presented by Zhou et al. [13], who observed that the drag forces increase with the microfriction coefficient between the granules and thus the macrofriction angle (since the latter depends significantly on the former). Observe also in Fig. 5(c) that the rate of increase of the dimensionless drag force increases with increasing macrofriction angle. This is due to the change in the shape of the Mohr-Coulomb yield surface, and the corresponding increase in the resistance of the granular matter to shearing at higher confining stresses. Lastly, the drag force appears to be rather insensitive to the value of the dilation angle, which merely controls the
022204-4
15
-10
10
-5
5
0
0
Critical state drag force, ND=FD/γHDL
0 -30
5
10
H/D
15
ND (Φ = 25 ) ND (Φ = 30o) ND (Φ = 35o) NL (Φ = 25o) NL (Φ = 30o) NL (Φ = 35o) o
-25 -20 -15
30
(c)
25 20 15
-5
5
0 10
H/D
15
-15
20
0
25 20 15
-10
10
-5
5
0
0 0
10
5
-20
20
-10
0
-25
30
(b)
5
-30
10
15
H/D
ND (Ψ = 5 ) ND (Ψ = 10o) ND (Ψ = 15o) NL (Ψ = 5o) NL (Ψ = 10o) NL (Ψ = 15o) o
-25 -20 -15
20 30
(d)
25 20 15
-10
10
-5
5
0
0 0
5
10
15
H/D
Critical state lift force, NL=FL/γHDL
-15
20
ND (μ = tan(0.25Φ)) ND (μ = tan(0.5Φ)) ND (μ = tan(0.75Φ)) NL (μ = tan(0.25Φ)) NL (μ = tan(0.5Φ)) NL (μ = tan(0.75Φ))
Critical state lift force, NL=FL/γHDL
-20
25
-30
Critical state drag force, ND=FD/γHDL
-25
30
(a)
Critical state drag force, ND=FD/γHDL
ND (ρ = 1.0 g/cm3) ND (ρ = 2.0 g/cm3) ND (ρ = 3.0 g/cm3) NL (ρ = 1.0 g/cm3) NL (ρ = 2.0 g/cm3) NL (ρ = 3.0 g/cm3)
Critical state lift force, NL=FL/γHDL
-30
PHYSICAL REVIEW E 91, 022204 (2015)
Critical state lift force, NL=FL/γHDL
Critical state drag force, ND=FD/γHDL
NUMERICAL INVESTIGATION OF THE CYLINDER . . .
20
FIG. 5. (Color online) Variation of the dimensionless critical-state drag force ND and lift force NL with the embedment depth. (a) Effect of density ρ with φ = 30◦ , ψ = 10◦ , μ = tan(0.5φ); (b) effect of surface roughness μ with ρ = 2.0 g/cm3 , φ = 30◦ , ψ = 10◦ ; (c) effect of friction angle φ with ρ = 2.0 g/cm3 , ψ = 10◦ , μ = tan(0.5φ); (d) effect of dilation angle ψ with ρ = 2.0 g/cm3 , φ = 30◦ , μ = tan(0.5φ).
-15
IV. DRAG AND LIFT FORCES DEVELOPING UNDER LARGE HORIZONTAL DISPLACEMENTS
15 ND NL
-10
In the previous section we investigated critical-state drag and lift forces at the instant where a complete failure surface is formed within the granular matter. Now we increase gradually the applied drag displacement on the cylinder ux from about 0.02D to 15D, to explore the evolution of the drag and lift forces at very large deformations by means of the PFEM. Initially we consider the cylinder being dragged at various embedment depths in loose Cornell filter sand (Table I), with a zero dilation angle, i.e., deforming at constant volume upon shearing. Dimensionless force-displacement curves for the analysis for embedment depth H = 4D are plotted in Fig. 6. Interestingly, the dimensionless drag force increases from the
10
Steady state
Critical State
-5
5
0
0 0
5
ux/D
10
Critical state lift force, NL=FL/γHDL
critical-state value ND = 7.2 (see also Fig. 3) to ND = 8, when a steady state is achieved at approximately ux = 10D. In contrast, the dimensionless lift force NL remains constant. To shed some light on this, we plot in Fig. 7 the evolution of displacement increment contours and of the deformed material free surface, with increasing drag displacement. This shows
Critical state drag force, ND=FD/γHDL
magnitude of plastic volume expansion upon shearing. The dilation angle does not greatly influence the shear resistance of the granular material in unrestrained problems, when volume expansion is not suppressed by the boundary conditions. We ensure this here by controlling the applied horizontal displacement, and therefore the results are somewhat expected. This is not the case though for the lift force, which increases with increasing dilation angle and associated volume expansion of the granular matter within the failure surface, which drifts together with the cylinder towards the free surface.
15
FIG. 6. (Color online) Variation of the dimensionless drag force ND and lift force NL with the applied displacement. Results correspond to cylinder embedment H = 4D and material dilation angle ψ = 0◦ .
022204-5
ZHANG, SHENG, KOURETZIS, KRABBENHOFT, AND SLOAN
(a)
PHYSICAL REVIEW E 91, 022204 (2015)
ux = 0.1D
(a) Displacement increment (cm) 0.1
Displacement Increment (cm) 0.1
ux = 0.1D
0
0
0.2
0.2
(b)
ux = 0.8D
(b)
(a)
ux = 0.8D
ux = 10D
(c)
ux = 15D
(d)
(c)
ux = 10D (d) FIG. 8. (Color online) Evolution of displacement increment contours illustrating flow of granules within the failure surface with increasing drag displacement ux : (a) ux = 0.1D, (b) ux = 1D, (c) ux = 10D, (d) ux = 15D. Results correspond to cylinder embedment H = 8D and material dilation angle ψ = 8◦ .
ux = 15D FIG. 7. (Color online) Evolution of displacement increment contours illustrating flow of granules within the failure surface with increasing drag displacement ux : (a) ux = 0.1D, (b) ux = 1D, (c) ux = 10D, (d) ux = 15D. Results correspond to cylinder embedment H = 4D and material dilation angle ψ = 0◦ .
that the surface of the granular bed deforms to accommodate the movement of the cylinder, and a heap is formed as granular material is accumulated upstream of the cylinder. Formation of the heap increases the embedment depth in the vicinity of the cylinder [Figs. 7(a)–7(c)] and, therefore, the dimensionless drag force, which is proportional to the depth (see also Fig. 3). At large horizontal displacements (ux > 10D) a steady state is reached, as the volume of the heap is limited by the shear strength of the material [Figs. 7(c) and 7(d)]. Unlike the drag force, the dimensionless lift force NL does not depend on the embedment depth, as discussed earlier (see also Fig. 3), so the steady-state lift force is equal to the critical-state lift force. At larger embedment depths, where a local failure mode will prevail and the failure surface does not propagate to the surface of the granular bed, the steady-state drag force will also be equal to the corresponding critical-state value. The above discussion focused on a material that does not exhibit volume expansion during dragging (ψ = 0◦ ). It is worth investigating the effect of a small amount of dilation (ψ < 8◦ ) on the drag and lift forces, albeit considering con-
stant material dilation with the rigid, plastic failure criterion employed suggests that the granular matter will continue to dilate indefinitely, as long as shearing continues. This type of behavior is somewhat fictitious; in reality most granular materials will reach a state where they will deform at constant volume upon further shearing (ψ = 0◦ ). However, insights gained from this analysis contribute to the understanding of the mechanisms governing cylinder-granular matter interaction at large horizontal displacements. Figure 8 is similar to Fig. 7, but resulted from an analysis for dilation angle ψ = 8◦ . Volume expansion upon continuous shearing results in the formation of an unrealistic heap. This is also depicted in Fig. 9, presenting the variation of the dimensionless drag and lift forces with increasing horizontal displacement of the cylinder. Drag forces increase following the increase of the embedment of the cylinder due to the accumulation of the expanding granular matter. The rate of this increase is a function of the value of the dilation angle; higher dilation angles will result in larger heap volumes, and thus higher increase rates. This is analogous to the linear increase of the deviatoric stress with continued undrained shearing in triaxial test simulations with the Mohr-Coulomb model and a dilation angle ψ > 0◦ . On the other hand, dimensionless lift forces remain constant and equal to the critical-state value, which suggests that the mechanism for the development of the lift forces is related to the gradient of the gravitational
022204-6
NUMERICAL INVESTIGATION OF THE CYLINDER . . .
20
ψ=0 ψ=4 ψ=8
-15
15
ND
-10
10
NL
-5
5
0
0 0
5
ux/D
10
APPENDIX: SOLUTION ALGORITHM
Critical state lift force, NL=FL/γHDL
Critical state drag force, ND=FD/γHDL
-20
PHYSICAL REVIEW E 91, 022204 (2015)
1. Time discretization
The momentum conservation equations (1) can be discretized using the θ method as vn+1 − vn , (A1)
t un+1 − un , (A2) θ2 vn+1 + (1 − θ2 )vn =
t where v are velocities, subscripts n and n + 1 refer to the known and new, unknown states, and t = tn+1 − tn is the time step. Rearranging the above equations leads to ∇ T [θ1 σ n+1 + (1 − θ1 )σ n ] + b = ρ
15
u ∇ T σ n+1 + b˜ = ρ˜ 2 ,
t 1 u − (1 − θ2 )vn , vn+1 = θ2 t
FIG. 9. (Color online) Variation of the dimensionless drag force ND and lift force NL developing on the cylinder which moves in granular matter with different dilation angles. Results correspond to cylinder embedment H = 8D.
(A3) (A4)
where u = un+1 − un and ρ , θ1 θ2 1 − θ1 T vn 1 + ∇ σ n, b˜ = b + ρ˜ θ1
t θ1
stresses, and drifting of the cylinder towards the free surface due to volume expansion.
ρ˜ =
V. CONCLUSIONS
Use of the recently developed particle finite element method (PFEM) allowed us to study the mechanisms governing horizontal movement of a rigid cylinder inside a bed of granular matter. This method overcomes the technical restrictions of simple physical models and the shortcomings of classical numerical approaches. Some interesting insights emerging from this study include the numerical documentation of the development of a lift force on the cylinder. After a transition embedment depth of about H/D = 2.2 (which is related to the shape of the failure surface but not to the mechanical properties of the granular matter), the dimensionless lift force becomes independent of the embedment depth and the magnitude of the applied cylinder horizontal drag, and its value is only a function of the shear strength and volume expansion potential of the granular matter. On the other hand, the dimensionless drag force increases linearly with the embedment of the cylinder after the transition embedment depth H/D = 2.2, in line with published experimental results, and its critical-state value is a function of the macrofriction angle of the material. As the applied cylinder horizontal drag increases beyond the critical displacement, which is the horizontal displacement required for the complete failure surface to develop within the granular matter, the dimensionless drag force rises up to a steady-state value. This phenomenon is attributed to the accumulation of granular material above the cylinder, which might take place even when the granular material does not exhibit a tendency for volume expansion upon continuous shearing. The particle finite element method (PFEM) can be further adapted to incorporate more complex aspects of granular material behavior such as a material’s volume expansion or contraction being dependent on the applied shearing, or even extended to three dimensions, to treat the nonsymmetric problem of stirring of granular matter.
(A5) (A6)
The natural boundary conditions (2) are approximated in an analogous manner leading to NT σ n+1 = ˜t, on S,
(A7)
where ˜t =
1 1 − θ1 T t− N σ n. θ1 θ1
(A8)
Following [20], the above problem can be stated in terms of a min-max problem: min max
u (σ ,r)n+1
subject to
˜ ˜ σ n+1 ,∇( u)V − b, u V − t, uS − 12 t 2 rn+1 ,ρ˜ −1 rn+1 V + rn+1 , uV
(A9)
F (σ n+1 ) 0
on the basis of the Hellinger-Reissner variational principle. In the above, rn+1 are a set of variables interpreted as dynamic forces [20], and the notation x,yA = xT ydA (A10) A
is utilized. The equivalence between the optimization problem (A9) and the governing equations at hand is established by demonstrating that the Euler-Lagrange equations associated with (A9) indeed reproduce the governing equations [20]. 2. Spatial discretization
Using the standard finite element notation, the following approximations,
022204-7
σ (x) ≈ Nσ σˆ ,
(A11)
r(x) ≈ Nr rˆ ,
(A12)
ˆ u(x) ≈ Nu u,
(A13)
ˆ n∇u ≈ Bu u,
(A14)
ZHANG, SHENG, KOURETZIS, KRABBENHOFT, AND SLOAN
PHYSICAL REVIEW E 91, 022204 (2015)
are introduced for the state variables, where σˆ , rˆ , and uˆ are the nodal variables, N matrices contain the shape functions, and Bu = ∇Nu . Substituting the above approximations into the variational principle (A9) results in the following discrete principle:
At this stage, the contact constraints (6) are imposed on all potential contact nodes (i.e., mesh nodes located on the boundaries), which leads to a final problem of the type
min max
σˆ Tn+1 B uˆ
subject to
1 − t 2 rˆ Tn+1 Dˆrn+1 + rˆ Tn+1 A uˆ 2 j F σˆ n+1 0, j = 1, . . . ,nσ
uˆ (σˆ ,ˆr)n+1
− f uˆ
V
V
subject to (A15)
V
S
NTu ˜tdS,
NTr ρ˜ −1 Nσ dV ,
(A17) (A18)
A= V
NTu Nr dV .
(A19)
Then, solving the minimization part of (A15) gives a maximization problem as maximize (σˆ ,ˆr)n+1
subject to
1 − t 2 rˆ Tn+1 Dˆrn+1 2 BT σˆ n+1 + AT rˆ n+1 = f j F σˆ n+1 0, j = 1, . . . ,nσ .
nc 1 − t 2 rˆ Tn+1 Dˆrn+1 − g0j pj 2 j =1 BT σˆ n+1 + AT rˆ n+1 + ET ρ = f F σˆ in+1 0, i = 1, . . . ,nσ
pj = −nT ρj ,
˜ + NTu bdV
D=
(σˆ ,ˆr,pj )n+1
T
where nσ is the number of Gauss integration points, and NTσ Bu dV , (A16) B= f=
maximize
(A20)
[1] R. Candelier and O. Dauchot, Phys. Rev. Lett. 103, 128001 (2009). [2] N. Daiyan, S. Kenny, R. Phillips, and R. Popescu, Can. Geotech. J. 48, 1683 (2011). [3] T. Hsu, Y. Chen, and W. Hung, J. Transp. Eng. 132, 175 (2006). [4] T.-W. Hsu, Can. Geotech. J. 33, 180 (1996). [5] K. Wieghardt, Annu. Rev. Fluid Mech. 7, 89 (1975). [6] A. Seguin, Y. Bertho, F. Martinez, J. Crassous, and P. Gondret, Phys. Rev. E 87, 012201 (2013). [7] Y. Amarouchene, J. F. Boudet, and H. Kellay, Phys. Rev. Lett. 86, 4286 (2001). [8] D. Chehata, R. Zenit, and C. R. Wassgren, Phys. Fluids (1994present) 15, 1622 (2003). [9] F. Guillard, Y. Forterre, and O. Pouliquen, Phys. Rev. Lett. 110, 138303 (2013). [10] F. Guillard, Y. Forterre, and O. Pouliquen, Phys. Fluids 26, 043301 (2014). [11] C. Trautmann and T. O’Rourke, J. Geotech. Eng. 111, 1077 (1985). [12] Y. Ding, N. Gravish, and D. I. Goldman, Phys. Rev. Lett. 106, 028001 (2011).
(A21)
j = 1, . . . ,nc
qj = −nˆ T ρj |qj | − μpj 0, where ρ = (ρ1 ,ρ2 )T are the nodal forces, n = (n1 ,n2 )T and nˆ = (−n2 ,n1 )T are the normal and the tangential of the rigid boundary, E is an index matrix of zeros and ones, and nc is the number of potential contacts. The above problem can be transformed into a standard form of the second-order cone program (SOCP) and then solved using the high-performance optimization solver MOSEK [27] with a default setup for the convergence tolerance. The transformation of (A21) into SOCP standard form is straightforward and has been documented in Refs. [18,19]. In the course of solving the problem, the kinematic variables (displacement increments and plastic multipliers) are recovered as the dual variables, or Lagrange multipliers, associated with the discrete equilibrium constraints. To date, the above scheme has been successfully applied to a variety of challenging problems including limit analysis [28], elastoplastic static analysis [18,19], rigid, plastic dynamic analysis with large deformations [20,21,24], and granular contact dynamic analysis [29].
[13] F. Zhou, S. G. Advani, and E. D. Wetzel, Phys. Fluids 19, 013301 (2007). [14] F. Q. Potiguar and Y. Ding, Phys. Rev. E 88, 012204 (2013). [15] J. Jung, T. O’Rourke, and N. Olson, J. Geotech. Geoenviron. Eng. 139, 2028 (2013). [16] G. P. Kouretzis, D. Sheng, and S. W. Sloan, Comput. Geotech. 54, 53 (2013). [17] S. Yimsiri, K. Soga, K. Yoshizaki, G. Dasari, and T. O’Rourke, J. Geotech. Geoenviron. Eng. 130, 830 (2004). [18] K. Krabbenhoft and A. V. Lyamin, Comput. Methods Appl. Mech. Eng. 209-212, 239 (2012). [19] K. Krabbenhøft, A. V. Lyamin, and S. W. Sloan, Int. J. Solids Struct. 44, 1533 (2007). [20] X. Zhang, K. Krabbenhoft, D. M. Pedroso, A. V. Lyamin, D. Sheng, M. V. da Silva, and D. Wang, Comput. Geotech. 54, 133 (2013). [21] X. Zhang, K. Krabbenhoft, and D. Sheng, Granular Matter 16, 609 (2014). [22] S. R. Idelsohn, E. O˜nate, and F. D. Pin, Int. J. Numer. Methods Eng. 61, 964 (2004). [23] E. O˜nate, S. R. Idelsohn, F. Del Pin, and R. Aubry, Int. J. Comput. Methods 01, 267 (2004).
022204-8
NUMERICAL INVESTIGATION OF THE CYLINDER . . .
PHYSICAL REVIEW E 91, 022204 (2015)
[24] X. Zhang, K. Krabbenhoft, D. Sheng, and W. Li, Comput. Mech. 55, 167 (2015). [25] P. Wriggers, Computational Contact Mechanics (Wiley, Chichester, UK, 2002). [26] G. P. Kouretzis, K. Krabbenhøft, D. Sheng, and S. W. Sloan, Can. Geotech. J. 51, 1087 (2014).
[27] E. D. Andersen, C. Roos, and T. Terlaky, Math. Program. 95, 249 (2003). [28] K. Krabbenhoft and L. Damkilde, Int. J. Numer. Methods Eng. 56, 165 (2003). [29] K. Krabbenhoft, J. Huang, M. V. da Silva, and A. V. Lyamin, Granular Matter 14, 607 (2012).
022204-9
Numerical investigation of the cylinder movement in granular matter Xue Zhang,* Daichao Sheng, George P. Kouretzis, Kristian Krabbenhoft, and Scott W. Sloan ARC Centre of Excellence for Geotechnical Science and Engineering, The University of Newcastle, Newcastle, Australia (Received 12 November 2014; revised manuscript received 6 January 2015; published 13 February 2015) We investigate numerically the mechanisms governing horizontal dragging of a rigid cylinder buried inside granular matter, with particular emphasis on enumerating drag and lift forces that resist cylinder movement. The recently proposed particle finite element method is employed, which combines the robustness of classical continuum mechanics formulations in terms of representing complex aspects of the material constitutive behavior, with the effectiveness of discrete element methods in simulating ultralarge deformation problems. The investigation focuses on the effect of embedment depth, cylinder roughness, granular matter macromechanical properties, and of the magnitude of the cylinder’s horizontal displacement on the amplitude of the resisting forces, which are discussed in light of published experimental data. Interpretation of the results provides insight on how the material flow around the cylinder affects the developing resistance, and a mechanism is proposed to describe the development of a steady-state drag force at large horizontal movements of the cylinder. DOI: 10.1103/PhysRevE.91.022204
PACS number(s): 45.70.Mg, 83.80.Fg, 46.15.−x
I. INTRODUCTION
The movement of an object inside a granular material is of relevance in many areas of engineering and applied science. Typical examples include soil-buried pipeline interaction problems in geotechnical engineering, stirring of granular media in pharmaceutical applications, and the locomotion of living organisms and robots buried in sand. Over the past decades, a number of experimental studies have been conducted to investigate the mechanics governing the movement of a cylinder fully immersed in granular matter. The experiments simulated dragging of a cylinder buried inside granular matter under gravitational stresses [1–6], or investigated the response of a fixed cylinder crossing the flow path of a sliding granular mass [7–10]. Despite the considerable number of studies devoted to this problem, several aspects related to the forces developing on the cylinder during relative movements still remain unclear, which is ascribable to certain limitations of the physical models used. One of particular interest is the influence of the embedment depth on the drag and lift forces developing during the movement of a cylinder on the horizontal plane. Trautmann and O’Rourke [11] and Ding et al. [12] simulated horizontal dragging of a cylinder in sand and glass beads, respectively, and measured drag and lift forces that were proportional to the embedment depth. On the contrary, Guillard et al. [9,10] reported that both drag and lift forces developing on the cylinder during experimental tests mimicking the stirring process of granular media eventually became depth independent, as the depth of embedment increased. Nevertheless Guillard et al. [9,10] refer mainly to the residual (steady-state) drag and lift forces, after “several rotations” of the cylinder, when arching effects could alter the gravitational stress field around the cylinder; the peak forces developing during the first half rotation of the cylinder increase with increasing embedment of the cylinder during stirring too. Another aspect of the problem not investigated deeply in the literature is how the macromechanical material properties,
*
[email protected]
1539-3755/2015/91(2)/022204(9)
associated with the granulometry as well as the volume density of the material (the ratio of the volume of solids over the total volume), contribute to and affect the forces developing to resist cylinder movement. Computational methods have a significant advantage in this respect, as they allow multiple “numerical experiments” to be performed to investigate different setups of the problem parametrically. Different methods were employed in the past to investigate this particular problem, with the bulk of the researchers adopting the discrete element method (DEM). Zhou et al. [13] used the DEM to simulate horizontal dragging of a cylinder through a confined granular bed, and found that the drag force is independent of the considered friction coefficient at the particle-cylinder interface, the roughness of the container, and the velocity of the cylinder (at relatively low speeds); it chiefly depends on the diameter of the cylinder, the volume density, the size of the particles, and the friction coefficient of the particles’ surface. Guillard et al. [10] also employed the DEM to study the mechanism triggering the development of lift forces, and attributed it to the gravitational pressure gradient resulting in nonsymmetric resistance to cylinder movement. Potiguar et al. [14] examined inertial effects on the resistance developing on a horizontally translating cylinder using the DEM too. They found that the drag force increases monotonically with the velocity and the diameter of the cylinder, unlike the lift force which appears to be less sensitive to the increase of the translation velocity. Potiguar et al. [14] attributed this feature to the decrease of granules being in contact with the cylinder at high speeds. Other research groups [2,15–17] working on pipeline engineering problems have followed the continuum modeling approach to simulate the behavior of pipelines buried in sand during permanent ground surface deformations, and used the classical displacement finite element method (FEM) and the finite element limit analysis method [16]. The above studies, which model the granular matter as a continuum and focus on the drag force only, do not discuss the development of the lift force on a pipeline. Apart from that, the range of application of classical continuum models is restricted to relatively small movements of the cylinder; for large displacements the granules will collapse into the cavity formed downstream of the cylinder, leading to severe mesh
022204-1
©2015 American Physical Society
ZHANG, SHENG, KOURETZIS, KRABBENHOFT, AND SLOAN
distortion and convergence problems, as discussed in detail by Yimsiri et al. [17]. Consequently, the ultimate drag force obtained from continuum models is the one at, or slightly after the “critical” displacement of the cylinder, defined herein as the horizontal displacement required for the complete failure surface to develop within the granular matter, when the reaction forces on the cylinder become independent of the elastic deformation properties of the granules. In reality, as the displacement of the cylinder increases beyond its critical value, the reaction forces may continue to increase due to the evolution of the failure surface, a phenomenon that cannot be addressed with classical continuum modeling formulations. In the following we refer to the state where the horizontal displacement has reached its critical value as the critical state, and the associated reaction forces as “critical forces.” In this study we will attempt to address certain limitations of the above-mentioned numerical studies, by employing the conic programming formulation proposed by Krabbenhoft et al. [18,19] and extended by Zhang et al. [20,21], under the framework of the recently developed continuum particle finite element method (PFEM) [20–24]. Compared to the DEM, the PFEM is a continuum approach that can make use of available continuum models, for which the calibration of the material parameters is less complicated, and the simulation is less time consuming. Besides that, another merit of this continuum approach is that simulations are not subjected to the typical limitations of the finite element method associated with mesh distortion and frictional contact issues. The basic idea [20,23] behind the PFEM is that mesh nodes are treated as a cloud of particles that can move freely, and even separate from the domain to which they originally belong. During a typical time step, the PFEM algorithm first identifies the material domain based on the position of the particles, and regenerates the mesh accordingly. Then a standard finite element procedure is employed to solve the governing equations on the new mesh. As a consequence, the PFEM combines the flexibility of mesh-free particle approaches in handling arbitrarily large deformation problems with the solid mathematical foundation of the traditional finite element method. This allows for the robust representation of complex aspects of the macromechanical behavior of the material, via implementation of appropriate constitutive laws. Following a short description of the problem and the adopted formulation, the presentation focuses on the drag and lift forces developing during the dragging of a fully buried cylinder horizontally through granular matter. We discuss the effect of burial depth; macromechanical properties of the granular material; roughness of the cylinder’s external surface; and of the magnitude of the horizontal displacements, on the amplitude of the resisting forces, and the flow mechanisms governing their development. Results from this exercise are compared against published data, and insights are discussed in view of existing experimental and numerical studies. II. PROBLEM DESCRIPTION
Figure 1 portrays the problem of a cylinder of diameter D which is fully buried in a granular bed with a free surface, at a depth H > 0.5D The cylinder is dragged along the horizontal plane by applying a horizontal displacement ux . Assuming
PHYSICAL REVIEW E 91, 022204 (2015)
Cylinder
D
H
ux
y x
FIG. 1. Schematic representation of horizontal dragging of a cylinder in granular matter.
the length of the cylinder is large compared to its diameter, the problem can be treated as plane strain when solving the governing equations that follow. The momentum conservation equations for the granular matter are written as ¨ in V , ∇ T σ + b = ρ u,
(1)
with the natural boundary conditions being NT σ = t, on S,
(2)
where σ are the stresses, b are body forces, t are tractions, ρ is the density of the granular matter, V is the domain under consideration with its boundary represented by S, and N corresponds to the component of the outward normal vector of the natural boundaries. A superposed dot denotes differentiation with respect to time, and ∇ is the linear operator relating strains ε and displacements u as ε = ∇u.
(3)
The cylinder is represented as a weightless rigid body, thus ignoring section deformation effects on the resistance forces. The mechanical behavior of the granular matter is represented by a rigid, perfectly-plastic constitutive model as F (σ ) 0, ˙ σ G(σ ), ε˙ = λ∇ ˙ (σ ) = 0, λF
(4)
λ˙ 0,
where F is the yield function, G is the plastic potential, and λ˙ is the plastic multiplier. For F = G, the flow rule is associated; otherwise it is nonassociated flow which is typical for granular materials. In this paper, both F and G obey the Mohr-Coulomb failure criterion: 2 + (σ F = (σxx − σyy )2 + 4σxy xx + σyy ) sin ϕ − 2c cos ϕ, 2 + (σ G = (σxx − σyy )2 + 4σxy (5) xx + σyy ) sin ψ, where φ is the friction angle, c is the cohesion, and ψ is the dilation angle, with φ = 0 implying that plastic flow takes place without any volume change. The assumed constitutive behavior ignores any elastic, reversible behavior of the granular material, and any arbitrary small displacement imposed on the cylinder will result in the cylinder-granular matter system reaching the critical state, defined earlier. This first-order approach simplifies the estimation of the critical drag and lift forces developing on the cylinder, which are of
022204-2
NUMERICAL INVESTIGATION OF THE CYLINDER . . .
PHYSICAL REVIEW E 91, 022204 (2015)
particular interest here, and its effectiveness has been proved in various similar studies. Albeit we utilize a simple rigid, perfectly-plastic model here, the consideration of elasticity, strain softening, and strain rate dependence is possible using the adopted framework. The contact between the cylinder’s external surface and its surrounding granular matter is simulated as frictional, obeying the classical Coulomb friction law [25]: (6)
where gN is the gap between the granular matter and the cylinder surface, p is the contact pressure which is positive corresponding to compression, q is the tangential stress, and μ is the friction coefficient. The first line of the constraints shown in Eq. (6) is known as the Hertz-Signorini-Moreau condition for frictionless contact, the meaning of which is obvious and can be expressed as follows: When the cylinder’s external surface is in contact with the granular matter one has the conditions that the gap gN = 0 and the normal force p 0, whereas the relations gN > 0 and p = 0 hold if there is a gap between the granular matter and the surface. The second line of (6) is a constraint on the tangential force which is governed by the standard Coulomb friction law [25]. Various numerical schemes can be used to solve the above governing equations. In the present work we treat them as a mathematical programming problem. In other words, the governing equations are discretized and then cast in terms of a standard optimization program which can be solved using a high-performance optimizer. The solution scheme utilized here was originally proposed by Krabbenhoft et al. [19] and Zhang et al. [20] and is documented briefly in the Appendix for the sake of completeness. III. CRITICAL-STATE DRAG AND LIFT FORCES
In this section we investigate the effect of the embedment depth, the macromechanical properties of the granular material, and the roughness of the cylinder on the criticalstate reaction forces due to an arbitrarily small horizontal displacement of the cylinder. Focusing first on the effect of shallow-to-moderate embedment depths (0.6D < H < 20D), we perform a set of benchmark analyses with fixed material properties to replicate experiments performed by Trautmann and O’Rourke [11]. These are outlined in Table I, and correspond to loose Cornell filter sand [11]. The diameter of the cylinder is D = 10.2 cm as per the tests of Trautmann and O’Rourke [11] and in the simulation both the width and the length of the granular bed are large enough
FIG. 2. Finite element meshes for the simulation.
to eliminate boundary effects. The cylinder-sand interface is assumed to be purely frictional with a friction coefficient μ = tan(0.5φ), in line with similar studies [16,17,26]. Figure 2 shows the mesh discretization of the computational domain. As can be seen, a fine mesh is utilized for the area close to the cylinder to ensure a suitable accuracy while a coarser one is used elsewhere to reduce the computational time of the modeling. In the simulation, a quadratic displacement–linear stress triangular element is utilized. In such an element, the displacement interpolation points are located at the corner and midside nodes, while the stress interpolation points are fixed at (λj−1 ,λj ,λj+1 ) = ( 16 , 46 , 16 ), j = 1, 2, 3, with λj being the area coordinates. Further details of the properties of this element can be found in Ref. [19]. Figure 3 presents the dimensionless critical-state drag force ND = FD /γ H DL (L is the length of the cylinder which is equal to L = 1 m in the plane strain model), and the dimensionless critical-state lift force NL = FL /γ H DL on the cylinder during horizontal dragging at different depths. The results indicate that the effect of the embedment depth on the developing forces is not monotonic, and two zones (A and B) of diverse behavior can be identified in Fig. 3. The criticalstate drag force remains practically constant as the cylinder is being dragged close to the surface bed (zone A, 0.6 < H/D < 2.2), whereas it increases linearly with H/D in zone B (2.2 < H/D < 20). Notably, due to the lack of experiments conducted for the case of the cylinder embedment located in -25 -20
20
ND (This study)
15
NL (This study)
-10
10
-5
5
0
0 0
14.8 kN/m3 31◦ 0◦ 0 kPa
25
B ND (Experimental data [11])
-15
TABLE I. Granular material properties considered in the benchmark analyses to investigate the effect of cylinder embedment on the developing drag and lift forces. Physical properties Unit weight, γ Critical friction angle, φ Dilation angle, ψ Cohesion, c
A
5
10
H/D
15
Critical state lift force, NL=FL/γHDL
|q| − μp 0
Critical state drag force, ND=FD/γHDL
gN 0, p 0, p gN = 0,
20
FIG. 3. (Color online) Variation of the dimensionless criticalstate drag force ND and the lift force NL with the embedment depth for the benchmark problem. The sign of the forces follows the coordinate system of Fig. 1. 022204-3
ZHANG, SHENG, KOURETZIS, KRABBENHOFT, AND SLOAN
(a)
(b)
Displacement increment (cm) 0.1 0
0.2
FIG. 4. (Color online) Critical-state flow patterns developing during horizontal translation of a cylinder buried at depths (a) H = 2D and (b) H = 12D.
zone A in Ref. [11], such a zone (which is identified in our studies) was not observed by Trautmann and O’Rourke [11]; nevertheless, the finding for zone B in our studies is compatible with those in Ref. [11], and a quantitative comparison of the computed and measured drag forces (Fig. 3) endorses the validity of our numerical approach. The dimensionless lift force, which was not measured by Trautmann and O’Rourke due to limitations of the experimental setup they developed, appears to peak when the cover of the cylinder crown is a minimum (H/D = 0.6), and decreases linearly in zone A. The dimensionless lift force developed remains constant throughout zone B, which implies that the absolute value of the lift force increases with the embedment depth. This behavior is in line with the lift measurements obtained by Ding et al. [12] while dragging a cylinder horizontally in granular matter at depths H < 6.5D. Guillard et al. [9,10], on the other hand, reported that drag and lift forces reached a plateau at a depth of about H = 10D, and remained constant as the embedment of the cylinder increased. At first reading this may seem contradictory but, as mentioned in the Introduction, Guillard et al. [9,10] experimentally modeled rotation (instead of dragging) of the cylinder on the horizontal plane. In addition they refer mainly to the steady-state rather than the critical-state forces, when stress redistribution due to arching is possible in the stirring problem. Figure 4 illustrates indicative flow patterns of granular matter, which define the failure surface around the cylinder, for different embedment depths. The transition from a surface failure mode [Fig. 4(a)] to a deep, localized failure mode [Fig. 4(b)] is compatible with the different zones identified
PHYSICAL REVIEW E 91, 022204 (2015)
in Fig. 3. Interestingly, the flow pattern in Fig. 4(b) is not symmetric, albeit the geometry of the problem is symmetric and the failure surface does not propagate up to the granular bed. However, one must bear in mind that the gradient of the gravitational stresses, which gives rise to the lift force during horizontal movement, results in nonsymmetric flow of the granular matter around the pipe. Indeed, if instead of a gravitational stress field we simulate cylinder movement under an isotropic stress field, the flow patter becomes symmetric with respect to the horizontal axis, and the lift force vanishes. Next we investigate parametrically the sensitivity of the critical-state reaction forces to the macromechanical properties of the granular matter (i.e., density, friction angle, dilation angle) as well as to the roughness of the cylinder’s external surface. In these simulations we keep only the diameter of the cylinder constant, and vary all the other parameters of the problem; the examined ranges are depicted in Fig. 5. A first, perhaps expected, observation from Fig. 5 is that the dimensionless reaction forces do not depend on the density of the granular matter, something that proves that the dimensionless terms are sufficient to describe the problem. It is worth mentioning here that the parametric investigation of the sensitivity of the reaction forces to the macrodensity of the material was performed while keeping the rest of the macromechanical properties constant. For an assembly of granular matter, the change of density is always associated with the change of the volume density (or void ratio), and thus of the shear strength of the material (for instance, the friction angle). Figure 5(b) illustrates the effect of the surface roughness of the cylinder on the induced drag and lift forces. It appears that the effect of the cylinder roughness is rather trivial, in line with the experimental findings of Trautmann and O’Rourke [11] where cylinders covered with coarse sandpaper, polyethylene plastic, or coated with gear oil all exhibited very similar behavior when dragged laterally in sand. The parametric investigations of Yimsiri et al. [17] (who used the finite element method) and Zhou et al. [13] (who used the discrete element method) which studied the effect of contact friction on the drag force, resulted in similar conclusions. On the contrary, it appears that both the critical-state drag and lift forces clearly depend on the macrofriction angle of the material [Fig. 5(c)]. The transition embedment depth defining the two distinct zones of behavior A and B first observed in Fig. 3 is not altered, and the dimensionless lift force remains constant at embedment depths of about H > 2.2D. Note that the drag force increasing with the macrofriction angle echoes the finding in the DEM simulations presented by Zhou et al. [13], who observed that the drag forces increase with the microfriction coefficient between the granules and thus the macrofriction angle (since the latter depends significantly on the former). Observe also in Fig. 5(c) that the rate of increase of the dimensionless drag force increases with increasing macrofriction angle. This is due to the change in the shape of the Mohr-Coulomb yield surface, and the corresponding increase in the resistance of the granular matter to shearing at higher confining stresses. Lastly, the drag force appears to be rather insensitive to the value of the dilation angle, which merely controls the
022204-4
15
-10
10
-5
5
0
0
Critical state drag force, ND=FD/γHDL
0 -30
5
10
H/D
15
ND (Φ = 25 ) ND (Φ = 30o) ND (Φ = 35o) NL (Φ = 25o) NL (Φ = 30o) NL (Φ = 35o) o
-25 -20 -15
30
(c)
25 20 15
-5
5
0 10
H/D
15
-15
20
0
25 20 15
-10
10
-5
5
0
0 0
10
5
-20
20
-10
0
-25
30
(b)
5
-30
10
15
H/D
ND (Ψ = 5 ) ND (Ψ = 10o) ND (Ψ = 15o) NL (Ψ = 5o) NL (Ψ = 10o) NL (Ψ = 15o) o
-25 -20 -15
20 30
(d)
25 20 15
-10
10
-5
5
0
0 0
5
10
15
H/D
Critical state lift force, NL=FL/γHDL
-15
20
ND (μ = tan(0.25Φ)) ND (μ = tan(0.5Φ)) ND (μ = tan(0.75Φ)) NL (μ = tan(0.25Φ)) NL (μ = tan(0.5Φ)) NL (μ = tan(0.75Φ))
Critical state lift force, NL=FL/γHDL
-20
25
-30
Critical state drag force, ND=FD/γHDL
-25
30
(a)
Critical state drag force, ND=FD/γHDL
ND (ρ = 1.0 g/cm3) ND (ρ = 2.0 g/cm3) ND (ρ = 3.0 g/cm3) NL (ρ = 1.0 g/cm3) NL (ρ = 2.0 g/cm3) NL (ρ = 3.0 g/cm3)
Critical state lift force, NL=FL/γHDL
-30
PHYSICAL REVIEW E 91, 022204 (2015)
Critical state lift force, NL=FL/γHDL
Critical state drag force, ND=FD/γHDL
NUMERICAL INVESTIGATION OF THE CYLINDER . . .
20
FIG. 5. (Color online) Variation of the dimensionless critical-state drag force ND and lift force NL with the embedment depth. (a) Effect of density ρ with φ = 30◦ , ψ = 10◦ , μ = tan(0.5φ); (b) effect of surface roughness μ with ρ = 2.0 g/cm3 , φ = 30◦ , ψ = 10◦ ; (c) effect of friction angle φ with ρ = 2.0 g/cm3 , ψ = 10◦ , μ = tan(0.5φ); (d) effect of dilation angle ψ with ρ = 2.0 g/cm3 , φ = 30◦ , μ = tan(0.5φ).
-15
IV. DRAG AND LIFT FORCES DEVELOPING UNDER LARGE HORIZONTAL DISPLACEMENTS
15 ND NL
-10
In the previous section we investigated critical-state drag and lift forces at the instant where a complete failure surface is formed within the granular matter. Now we increase gradually the applied drag displacement on the cylinder ux from about 0.02D to 15D, to explore the evolution of the drag and lift forces at very large deformations by means of the PFEM. Initially we consider the cylinder being dragged at various embedment depths in loose Cornell filter sand (Table I), with a zero dilation angle, i.e., deforming at constant volume upon shearing. Dimensionless force-displacement curves for the analysis for embedment depth H = 4D are plotted in Fig. 6. Interestingly, the dimensionless drag force increases from the
10
Steady state
Critical State
-5
5
0
0 0
5
ux/D
10
Critical state lift force, NL=FL/γHDL
critical-state value ND = 7.2 (see also Fig. 3) to ND = 8, when a steady state is achieved at approximately ux = 10D. In contrast, the dimensionless lift force NL remains constant. To shed some light on this, we plot in Fig. 7 the evolution of displacement increment contours and of the deformed material free surface, with increasing drag displacement. This shows
Critical state drag force, ND=FD/γHDL
magnitude of plastic volume expansion upon shearing. The dilation angle does not greatly influence the shear resistance of the granular material in unrestrained problems, when volume expansion is not suppressed by the boundary conditions. We ensure this here by controlling the applied horizontal displacement, and therefore the results are somewhat expected. This is not the case though for the lift force, which increases with increasing dilation angle and associated volume expansion of the granular matter within the failure surface, which drifts together with the cylinder towards the free surface.
15
FIG. 6. (Color online) Variation of the dimensionless drag force ND and lift force NL with the applied displacement. Results correspond to cylinder embedment H = 4D and material dilation angle ψ = 0◦ .
022204-5
ZHANG, SHENG, KOURETZIS, KRABBENHOFT, AND SLOAN
(a)
PHYSICAL REVIEW E 91, 022204 (2015)
ux = 0.1D
(a) Displacement increment (cm) 0.1
Displacement Increment (cm) 0.1
ux = 0.1D
0
0
0.2
0.2
(b)
ux = 0.8D
(b)
(a)
ux = 0.8D
ux = 10D
(c)
ux = 15D
(d)
(c)
ux = 10D (d) FIG. 8. (Color online) Evolution of displacement increment contours illustrating flow of granules within the failure surface with increasing drag displacement ux : (a) ux = 0.1D, (b) ux = 1D, (c) ux = 10D, (d) ux = 15D. Results correspond to cylinder embedment H = 8D and material dilation angle ψ = 8◦ .
ux = 15D FIG. 7. (Color online) Evolution of displacement increment contours illustrating flow of granules within the failure surface with increasing drag displacement ux : (a) ux = 0.1D, (b) ux = 1D, (c) ux = 10D, (d) ux = 15D. Results correspond to cylinder embedment H = 4D and material dilation angle ψ = 0◦ .
that the surface of the granular bed deforms to accommodate the movement of the cylinder, and a heap is formed as granular material is accumulated upstream of the cylinder. Formation of the heap increases the embedment depth in the vicinity of the cylinder [Figs. 7(a)–7(c)] and, therefore, the dimensionless drag force, which is proportional to the depth (see also Fig. 3). At large horizontal displacements (ux > 10D) a steady state is reached, as the volume of the heap is limited by the shear strength of the material [Figs. 7(c) and 7(d)]. Unlike the drag force, the dimensionless lift force NL does not depend on the embedment depth, as discussed earlier (see also Fig. 3), so the steady-state lift force is equal to the critical-state lift force. At larger embedment depths, where a local failure mode will prevail and the failure surface does not propagate to the surface of the granular bed, the steady-state drag force will also be equal to the corresponding critical-state value. The above discussion focused on a material that does not exhibit volume expansion during dragging (ψ = 0◦ ). It is worth investigating the effect of a small amount of dilation (ψ < 8◦ ) on the drag and lift forces, albeit considering con-
stant material dilation with the rigid, plastic failure criterion employed suggests that the granular matter will continue to dilate indefinitely, as long as shearing continues. This type of behavior is somewhat fictitious; in reality most granular materials will reach a state where they will deform at constant volume upon further shearing (ψ = 0◦ ). However, insights gained from this analysis contribute to the understanding of the mechanisms governing cylinder-granular matter interaction at large horizontal displacements. Figure 8 is similar to Fig. 7, but resulted from an analysis for dilation angle ψ = 8◦ . Volume expansion upon continuous shearing results in the formation of an unrealistic heap. This is also depicted in Fig. 9, presenting the variation of the dimensionless drag and lift forces with increasing horizontal displacement of the cylinder. Drag forces increase following the increase of the embedment of the cylinder due to the accumulation of the expanding granular matter. The rate of this increase is a function of the value of the dilation angle; higher dilation angles will result in larger heap volumes, and thus higher increase rates. This is analogous to the linear increase of the deviatoric stress with continued undrained shearing in triaxial test simulations with the Mohr-Coulomb model and a dilation angle ψ > 0◦ . On the other hand, dimensionless lift forces remain constant and equal to the critical-state value, which suggests that the mechanism for the development of the lift forces is related to the gradient of the gravitational
022204-6
NUMERICAL INVESTIGATION OF THE CYLINDER . . .
20
ψ=0 ψ=4 ψ=8
-15
15
ND
-10
10
NL
-5
5
0
0 0
5
ux/D
10
APPENDIX: SOLUTION ALGORITHM
Critical state lift force, NL=FL/γHDL
Critical state drag force, ND=FD/γHDL
-20
PHYSICAL REVIEW E 91, 022204 (2015)
1. Time discretization
The momentum conservation equations (1) can be discretized using the θ method as vn+1 − vn , (A1)
t un+1 − un , (A2) θ2 vn+1 + (1 − θ2 )vn =
t where v are velocities, subscripts n and n + 1 refer to the known and new, unknown states, and t = tn+1 − tn is the time step. Rearranging the above equations leads to ∇ T [θ1 σ n+1 + (1 − θ1 )σ n ] + b = ρ
15
u ∇ T σ n+1 + b˜ = ρ˜ 2 ,
t 1 u − (1 − θ2 )vn , vn+1 = θ2 t
FIG. 9. (Color online) Variation of the dimensionless drag force ND and lift force NL developing on the cylinder which moves in granular matter with different dilation angles. Results correspond to cylinder embedment H = 8D.
(A3) (A4)
where u = un+1 − un and ρ , θ1 θ2 1 − θ1 T vn 1 + ∇ σ n, b˜ = b + ρ˜ θ1
t θ1
stresses, and drifting of the cylinder towards the free surface due to volume expansion.
ρ˜ =
V. CONCLUSIONS
Use of the recently developed particle finite element method (PFEM) allowed us to study the mechanisms governing horizontal movement of a rigid cylinder inside a bed of granular matter. This method overcomes the technical restrictions of simple physical models and the shortcomings of classical numerical approaches. Some interesting insights emerging from this study include the numerical documentation of the development of a lift force on the cylinder. After a transition embedment depth of about H/D = 2.2 (which is related to the shape of the failure surface but not to the mechanical properties of the granular matter), the dimensionless lift force becomes independent of the embedment depth and the magnitude of the applied cylinder horizontal drag, and its value is only a function of the shear strength and volume expansion potential of the granular matter. On the other hand, the dimensionless drag force increases linearly with the embedment of the cylinder after the transition embedment depth H/D = 2.2, in line with published experimental results, and its critical-state value is a function of the macrofriction angle of the material. As the applied cylinder horizontal drag increases beyond the critical displacement, which is the horizontal displacement required for the complete failure surface to develop within the granular matter, the dimensionless drag force rises up to a steady-state value. This phenomenon is attributed to the accumulation of granular material above the cylinder, which might take place even when the granular material does not exhibit a tendency for volume expansion upon continuous shearing. The particle finite element method (PFEM) can be further adapted to incorporate more complex aspects of granular material behavior such as a material’s volume expansion or contraction being dependent on the applied shearing, or even extended to three dimensions, to treat the nonsymmetric problem of stirring of granular matter.
(A5) (A6)
The natural boundary conditions (2) are approximated in an analogous manner leading to NT σ n+1 = ˜t, on S,
(A7)
where ˜t =
1 1 − θ1 T t− N σ n. θ1 θ1
(A8)
Following [20], the above problem can be stated in terms of a min-max problem: min max
u (σ ,r)n+1
subject to
˜ ˜ σ n+1 ,∇( u)V − b, u V − t, uS − 12 t 2 rn+1 ,ρ˜ −1 rn+1 V + rn+1 , uV
(A9)
F (σ n+1 ) 0
on the basis of the Hellinger-Reissner variational principle. In the above, rn+1 are a set of variables interpreted as dynamic forces [20], and the notation x,yA = xT ydA (A10) A
is utilized. The equivalence between the optimization problem (A9) and the governing equations at hand is established by demonstrating that the Euler-Lagrange equations associated with (A9) indeed reproduce the governing equations [20]. 2. Spatial discretization
Using the standard finite element notation, the following approximations,
022204-7
σ (x) ≈ Nσ σˆ ,
(A11)
r(x) ≈ Nr rˆ ,
(A12)
ˆ u(x) ≈ Nu u,
(A13)
ˆ n∇u ≈ Bu u,
(A14)
ZHANG, SHENG, KOURETZIS, KRABBENHOFT, AND SLOAN
PHYSICAL REVIEW E 91, 022204 (2015)
are introduced for the state variables, where σˆ , rˆ , and uˆ are the nodal variables, N matrices contain the shape functions, and Bu = ∇Nu . Substituting the above approximations into the variational principle (A9) results in the following discrete principle:
At this stage, the contact constraints (6) are imposed on all potential contact nodes (i.e., mesh nodes located on the boundaries), which leads to a final problem of the type
min max
σˆ Tn+1 B uˆ
subject to
1 − t 2 rˆ Tn+1 Dˆrn+1 + rˆ Tn+1 A uˆ 2 j F σˆ n+1 0, j = 1, . . . ,nσ
uˆ (σˆ ,ˆr)n+1
− f uˆ
V
V
subject to (A15)
V
S
NTu ˜tdS,
NTr ρ˜ −1 Nσ dV ,
(A17) (A18)
A= V
NTu Nr dV .
(A19)
Then, solving the minimization part of (A15) gives a maximization problem as maximize (σˆ ,ˆr)n+1
subject to
1 − t 2 rˆ Tn+1 Dˆrn+1 2 BT σˆ n+1 + AT rˆ n+1 = f j F σˆ n+1 0, j = 1, . . . ,nσ .
nc 1 − t 2 rˆ Tn+1 Dˆrn+1 − g0j pj 2 j =1 BT σˆ n+1 + AT rˆ n+1 + ET ρ = f F σˆ in+1 0, i = 1, . . . ,nσ
pj = −nT ρj ,
˜ + NTu bdV
D=
(σˆ ,ˆr,pj )n+1
T
where nσ is the number of Gauss integration points, and NTσ Bu dV , (A16) B= f=
maximize
(A20)
[1] R. Candelier and O. Dauchot, Phys. Rev. Lett. 103, 128001 (2009). [2] N. Daiyan, S. Kenny, R. Phillips, and R. Popescu, Can. Geotech. J. 48, 1683 (2011). [3] T. Hsu, Y. Chen, and W. Hung, J. Transp. Eng. 132, 175 (2006). [4] T.-W. Hsu, Can. Geotech. J. 33, 180 (1996). [5] K. Wieghardt, Annu. Rev. Fluid Mech. 7, 89 (1975). [6] A. Seguin, Y. Bertho, F. Martinez, J. Crassous, and P. Gondret, Phys. Rev. E 87, 012201 (2013). [7] Y. Amarouchene, J. F. Boudet, and H. Kellay, Phys. Rev. Lett. 86, 4286 (2001). [8] D. Chehata, R. Zenit, and C. R. Wassgren, Phys. Fluids (1994present) 15, 1622 (2003). [9] F. Guillard, Y. Forterre, and O. Pouliquen, Phys. Rev. Lett. 110, 138303 (2013). [10] F. Guillard, Y. Forterre, and O. Pouliquen, Phys. Fluids 26, 043301 (2014). [11] C. Trautmann and T. O’Rourke, J. Geotech. Eng. 111, 1077 (1985). [12] Y. Ding, N. Gravish, and D. I. Goldman, Phys. Rev. Lett. 106, 028001 (2011).
(A21)
j = 1, . . . ,nc
qj = −nˆ T ρj |qj | − μpj 0, where ρ = (ρ1 ,ρ2 )T are the nodal forces, n = (n1 ,n2 )T and nˆ = (−n2 ,n1 )T are the normal and the tangential of the rigid boundary, E is an index matrix of zeros and ones, and nc is the number of potential contacts. The above problem can be transformed into a standard form of the second-order cone program (SOCP) and then solved using the high-performance optimization solver MOSEK [27] with a default setup for the convergence tolerance. The transformation of (A21) into SOCP standard form is straightforward and has been documented in Refs. [18,19]. In the course of solving the problem, the kinematic variables (displacement increments and plastic multipliers) are recovered as the dual variables, or Lagrange multipliers, associated with the discrete equilibrium constraints. To date, the above scheme has been successfully applied to a variety of challenging problems including limit analysis [28], elastoplastic static analysis [18,19], rigid, plastic dynamic analysis with large deformations [20,21,24], and granular contact dynamic analysis [29].
[13] F. Zhou, S. G. Advani, and E. D. Wetzel, Phys. Fluids 19, 013301 (2007). [14] F. Q. Potiguar and Y. Ding, Phys. Rev. E 88, 012204 (2013). [15] J. Jung, T. O’Rourke, and N. Olson, J. Geotech. Geoenviron. Eng. 139, 2028 (2013). [16] G. P. Kouretzis, D. Sheng, and S. W. Sloan, Comput. Geotech. 54, 53 (2013). [17] S. Yimsiri, K. Soga, K. Yoshizaki, G. Dasari, and T. O’Rourke, J. Geotech. Geoenviron. Eng. 130, 830 (2004). [18] K. Krabbenhoft and A. V. Lyamin, Comput. Methods Appl. Mech. Eng. 209-212, 239 (2012). [19] K. Krabbenhøft, A. V. Lyamin, and S. W. Sloan, Int. J. Solids Struct. 44, 1533 (2007). [20] X. Zhang, K. Krabbenhoft, D. M. Pedroso, A. V. Lyamin, D. Sheng, M. V. da Silva, and D. Wang, Comput. Geotech. 54, 133 (2013). [21] X. Zhang, K. Krabbenhoft, and D. Sheng, Granular Matter 16, 609 (2014). [22] S. R. Idelsohn, E. O˜nate, and F. D. Pin, Int. J. Numer. Methods Eng. 61, 964 (2004). [23] E. O˜nate, S. R. Idelsohn, F. Del Pin, and R. Aubry, Int. J. Comput. Methods 01, 267 (2004).
022204-8
NUMERICAL INVESTIGATION OF THE CYLINDER . . .
PHYSICAL REVIEW E 91, 022204 (2015)
[24] X. Zhang, K. Krabbenhoft, D. Sheng, and W. Li, Comput. Mech. 55, 167 (2015). [25] P. Wriggers, Computational Contact Mechanics (Wiley, Chichester, UK, 2002). [26] G. P. Kouretzis, K. Krabbenhøft, D. Sheng, and S. W. Sloan, Can. Geotech. J. 51, 1087 (2014).
[27] E. D. Andersen, C. Roos, and T. Terlaky, Math. Program. 95, 249 (2003). [28] K. Krabbenhoft and L. Damkilde, Int. J. Numer. Methods Eng. 56, 165 (2003). [29] K. Krabbenhoft, J. Huang, M. V. da Silva, and A. V. Lyamin, Granular Matter 14, 607 (2012).
022204-9
