- Email: [email protected]
Particle finite element analysis of large deformation and granular flow problems
Computers and Geotechnics 54 (2013) 133–142
Contents lists available at SciVerse ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Particle finite element analysis of large deformation and granular flow problems X. Zhang a, K. Krabbenhoft a,⇑, D.M. Pedroso b, A.V. Lyamin a, D. Sheng a, M. Vicente da Silva c, D. Wang d a
Centre for Geotechnical Science and Engineering, University of Newcastle, NSW, Australia Department of Civil Engineering, University of Queenland, QLD, Australia c Departamento de Engenharia Civil, Universidade Nova de Lisboa, Portugal d Centre for Geotechnical Science and Engineering, University of Western Australia, WA, Australia b
a r t i c l e
i n f o
Article history: Received 16 April 2013 Received in revised form 1 July 2013 Accepted 2 July 2013 Available online 31 July 2013 Keywords: Particle Finite Element Method PFEM Large deformations
a b s t r a c t A version of the Particle Finite Element Method applicable to geomechanics applications is presented. A simple rigid-plastic material model is adopted and the governing equations are cast in terms of a variational principle which facilitates a straightforward solution via mathematical programming techniques. In addition, frictional contact between rigid and deformable solids is accounted for using an approach previously developed for discrete element simulations. The capabilities of the scheme is demonstrated on a range of quasi-static and dynamic problems involving very large deformations. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Although most geotechnical structures operate in the small deformation range, there are numerous problems within geotechnical engineering and related fields that call for methods where changes to the problem geometry as a result of deformations are taken into account. These include problems of landslides and debris flow, penetration of various devices such as cones and torpedo anchors into the ground, and the interaction, during installation or under operating conditions, of off-shore oil and gas infrastructure with the seabed. For some of these problems, the deformation pattern resembles a fluid flow more than a solid undergoing large distortions. In the framework on the standard finite element method, such problems give rise to two fundamental challenges. The first one relates to geometry in the sense that the magnitude of the deformations is bound to lead not only to severe mesh distortion, but also to situations where the boundaries of the problem change from one time step to the next. Of these two separate but related issues, the former has received by far the most attention. Indeed, for many problems the original boundaries are maintained even after relatively large distortions. The perhaps most common approach to avoiding or alleviating mesh distortion is the Arbitrary Eulerian Lagrangian (ALE) method. This method utilizes the respective advantages of pure Eulerian and pure Lagrangian formulations and has been used ⇑ Corresponding author. E-mail address: [email protected] (K. Krabbenhoft). 0266-352X/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2013.07.001
quite successfully in geotechnical applications [1,2] and other solid as well as fluid mechanics problems [3,4]. Another popular method for geotechnical applications is the so-called Remeshing and Interpolation Technique with Small Strain (RITSS) technique proposed by Hu and Randolph [5,6]. While both the ALE and the RITSS have been used to solve problems involving relatively large deformations, they both have shortcomings in the case where the original boundaries change in the course of the deformation process, for example in the case where an initially contiguous solid separates into two or more parts as a result of external actions. The second challenge, which in many ways is the more serious one (though it remains much less explored), is that of solving the governing equations – comprising momentum balance, strain–displacement relations and constitutive relations that usually are highly nonlinear and may give rise to ill-posedness, localization of deformations, etc. Indeed, it is well known that even small deformation problems are hard to deal with for constitutive models that involve nonassociated flow rules [7]. In this paper a new scheme that addresses both of the fundamental challenges described above is presented. The scheme is applicable to general large deformation problems with no real limitations on the magnitude of the deformations. In other words, both problems where the deformation patterns resemble fluid flows and those that merely involve the deformation of solids slightly beyond the small deformation limit can be handled. More specifically, issues related to geometry are handled by means of the Particle Finite Element Method (PFEM) [8–10] while the solution of the governing equations is addressed by means of variational and
134
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
mathematical programming methods that have their origin in computational limit analysis [11–16] but since have been applied to a wide range of other problems including elastoplasticity [17,18,7] and discrete element type analysis [19–21]. Other issues dealt with include dynamics and frictional contact. The paper is organized as follows. In Section 2, the fundamentals of the PFEM are briefly described with emphasis on the alpha-shape method used for identifying solid and void domains on the basis of a cloud of points. Section 3 details the governing equations and their variational formulation. Next, in Section 4, the discretization and solution of the governing equations are described before the treatment of contact is detailed in Section 5. Finally, in Section 6, a number of examples demonstrating the capabilities of the scheme are presented before conclusions are drawn in Section 7. While all aspects of the new scheme are applicable to the general three-dimensional setting, the examples given in this paper are limited to two dimensions assuming plane strain. Standard matrix notation is used throughout with bold upper and lower case letters denoting matrices and vectors respectively. 2. Particle Finite Element Method The Particle Finite Element Method (PFEM) is, despite its name, a mesh based continuum method. First developed for fluid dynamics applications [8–10], the PFEM makes use of a Lagrangian description to account for the motion of nodes of the finite element mesh. A key feature of the method is that nodes are viewed as free ‘particles’ that can separate from the solid to which they originally belong. On the basis of the resulting cloud of points, the solid and void domains are identified and a standard finite element discretization used to advance the simulation by a given time step. More specifically, considering a time step tn ? tn+1, the steps of the PFEM are as follows (see also Fig. 1): 0. A cloud of particles, Cn, is given at time tn. 1a. On the basis of Cn, identify the computational domain, Vn. 1b. Mesh the domain and discretize the governing equations on Mn. 1c. Map the state variables (velocities, stresses, etc.) from the old mesh, Mn1, to the new mesh, Mn. 2a. Solve the discrete governing equations to obtain the displacement of the nodes.
Fig. 1. Steps in the Particle Finite Element Method.
2b. Update the positions of the nodes to arrive at Cn+1 and repeat. 2.1. Alpha-shape method The critical issue in the steps above is the identification of the computational domain on the basis of a cloud of points. For the general case, there is no unique solution to this problem. The solution originally proposed by Idelsohn et al. [8] and subsequently adopted as a standard feature of the PFEM was to use the alphashape method previously developed by Edelsbrunner and Mucke [22] for computer graphics applications. The basic principle of the alpha-shape method is as follows. Consider a cloud of points with a characteristic spacing h. Then for some predefined value of a parameter a, all nodes on an empty sphere with a radius greater than ah are considered boundary nodes. In other words: for each point in the domain examine whether it is possible to place a sphere with radius ah such that it contains only that point. If possible, the point is a boundary point and if not, i.e. if the sphere inevitably will contain more than the one point, it is an internal point. While a number of algorithms for recognizing boundaries by means of the alpha-shape method are available, another, and in many ways more straightforward possibility, has been proposed by Cremonesi et al. [23]. The steps in this scheme are as follows. Consider the cloud of particles as shown in Fig. 2(a). A Delaunay triangulation is first performed to generate the convex domain shown in Fig. 2(b). Next, the radius of the circumcircle of each triangle is examined and triangles with a circumcircle radius greater than ah are deleted. For the example at hand, this leads to the final configuration shown in Fig. 2(c). As shown by Cremonesi et al., this procedure is equivalent to the original alpha-shape approach owing to the property that the circumcircles of all triangles generated by the Delaunay triangulation are empty. 2.1.1. Choice of a It is clear that the alpha-shape method involves an element of subjectivity. Indeed, the resulting domain is a direct function of the value of the parameter a. This is illustrated in Fig. 3. A value of a = 1.2 here produces a set of boundaries that in many cases would be deemed reasonable. In fact, any value of a in the interval 0.9 6 a 6 1.3 produces this set of boundaries. Decreasing a below 0.9 leads first to internal voids (a = 0.5) and subsequently to a disintegration of the external boundaries as well (a = 0.4). On the other hand, increasing a above 1.3 leads first to a coalescence of the two distinct sets of points (a = 1.5) and then, for large values of a, to a solid defined by the convex hull inscribing the cloud of points. The conclusion of this example, that a value of a slightly greater than 1 is appropriate, is consistent with experience from application to actual physical problems. For example, for a wide range of coupled fluid–solid interaction problems, Onate et al. [24] conclude that a value of a in the range of 1.3–1.5 is appropriate. For
Fig. 2. Boundary recognition via the scheme of Cremonesi et al. [23]: cloud of points (a), Delaunay triangulation (b), after deletion of triangles with circumcircle greater than ah (c).
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
135
3.2. Constitutive equations The constitutive relations assumed in the following are those of rigid-perfect plasticity. That is, no deformations take place below yield while, once the yield condition is satisfied, plastic deformation occur. In the static case, these are in principle unlimited while in the dynamic case, the extent of plastic straining is limited by the momentum balance equations. Moreover, the flow rule is assumed nonassociated. The constitutive conditions are thus given by
FðrÞ 6 0 e_ ¼ k_ rr GðrÞ _ rÞ ¼ 0; kFð
ð3Þ k_ P 0
where F is the yield function, G is the plastic potential, and k_ is a plastic multiplier. In this paper, both the yield conditions and plastic potential are given in terms of the Mohr–Coulomb relations which under plane strain conditions can be expressed as
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðr11 r33 Þ2 þ 4r213 þ ðr11 þ r33 Þ sin / 2c cos / qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi G ¼ ðr11 r33 Þ2 þ 4r213 þ ðr11 þ r33 Þ sin w F¼
Fig. 3. Boundary recognition via the alpha-shape method.
the examples of the present paper which cover static and dynamic deformation processes of purely cohesive and pure frictional solids, we find the appropriate range of a to be approximately 1.4–1.6. Thus, while the ‘optimal’ value of a is problem dependent, the range of possible values is in practice rather limited and requires a relatively minor effort to establish by simple trial and error. Nevertheless, regardless of the choice of a, the total volume, and thereby, for incompressible materials, the total mass, is bound to fluctuate somewhat in the course of a typical simulation. However, as will be demonstrated via examples, this error is rather small. In principle, it could be controlled in various ways, for example by adjusting a in the course of the simulations. At present no general scheme is available for this purpose although a particular scheme, applicable to fluid dynamics recently has been proposed [25].
3. Governing equations In the following, the governing equations of rate-independent elastoplasticity are briefly summarized. 3.1. Momentum balance and strain–displacement relations Assuming infinitesimal deformations, the strain–displacement relations are given by
e ¼ $u
ð1Þ
where u = (u1, u2, u3)T are the displacements, e = (e11, e22, e33, 2e12, 2e23, 2e31)T are the strains, and $ is the usual linear strain–displacement operator. The momentum conservation equations for a body subjected to a time-varying motion are given by
€; $T r þ b ¼ qu T
N r ¼ t;
in V
on S
ð2Þ
where r = (r11, r22, r33, r12, r23, r31)T are the stresses, V is the domain under consideration, S is its boundary, b are body forces, t are tractions, and N ¼ $ðnxT Þ with n = (n1, n2, n3)T being the outward normal to the boundary and x = (x1, x2, x3)T being the spatial coordinate. The material density is given by q and a superposed dot denotes differentiation with respect to time.
ð4Þ
where / is the friction angle, c is the cohesion, and w is the dilation angle (with w = 0 implying incompressibility). While this constitutive model is very simple, it is for many problems involving very large deformations probably quite appropriate. Indeed, the model verifies the conditions at the critical state and with the initial transient phase leading up to this state thus being ignored, it may be thought of as a reasonably accurate first-order model. Besides failing to capture the conditions prior to fully developed plastic flow, the only obvious flaw in the model is the absence of rate effects. While such effects usually can be ignored at the microscopic level, the processes at this level nevertheless manifest themselves in a what, at the macroscopic level, appears to be a rate effect, essentially implying an increase in friction angle with shear strain rate [26,27]. While such an enhancement is entirely possible within the general framework to be presented, it will not be considered in this paper. Finally, we note that the incremental response in each time step is based on infinitesimal strain theory. This may lead to several errors of which the generation of strains as a result of rigid body motions is the most serious one. However, based on standard convergence analysis, it appears that this and related errors are relatively minor for the kind of time steps used in typical simulations. As such, the price to pay for the convenience of being able to operate with usual infinitesimal strain theory appears to be very small. 4. Discretization and solution In this section, the discretization and solution of the governing equations are considered. The temporal discretization is first detailed after which we derive a time-discrete variational principle that lends itself to a straightforward spatial discretization resulting in a nonlinear mathematical program. We finally comment briefly on the solution of this program. 4.1. Time discretization Introducing velocities, v, the first part of the momentum balance equations may be written in terms of two coupled first-order equations:
$T r þ b ¼ qv_
v ¼ u_
ð5Þ
136
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
These equations are then approximated by means of the hmethod to yield:
$T ½h1 rnþ1 þ ð1 h1 Þrn þ b ¼ q h2 v nþ1 þ ð1 h2 Þv n ¼
v nþ1 v n
unþ1 un Dt
Dt
ð6Þ
where subscripts n and n + 1 refer to the known and new, unknown, states respectively and Dt = tn+1 tn is the time step. Rearranging gives the following equation in the stresses and displacements
Du ~ 2 $ rnþ1 þ b~ ¼ q Dt T
h1 h2 1 v 1 h1 T ~ ~ nþ b¼ bþq $ rn h1 Dt h1
ð8Þ
After the displacements at tn+1 have been determined, the velocities are computed from
1 Du ð1 h2 Þv n h2 Dt
ð9Þ
The static boundary conditions are approximated in an analogous manner leading to
N T ½h1 rnþ1 þ ð1 h1 Þrn ¼ t;
on S
ð10Þ
or:
N T rnþ1 ¼ ~t;
on S
ð11Þ
where
~t ¼ 1 t 1 h1 N T rn : h1 h1
ð12Þ
The above time integration scheme is unconditionally stable provided that h1 P 12 and h2 P 12. For h1 ¼ h2 ¼ 12 it coincides with the Newmark average acceleration scheme while, for h1 = 0, a fully explicit and conditionally stable scheme is obtained. In the following, we will only consider implicit and unconditionally stable schemes. Regarding the constitutive model, we require that all relations are satisfied at time tn+1:
Fðrnþ1 Þ 6 0
$ðDuÞ ¼ Dkrr Gðrnþ1 Þ DkFðrnþ1 Þ ¼ 0; Dk P 0:
ð13Þ
4.2. Variational formulation We are now in a position to postulate that the above time discrete governing equations can be stated in terms of the following time-discrete optimization problem, or variational principle,
min Du
max
ðr;rÞnþ1
~ 1 r nþ1 iV þ hr nþ1 ; DuiV DkFðrnþ1 Þ 1=2Dt 2 hr nþ1 ; q
~ Dui h~t; Dui hrnþ1 ; $ðDuÞiV hb; V S
~ 1 rnþ1 iV þ hr nþ1 ; DuiV 1=2Dt2 hr nþ1 ; q subject to Fðrnþ1 Þ 6 0
(
$T rnþ1 þ b~ ¼ rnþ1 ; in V N T rnþ1 ¼ ~t;
on S dJ 2 ~ 1 ¼ Dt q r nþ1 þ Du ¼ 0 dr nþ1 dJ ¼ $u Dkrr Fðrnþ1 Þ ¼ 0: drnþ1
4.3. Spatial discretization The principle Eq. (14) is discretized in space by postulating finite element approximations to the state variables r, r and u. Using standard finite element notation, we have
rðxÞ N r ðxÞr^ rðxÞ N r ðxÞr^ ^ ; $u Bu ðxÞu ^ uðxÞ N u ðxÞu
hx; yiA ¼
xT ydA
ð18Þ
^ are the nodal variables, the N matrices contain the ^; ^ r and u where r shape functions, and Bu ¼ $N u . In standard finite formulations, the stress shape functions in Nr are chosen as being continuous within the elements and discontinuous between elements while functions in Nu in addition are continuous between elements and usually one polynomial degree higher than those of Nr. The shape functions for the dynamics forces are taken to be identical to those of the displacements but are discontinuous between elements. Substituting the above approximations into the variational principle Eq. (14) leads to the following discrete principle
max
^ ;^rÞnþ1 ðr
^ Du
subject to
r^ Tnþ1 BDu^ f T Du^ 1=2^rTnþ1 D^rnþ1 þ ^rTnþ1 ADu^ ^ jnþ1 Þ 6 0; Fðr
j ¼ 1; . . . ; nr ð19Þ
where we have used the notation
Z
ð17Þ
Furthermore, by addition of appropriate penalty terms to the functional, the yield and plastic consistency conditions are recovered (see [17,28,7,18] for details). From the second set of conditions, we see that the variables rn+1 are the dynamic forces ~ Du=Dt 2 . The momentum balance equations are thus verir nþ1 ¼ q fied. However, from the last set of conditions we see that the flow rule is associated, i.e. the flow potential is equal to the yield potential F. This is a consequence of the variational formulation and is not possible to circumvent directly. In this paper we use the approach suggested in [7] of replacing the original yield function by an ‘effective’ one which, when used as a flow potential, yields the correct plastic strains. This approach was shown to produce reliable results for a range of static small deformation problems and in this paper we demonstrate that it is equally reliable for dynamic problems involving very large deformations.
min ð14Þ
ð16Þ
where Dk are Lagrange multipliers. Next, taking functional derivatives, the following equations appear:
dJ ¼ dDu
q
v nþ1 ¼
~ Dui h~t; Dui J ¼ hrnþ1 ; $ðDuÞiV hb; V S
ð7Þ
where Du = un+1 un and
q~ ¼
4.2.1. Euler–Lagrange equations The equivalence between the above optimization problem and the governing equations follows by showing that the Euler–Lagrange equations associated with Eq. (14) do indeed reproduce the governing equations. To this end, define the functional:
ð15Þ
where the yield function is imposed at a finite number of points nr and
A
and r is a set of variables whose significance will become apparent shortly.
B¼
Z V
N Tr Bu dV
ð20Þ
137
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
f ¼
Z V
D¼
Z V
A¼
Z V
~ þ N Tu bdV
Z S
N Tu ~tdS
ð21Þ
~ 1 N r dV N Tr q
ð22Þ
N Tu N r dV:
ð23Þ
The minimization part of Eq. (19) may be solved first to yield
^ nþ1 þ AT ^r nþ1 f ¼ 0 BT r
ð24Þ
which is included as a constraint in the remaining maximization part of the problem:
maximize ^ ;^r Þnþ1 ðr
1=2^rTnþ1 D^r nþ1
^ nþ1 þ AT ^rnþ1 ¼ f subject to BT r ^ jnþ1 Þ 6 0; Fðr
ð25Þ
j ¼ 1; . . . ; nr :
This is the problem solved in each time step. In the course of solving the problem, the displacement increments appear as the dual variables, or Lagrange multipliers, associated with the discrete equilibrium constraints. The above problem is similar in structure to those arising from limit analysis [13,14], static elastoplasticity [17,7,18], and granular contact dynamics [19–21].
Fig. 4. Contact specification.
and those of the rigid body. The edges in this triangulation then define the potential contacts. Considering each individual particle as a free particle at a distance g0 from the boundary and with potential contact forces p and q, a variational formulation of the corresponding frictional contact problem is given by [30,20,21]:
maximize
g0 p
subject to jqj lp 6 0; p P 0
where the standard Coulomb friction law has been used, l being the friction coefficient. This principle implies the following kinematics:
DuN ¼ g 0 k
ð27Þ
DuT ¼ k 4.4. Solution algorithm
ð26Þ
The optimization problems Eq. (25) generated by the scheme in each time step are solved by means of the second-order cone programming (SOCP) solver MOSEK [29]. Such solvers have become increasingly more powerful and robust in recent years and require a computational effort similar to conventional implicit Newton– Raphson based schemes. The SOCP standard form operates with a linear objective function [as opposed to the quadratic one in Eq. (25)] and with inequality constraints that are slightly different from those of the Mohr–Coulomb criterion. However, a transformation of Eq. (25) into SOCP standard form is straightforward as has been documented previously, e.g. [14,18].
where k P 0 is a Lagrange multiplier such that k(j q j lp) = 0 and DuN and DuT are the normal and tangential displacement increments respectively. This means that the sliding rule is associated in the sense that sliding is accompanied by a dilation of relative magnitude l. As in the continuum case, this dilation is not desirable and can rarely be justified physically. However, as discussed in detail in [30], the dilation arising at frictional contacts may be viewed as an artifact of the time discretization that gradually decreases as the time step is reduced. Considering all potential contacts and incorporating the above contact problem into the original problem Eq. (25) leads to a final problem of the type
5. Treatment of contact
maximize 1=2^rTnþ1 D^r nþ1
For general problems, proper treatment of frictional contact is an important component of the overall problem. In this paper, we will consider the contact between the deforming solid and rigid boundaries. Extension to general non-stationary rigid bodies is straightforward while the extension to contact with other deformable bodies, including self-contact, is less so although it in principle can be dealt with using the methodology described in the following. The treatment of the above mentioned types of contact can be accommodated within the general variational framework by using the recently proposed granular contact dynamics method [30,20,21]. Though originally conceived for particle dynamics of the discrete element type, the basic principles are equally applicable in a finite element context. The procedure is as follows. Consider a body discretized by finite elements as shown in Fig. 4. Initially, each boundary node of the system is viewed as a particle with zero radius. Next, the particles that are likely to come into contact with the rigid body (in Fig. 4 taken to be a single segment) are identified. This can be done in various ways as described in detail in [21]. The most general approach, and the one adopted in this paper, is to discretize the boundary of the rigid body and perform a Delaunay triangulation with respect to the boundary nodes of the deformable solid
^ ;^r Þnþ1 ðr
nc X g 0;j pj j¼1
^ nþ1 þ AT ^r nþ1 þ ET q ¼ f subject to B r T
^ jnþ1 Þ 6 0; Fðr pj ¼
nTj
j ¼ 1; . . . ; nr
ð28Þ
qj ; j ¼ 1; . . . nC
^ Tj qj ; qj ¼ n jqj j lpj 6 0;
j ¼ 1; . . . nC j ¼ 1; . . . nC T
with q = (q1, q2) being the nodal forces, n the normal of the rigid ^ ¼ ðn2 ; n1 ÞT . Furthermore, E is an index matrix of boundary and n zeros and ones and nC is the number of potential contacts. 6. Examples In the following, four examples demonstrating the capabilities of the new scheme are presented. For all examples, a reasonable value of a was first determined by trial and error for a coarse discretization (in both time and space) and for a very limited range of the parameter set involved in each problem. As discussed previously, this procedure is less than ideal but at present the only feasible one. Also, in practice, it is not particularly arduous and is usually completed in a fraction of the time that it takes to solve a given problem (using more refined temporal and spatial discretizations) with a fixed value of a. By this initial calibration, a value
138
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
of a = 1.4 was adopted for the first two examples, while a = 1.6 for the two last. In all examples, the gravitational acceleration is g = 9.8 m/s2. 6.1. Cylinder–soil interaction The first example concerns the interaction between a rigid weightless cylinder and a Tresca soil. The cylinder first penetrates vertically to a depth of 0.8 m after which it is displaced horizontally for a total distance of 2 m (see Fig. 5). This type of problem involving the combined horizontal and vertical motion of a rigid cylindrical object resembles that of a pipeline subjected to either external loads due to waves and current or to internal temperature and pressure loadings. The analysis conducted uses 280 displacement increments of magnitude equal to 0.01 m (80 for the penetration and 200 for the subsequent horizontal motion). The cylinder-soil interface is assumed purely frictional with a friction coefficient l = tan 30° ’ 0.577. Other parameters are as follows: soil density q = 1.6 g/cm3, undrained shear strength c = cu = 1 kPa, and cylinder diameter D = 1 m. The analysis is fully dynamic with a time step of 1s. Under these conditions, inertial effects are negligible and the simulation may be regarded as quasi-static. The results of the simulation in terms of load–displacement curves are shown in Fig. 6. During the initial penetration phase, the resultant horizontal force remains at zero while the vertical force increases approximately linearly up to around Fy = 18 kN/m. The cylinder is then displaced horizontally while being fixed in its vertical position of uy = 0.8 m. This change in motion brings about an immediate increase in horizontal force and a decrease in the vertical force. This is followed by a more gradual change of the forces until the motion becomes steady at around ux = 1.5 m. At this point, the resultant forces are Fx ’ 9 kN/m and Fy = 7.5 kN/ m respectively. Thus, an appreciable vertical force is required to prevent the cylinder moving towards the surface in response to the purely horizontal motion. These findings echo those of [31] where a vertically constrained cylinder was plowed through a granular material resulting in a similar ratio of vertical to horizontal force. Finally, we note that the fundamental change in geometry associated with the self-contact at the back of the cylinder between uH = 1.0 m and uH = 1.5 m is handled seamlessly and without resorting to any additional rules or procedures. 6.2. Accretionary wedge The next problem concerns the deformation of an accretionary wedge as shown in Fig. 7. This is a common model problem in structural geology used to study, among other things, mountain building via folding and thrusting events [32–35]. The wedge has a height of 8 cm and a length of 140 cm and is composed of a purely frictional material with / = 30° and w = 0 and a density of 2 g/cm3. The left vertical and bottom horizontal boundaries are rigid with friction coefficients of l = 30° and l = 15° respectively. The simulation is carried out by displacing the left wall horizontally with a velocity of 0.1 cm/s. A total of 1000 steps, with a constant time step of 0.16 s, were used to yield a total wall displacement of 16 cm. Snapshots of the deformations at different times are shown in Fig. 7. Contrary to what might be expected, the deformations do not proceed by evolution of one of more dominant shear bands. Rather, the location of the shear bands changes in a somewhat unpredictable manner: the classic V-band at the initial phases first propagates in the direction of loading (Du = 8.0 cm). Further loading then gives way to a new shear band much like the initial one (Du = 10 cm). This again gives way to another family of shear
Fig. 5. Interaction of a rigid cylinder with a Tresca soil.
bands further down the direction of loading. The end result (when viewed as a continuous process) is a seemingly random appearance of shear bands. This phenomenon has recently been studied in some detail by Mary et al. [36] using a semi-analytical slip line type model. The conclusion is here that the system displays deterministic chaos of a nature similar to more common chaotic dynamical system. Without repeating this analysis in the present context,
Normalized resultant force, F (kN/m)
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
139
20 15 FV
10 FH
5
Fig. 8. Granular column problem: initial column and final deposit.
0
0
−5 0
0.4
0.5
1.0
1.5
2.0 uH
0.8 uV
Displacment, u (m) Fig. 6. Horizontal and vertical resultant forces acting on the cylinder.
we do note that the behavior of the system displays a behavior consistent with that observed and analyzed by Mary et al.
6.3. Collapse of granular columns The next example concerns the collapse of granular columns. The setup is shown in Fig. 8. A column of granular material of width d0 and height h0 is allowed to collapse under the action of gravity by quickly removing the right wall. This type of experiment, which may be seen as a model landslide, has received significant attention in recent years, both from an experimental point of view [37–39] as well as in terms of reproducing the experiments using both particle methods [40–45,30,21] and continuum models [45–47]. The two-dimensional plane strain version of the experiment considered in the following is due to Lube et al. [38]. On the basis of a rather extensive experimental program, the following fits to
the height, h1, and diameter, d1, of the final deposit were determined:
h1 ¼ d0
a;
a 6 1:15 2
1:1a5 ; a > 1:15
8 a < 1:8 > 1:6a; d1 d0 < ¼ transition region; 1:8 6 a 6 2:8 > d0 2 : 2:2a3 ; a > 2:8
ð29Þ
ð30Þ
where a = h0/d0 is the initial aspect ratio. We note that there is a transition region at 1.8 6 a 6 2.8 for which no expressions for the final width were given. To validate the proposed scheme, a series of experiments involving columns with different aspect ratios were conducted. In all simulations, the following parameters were used: density q = 2 g/cm3, cohesion c = 0, friction angle / = 31° (as quoted by Lube et al.) or / = 35°, and dilation angle w = 0. The proceed from time psimulations ffiffiffiffiffiffiffiffiffiffi t = 0 and are terminated at t ¼ t= h0 =g ¼ 4:0. Fig. 9 shows three columns with aspect ratios a = 0.5, 1.0, and 7.0 at different times instants of the collapse process using a time step of Dt ¼ 0:01. The internal deformation patterns observed, in particular the apparent ‘buckling’ of internal columns, corresponds very well with what can be observed experimentally and in DEM type simulations [30]. As expected, and predicted from Eq. (29), the two shortest
Fig. 7. Accretionary wedge problem. Colors are proportional to the equivalent shear strain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
140
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
100 time steps
200 time steps
400 time steps
Fig. 9. Collapse of granular columns (dots represent the finite element nodes).
3.0
Normalized final height, h∞ /d0
the model of Andrade et al. [27], better fits to the experiments can be expected. To gauge the influence of the time step, the a = 1 simulation was run with different time steps. The resulting final deposits are shown in Fig. 11. As seen, the shape and dimensions of the final deposit are relatively insensitive to the time step. Indeed, the final 2
Volume change (%)
columns collapse only partially while the spread of the tallest one is much wider. A more complete comparison with the experimental fits of Lube et al. is shown in Fig. 10. We here see that the columns essentially lack resistance, even for a friction angle of 35°. That is, the final height is underpredicted while the final run-out distance is overpredicted, and with the deviation between experiment and simulation increasing with increasing initial aspect ratio. As discussed in Section 3.2, the key feature missing in the constitutive model adopted is a rate-dependence of the friction angle on the shear strain rate. With such an enhancement, for example following
Fig. 11. Convergence of final deposit as function of time step for a = 1 (dots represent the finite element nodes).
2.5
0
−2 100 time steps
−4
200 time steps 400 time steps
−6
2.0
0
1
2
3 −
Normalized time, t
1.5
Fig. 12. Volume change as function of time for a = 1.
1.0
φ = 31°
0.5 0
φ = 35° 0
2
4
6
8
10
12
Normalized final width, (d∞ -d0 )/d0
Aspect ratio, a = h0 /d0 12 100 cm
10 8 6 4
φ = 31°
2 0
35 cm
φ = 35° 0
2
4
6
8
10
12
8 cm
Aspect ratio, a = h0 /d0 50 cm Fig. 10. Normalized final height and width of granular columns as function of aspect ratio (dots represent the finite element nodes).
Fig. 13. Silo discharge problem.
4
141
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
0.0 s
3.5 s
1.0 s
4.0 s
0.0 s
3.0 s
6.0 s
7.0 s
5.0 s
2.0 s
4.5 s
8.0 s
Fig. 15. Discharge of rough walled silo (dots represent the finite element nodes).
Fig. 14. Discharge of smoothed walled silo (dots represent the finite element nodes).
7
6.4. Discharge of silo The final example concerns the discharge of a two-dimensional silo with the dimensions shown in Fig. 13. The material parameters are: density q = 1.5 g/cm3, cohesion c = 0, friction angle / = 31°, and dilation angle w = 0°. The solid mass confined by the silo wall is initially at rest after which, at time t = 0, the bottom wall is removed. Two different sets of simulations were run: assuming perfectly smooth walls and assuming perfectly rough walls (corresponding to a wall-solid friction coefficient of l = tan /). Snapshots of the discharge in the two cases are shown in Figs. 14 and 15. As can be seen, the wall-solid interface properties have a significant influence on the overall flow pattern. Again, these patterns are consistent with those observed in both experiments and particle simulations [48]. Moreover, the wall-solid interface conditions influences the total discharge time, from approximately 5 s for the smooth silo to almost twice that for the rough silo.
Δ t = 0.04s
Volume in silo (1000 cm3 )
height changes very little between the three different time steps while the run-out distance does increase somewhat between the first two simulations involving 100 and 200 time steps respectively. A similar trend is observed for other aspect ratios. Finally, we investigate the mass balance properties of the scheme as function of the time step. The results, shown in Fig. 12, reveal the loss of mass to be around 2–6% and decreasing as the time step is decreased (the number of time steps increased).
6
Δ t = 0.02s
5
Δ t = 0.01s Δ t = 0.005s
4 3 2 1 0
0
1
2
3
4
5
6
7
Time (s) Fig. 16. Silo volume versus time for different time steps.
The convergence of the results in terms the discharge versus time for different time steps is shown in Fig. 16 for the smooth walled silo. As in the previous example, some 500 time steps (corresponding to Dt = 0.01) appears to be sufficient to ensure a reasonable degree of accuracy. This trend is the same for the rough walled silo. 7. Conclusions The variant of the Particle Element Method has been developed and applied to a range of geomechanics problems involving very large deformations. A variational approach has been taken whereby
142
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
the discrete governing equations give rise to a second-order cone program that can be solved in an efficient and robust manner using available tools. Moreover, frictional contact is accounted for in a straightforward manner and does not give rise to any complications. A number of examples have been presented and while the results generally are encouraging, there are several issues that warrant further consideration. Firstly, regarding the alpha-shape method which is at the heart of the PFEM. As discussed, the parameter a is at present chosen on the basis of simple trial and error. Clearly, a more rigorous or at least automated procedure would be highly desirable. This could well involve a parameter that would vary throughout the simulation in contrast to the constant value currently used for the entire simulation. Secondly, while the simple rate-independent model used in this paper in many cases probably is quite adequate, it has its limitations. This is most apparent in the granular column collapse simulations where the deviation between simulation and experiment increases as the column height, and thereby the characteristic velocity, increases. A rate-dependent model, essentially implying an increase in friction angle with velocity, is here likely to improve the results.
References [1] Nazem M, Sheng DC, Carter JP, Sloan SW. Arbitrary Lagrangian–Eulerian method for large-strain consolidation problems. Int J Numer Anal Meth Geomech 2008;32:1023–50. [2] Nazem M, Carter JP, Airey D. Arbitrary Lagrangian–Eulerian method for dynamic analysis of geotechnical problems. Comput Geotech 2009;36:549–57. [3] Belytschko T, Liu WK, Moran B. Nonlinear finite elements for continua and structures. Wiley; 2000. [4] Donea J, Huerta A, Ponthot JP, Rodriguez-Ferran A. Arbitrary Lagrangian Eulerian method. In: Stein E, de Borst R, Hughes TJR, editors. Encyclopedia of Computational Mechanics, vol. 1; 2004. [5] Hu Y, Randolph MF. A practical numerical approach for large deformation problems in soil. Int J Numer Anal Meth Geomech 1998;22:327–50. [6] Yu L, Hu YX, Liu J, Randolph MF, Kong XJ. Numerical study of spudcan penetration in loose sand overlying clay. Comput Geotech 2012;46:1–12. [7] Krabbenhoft K, Karim MR, Lyamin AV, Sloan SW. Associated computational plasticity schemes for nonassociated frictional materials. Int J Numer Meth Eng 2012;89:1089–117. [8] Idelsohn SR, Onate E, Del Pin F. The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. Int J Numer Meth Eng 2004;61:964–89. [9] Onate E, Idelsohn SR, Del Pin F, Aubry R. The particle finite element method an overview. Int J Comput Meth 2004;1:267–307. [10] Onate E, Idelsohn SR, Celigueta MA, Rossi R. Advances int the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows. Comput Meth Appl Mech Eng 2008;197:1777–800. [11] Lyamin AV, Sloan SW. Lower bound limit analysis using non-linear programming. Int J Numer Meth Eng 2002;55:573–611. [12] Lyamin AV, Sloan SW. Upper bound limit analysis using linear finite elements and non-linear programming. Int J Numer Anal Meth Geomech 2002;26:181–216. [13] Krabbenhoft K, Damkilde L. A general nonlinear optimization algorithm for lower bound limit analysis. Int J Numer Meth Eng 2003;56:165–84. [14] Krabbenhoft K, Lyamin AV, Sloan SW. Formulation and solution of some plasticity problems as conic programs. Int J Solids Struct 2007;44:1533–49. [15] Vicente da Silva M, Antao AN. A nonlinear programming method approach for upper bound limit analysis. Int J Numer Meth Eng 2007;72:1192–218. [16] Vicente da Silva M, Antao AN. Upper bound limit analysis with a parallel mixed finite element formulation. Int J Solids Struct 2008;45:5788–804. [17] Krabbenhoft K, Lyamin AV, Sloan SW, Wriggers P. An interior-point method for elastoplasticity. Int J Numer Meth Eng 2007;69:592–626.
[18] Krabbenhoft K, Lyamin AV. Computational Cam clay plasticity using secondorder cone programming. Comput Meth Appl Mech Eng 2012;209– 212:239–49. [19] Krabbenhoft K, Lyamin AV, Huang J, Vicente da Silva M. Granular contact dynamics using mathematical programming methods. Comput Geotech 2012;43:165–76. [20] Krabbenhoft K, Huang J, Vicente da Silva M, Lyamin AV. Granular contact dynamics with particle elasticity. Granul Matter 2012;14:607–19. [21] Huang J, Vicente da Silva M, Krabbenhoft K. Three-dimensional granular contact dynamics with rolling resistance. Comput Geotech 2013;49:289–98. [22] Edelsbrunner H, Mucke EP. Three dimensional alpha shapes. ACM Trans Graph 1994;13(1):43–72. [23] Cremonesi M, Frangi A, Perego U. A lagrangian finite element approach for the analysis of fluid-structure interaction problems. Int J Numer Meth Eng 2010;84:610–30. [24] Onate E, Idelsohn SR, Celigueta MA, Rossi R, Marti J, Carbonell JM, et al. Advances in the particle finite element method (PFEM) for solving coupled problems in engineering. In: Particle-based methods fundamentals and applications computational methods in applied sciences, vol. 25. Springer; 2011. p. 47–54. [25] Ryzhakov P, Onate E, Rossi R, Idelsohn S. Improving mass conservation in simulation of incompressible flows. Int J Numer Meth Eng 2012;90:1435–51. [26] Forterre Y, Jop P, Pouliquen O. A constitutive law for dense granular flows. Nature 2006;441:727–30. [27] Andrade JE, Chen Q, Lee PH, Avila CF, Evans TM. A constitutive law for dense granular flows. J Mech Phys Solids 2012;60:1122–36. [28] Krabbenhoft K. A variational principle of elastoplasticity and its application to the modeling of frictional materials. Int J Solids Struct 2009;46:464–79. [29] Andersen ED, Roos C, Terlaky T. On implementing a primal–dual interior-point method for conic quadratic optimization. Math Program 2003;95:249–77. [30] Krabbenhoft K, Lyamin AV, Huang J, Vicente da Silva M. Granular contact dynamics using mathematical programming methods. Comput Geotech 2012;43:165–76. [31] Ding Y, Gravish N, Goldman DI. Drag induced lift in granular media. Phys Rev Lett 2011;106:028001. [32] Souloumiac P, Leroy YM, Maillot B, Krabbenhoft K. Predicting stress distributions in fold-and-thrust belts and accretionary wedges by optimization. J Geophys Res 2009;114:B09404. [33] Souloumiac P, Krabbenhoft K, Leroy YM, Maillot B. Failure in accretionary wedges with the maximum strength theorem: numerical algorithm and 2d validation. Comput Geosci 2010;14:793–811. [34] Cubas N. Prediction of thrusting sequences in accretionary wedges. J Geophys Res 2008;113:B12412. [35] Maillot B, Leroy Y. Kink-fold onset and development based on the maximum strength theorem. J Mech Phys Solids 2006;54:2030–59. [36] Mary BCL, Maillot B, Leroy YM. Deterministic chaos in frictional wedges revealed by convergence analysis. Int J Numer Anal Meth Geomech 2013. http://dx.doi.org/10.1002/nag.2177 [in press]. [37] Lube G, Huppert HE, Sparks RSJ, Hallworth MA. Axisymmetric collapses of granular columns. J Fluid Mech 2004;508:175–99. [38] Lube G, Huppert HE, Sparks RSJ, Freundt A. Collapses of two-dimensional granular columns. Phys Rev E 2005;72:041301. [39] Lajeunesse E, Mangeney-Castelneau A, Vilotte JP. Spreading of a granular mass on an horizontal plane. Phys Fluids 2004;16:2731–8. [40] Lacaze L, Phillips JC, Kerswell RR. Planar collapse of a granular column: experiments and discrete element simulations. Phys Fluids 2008;20:063302. [41] Lacaze L, Kerswell RR. Axisymmetric granular collapse: a transient 3D flow test of viscoplasticity. Phys Rev Lett 2009;102:108305. [42] Zenit R. Computer simulations of the collapse of a granular column. Phys Fluids 2005;17:031703. [43] Staron L, Hinch EJ. Study of the collapse of granular columns using twodimensional discrete-grain simulation. J Fluid Mech 2005;545:1–27. [44] Staron L, Hinch EJ. The spreading of a granular mass: role of grain properties and initial conditions. Granul Matter 2007;9:205–17. [45] Lagree PY, Staron L, Popinet S. The granular column collapse as a continuum: validity of a two-dimensional Navier–Stokes model with a l(i)-rheology. J Fluid Mech 2011;686:378–408. [46] Mangeney-Castelnau A, Bouchut F, Vilotte JP, Lajeunesse E, Aubertin A, Pirulli M. On the use of saint venant equations to simulate the spreading of a granular mass. J Geophys Res 2005;110:B09103. [47] Kerswell RR. Dam break with Coulomb friction: a model for granular slumping? Phys Fluids 2005;17:057101. [48] Gonzalez-Montellano C, Ayuga F, Ooi JY. Discrete element modelling of grain flow in a planar silo: influence of simulation parameters. Granul Matter 2011;13:149–58.
Contents lists available at SciVerse ScienceDirect
Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo
Particle finite element analysis of large deformation and granular flow problems X. Zhang a, K. Krabbenhoft a,⇑, D.M. Pedroso b, A.V. Lyamin a, D. Sheng a, M. Vicente da Silva c, D. Wang d a
Centre for Geotechnical Science and Engineering, University of Newcastle, NSW, Australia Department of Civil Engineering, University of Queenland, QLD, Australia c Departamento de Engenharia Civil, Universidade Nova de Lisboa, Portugal d Centre for Geotechnical Science and Engineering, University of Western Australia, WA, Australia b
a r t i c l e
i n f o
Article history: Received 16 April 2013 Received in revised form 1 July 2013 Accepted 2 July 2013 Available online 31 July 2013 Keywords: Particle Finite Element Method PFEM Large deformations
a b s t r a c t A version of the Particle Finite Element Method applicable to geomechanics applications is presented. A simple rigid-plastic material model is adopted and the governing equations are cast in terms of a variational principle which facilitates a straightforward solution via mathematical programming techniques. In addition, frictional contact between rigid and deformable solids is accounted for using an approach previously developed for discrete element simulations. The capabilities of the scheme is demonstrated on a range of quasi-static and dynamic problems involving very large deformations. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction Although most geotechnical structures operate in the small deformation range, there are numerous problems within geotechnical engineering and related fields that call for methods where changes to the problem geometry as a result of deformations are taken into account. These include problems of landslides and debris flow, penetration of various devices such as cones and torpedo anchors into the ground, and the interaction, during installation or under operating conditions, of off-shore oil and gas infrastructure with the seabed. For some of these problems, the deformation pattern resembles a fluid flow more than a solid undergoing large distortions. In the framework on the standard finite element method, such problems give rise to two fundamental challenges. The first one relates to geometry in the sense that the magnitude of the deformations is bound to lead not only to severe mesh distortion, but also to situations where the boundaries of the problem change from one time step to the next. Of these two separate but related issues, the former has received by far the most attention. Indeed, for many problems the original boundaries are maintained even after relatively large distortions. The perhaps most common approach to avoiding or alleviating mesh distortion is the Arbitrary Eulerian Lagrangian (ALE) method. This method utilizes the respective advantages of pure Eulerian and pure Lagrangian formulations and has been used ⇑ Corresponding author. E-mail address: [email protected] (K. Krabbenhoft). 0266-352X/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2013.07.001
quite successfully in geotechnical applications [1,2] and other solid as well as fluid mechanics problems [3,4]. Another popular method for geotechnical applications is the so-called Remeshing and Interpolation Technique with Small Strain (RITSS) technique proposed by Hu and Randolph [5,6]. While both the ALE and the RITSS have been used to solve problems involving relatively large deformations, they both have shortcomings in the case where the original boundaries change in the course of the deformation process, for example in the case where an initially contiguous solid separates into two or more parts as a result of external actions. The second challenge, which in many ways is the more serious one (though it remains much less explored), is that of solving the governing equations – comprising momentum balance, strain–displacement relations and constitutive relations that usually are highly nonlinear and may give rise to ill-posedness, localization of deformations, etc. Indeed, it is well known that even small deformation problems are hard to deal with for constitutive models that involve nonassociated flow rules [7]. In this paper a new scheme that addresses both of the fundamental challenges described above is presented. The scheme is applicable to general large deformation problems with no real limitations on the magnitude of the deformations. In other words, both problems where the deformation patterns resemble fluid flows and those that merely involve the deformation of solids slightly beyond the small deformation limit can be handled. More specifically, issues related to geometry are handled by means of the Particle Finite Element Method (PFEM) [8–10] while the solution of the governing equations is addressed by means of variational and
134
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
mathematical programming methods that have their origin in computational limit analysis [11–16] but since have been applied to a wide range of other problems including elastoplasticity [17,18,7] and discrete element type analysis [19–21]. Other issues dealt with include dynamics and frictional contact. The paper is organized as follows. In Section 2, the fundamentals of the PFEM are briefly described with emphasis on the alpha-shape method used for identifying solid and void domains on the basis of a cloud of points. Section 3 details the governing equations and their variational formulation. Next, in Section 4, the discretization and solution of the governing equations are described before the treatment of contact is detailed in Section 5. Finally, in Section 6, a number of examples demonstrating the capabilities of the scheme are presented before conclusions are drawn in Section 7. While all aspects of the new scheme are applicable to the general three-dimensional setting, the examples given in this paper are limited to two dimensions assuming plane strain. Standard matrix notation is used throughout with bold upper and lower case letters denoting matrices and vectors respectively. 2. Particle Finite Element Method The Particle Finite Element Method (PFEM) is, despite its name, a mesh based continuum method. First developed for fluid dynamics applications [8–10], the PFEM makes use of a Lagrangian description to account for the motion of nodes of the finite element mesh. A key feature of the method is that nodes are viewed as free ‘particles’ that can separate from the solid to which they originally belong. On the basis of the resulting cloud of points, the solid and void domains are identified and a standard finite element discretization used to advance the simulation by a given time step. More specifically, considering a time step tn ? tn+1, the steps of the PFEM are as follows (see also Fig. 1): 0. A cloud of particles, Cn, is given at time tn. 1a. On the basis of Cn, identify the computational domain, Vn. 1b. Mesh the domain and discretize the governing equations on Mn. 1c. Map the state variables (velocities, stresses, etc.) from the old mesh, Mn1, to the new mesh, Mn. 2a. Solve the discrete governing equations to obtain the displacement of the nodes.
Fig. 1. Steps in the Particle Finite Element Method.
2b. Update the positions of the nodes to arrive at Cn+1 and repeat. 2.1. Alpha-shape method The critical issue in the steps above is the identification of the computational domain on the basis of a cloud of points. For the general case, there is no unique solution to this problem. The solution originally proposed by Idelsohn et al. [8] and subsequently adopted as a standard feature of the PFEM was to use the alphashape method previously developed by Edelsbrunner and Mucke [22] for computer graphics applications. The basic principle of the alpha-shape method is as follows. Consider a cloud of points with a characteristic spacing h. Then for some predefined value of a parameter a, all nodes on an empty sphere with a radius greater than ah are considered boundary nodes. In other words: for each point in the domain examine whether it is possible to place a sphere with radius ah such that it contains only that point. If possible, the point is a boundary point and if not, i.e. if the sphere inevitably will contain more than the one point, it is an internal point. While a number of algorithms for recognizing boundaries by means of the alpha-shape method are available, another, and in many ways more straightforward possibility, has been proposed by Cremonesi et al. [23]. The steps in this scheme are as follows. Consider the cloud of particles as shown in Fig. 2(a). A Delaunay triangulation is first performed to generate the convex domain shown in Fig. 2(b). Next, the radius of the circumcircle of each triangle is examined and triangles with a circumcircle radius greater than ah are deleted. For the example at hand, this leads to the final configuration shown in Fig. 2(c). As shown by Cremonesi et al., this procedure is equivalent to the original alpha-shape approach owing to the property that the circumcircles of all triangles generated by the Delaunay triangulation are empty. 2.1.1. Choice of a It is clear that the alpha-shape method involves an element of subjectivity. Indeed, the resulting domain is a direct function of the value of the parameter a. This is illustrated in Fig. 3. A value of a = 1.2 here produces a set of boundaries that in many cases would be deemed reasonable. In fact, any value of a in the interval 0.9 6 a 6 1.3 produces this set of boundaries. Decreasing a below 0.9 leads first to internal voids (a = 0.5) and subsequently to a disintegration of the external boundaries as well (a = 0.4). On the other hand, increasing a above 1.3 leads first to a coalescence of the two distinct sets of points (a = 1.5) and then, for large values of a, to a solid defined by the convex hull inscribing the cloud of points. The conclusion of this example, that a value of a slightly greater than 1 is appropriate, is consistent with experience from application to actual physical problems. For example, for a wide range of coupled fluid–solid interaction problems, Onate et al. [24] conclude that a value of a in the range of 1.3–1.5 is appropriate. For
Fig. 2. Boundary recognition via the scheme of Cremonesi et al. [23]: cloud of points (a), Delaunay triangulation (b), after deletion of triangles with circumcircle greater than ah (c).
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
135
3.2. Constitutive equations The constitutive relations assumed in the following are those of rigid-perfect plasticity. That is, no deformations take place below yield while, once the yield condition is satisfied, plastic deformation occur. In the static case, these are in principle unlimited while in the dynamic case, the extent of plastic straining is limited by the momentum balance equations. Moreover, the flow rule is assumed nonassociated. The constitutive conditions are thus given by
FðrÞ 6 0 e_ ¼ k_ rr GðrÞ _ rÞ ¼ 0; kFð
ð3Þ k_ P 0
where F is the yield function, G is the plastic potential, and k_ is a plastic multiplier. In this paper, both the yield conditions and plastic potential are given in terms of the Mohr–Coulomb relations which under plane strain conditions can be expressed as
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðr11 r33 Þ2 þ 4r213 þ ðr11 þ r33 Þ sin / 2c cos / qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi G ¼ ðr11 r33 Þ2 þ 4r213 þ ðr11 þ r33 Þ sin w F¼
Fig. 3. Boundary recognition via the alpha-shape method.
the examples of the present paper which cover static and dynamic deformation processes of purely cohesive and pure frictional solids, we find the appropriate range of a to be approximately 1.4–1.6. Thus, while the ‘optimal’ value of a is problem dependent, the range of possible values is in practice rather limited and requires a relatively minor effort to establish by simple trial and error. Nevertheless, regardless of the choice of a, the total volume, and thereby, for incompressible materials, the total mass, is bound to fluctuate somewhat in the course of a typical simulation. However, as will be demonstrated via examples, this error is rather small. In principle, it could be controlled in various ways, for example by adjusting a in the course of the simulations. At present no general scheme is available for this purpose although a particular scheme, applicable to fluid dynamics recently has been proposed [25].
3. Governing equations In the following, the governing equations of rate-independent elastoplasticity are briefly summarized. 3.1. Momentum balance and strain–displacement relations Assuming infinitesimal deformations, the strain–displacement relations are given by
e ¼ $u
ð1Þ
where u = (u1, u2, u3)T are the displacements, e = (e11, e22, e33, 2e12, 2e23, 2e31)T are the strains, and $ is the usual linear strain–displacement operator. The momentum conservation equations for a body subjected to a time-varying motion are given by
€; $T r þ b ¼ qu T
N r ¼ t;
in V
on S
ð2Þ
where r = (r11, r22, r33, r12, r23, r31)T are the stresses, V is the domain under consideration, S is its boundary, b are body forces, t are tractions, and N ¼ $ðnxT Þ with n = (n1, n2, n3)T being the outward normal to the boundary and x = (x1, x2, x3)T being the spatial coordinate. The material density is given by q and a superposed dot denotes differentiation with respect to time.
ð4Þ
where / is the friction angle, c is the cohesion, and w is the dilation angle (with w = 0 implying incompressibility). While this constitutive model is very simple, it is for many problems involving very large deformations probably quite appropriate. Indeed, the model verifies the conditions at the critical state and with the initial transient phase leading up to this state thus being ignored, it may be thought of as a reasonably accurate first-order model. Besides failing to capture the conditions prior to fully developed plastic flow, the only obvious flaw in the model is the absence of rate effects. While such effects usually can be ignored at the microscopic level, the processes at this level nevertheless manifest themselves in a what, at the macroscopic level, appears to be a rate effect, essentially implying an increase in friction angle with shear strain rate [26,27]. While such an enhancement is entirely possible within the general framework to be presented, it will not be considered in this paper. Finally, we note that the incremental response in each time step is based on infinitesimal strain theory. This may lead to several errors of which the generation of strains as a result of rigid body motions is the most serious one. However, based on standard convergence analysis, it appears that this and related errors are relatively minor for the kind of time steps used in typical simulations. As such, the price to pay for the convenience of being able to operate with usual infinitesimal strain theory appears to be very small. 4. Discretization and solution In this section, the discretization and solution of the governing equations are considered. The temporal discretization is first detailed after which we derive a time-discrete variational principle that lends itself to a straightforward spatial discretization resulting in a nonlinear mathematical program. We finally comment briefly on the solution of this program. 4.1. Time discretization Introducing velocities, v, the first part of the momentum balance equations may be written in terms of two coupled first-order equations:
$T r þ b ¼ qv_
v ¼ u_
ð5Þ
136
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
These equations are then approximated by means of the hmethod to yield:
$T ½h1 rnþ1 þ ð1 h1 Þrn þ b ¼ q h2 v nþ1 þ ð1 h2 Þv n ¼
v nþ1 v n
unþ1 un Dt
Dt
ð6Þ
where subscripts n and n + 1 refer to the known and new, unknown, states respectively and Dt = tn+1 tn is the time step. Rearranging gives the following equation in the stresses and displacements
Du ~ 2 $ rnþ1 þ b~ ¼ q Dt T
h1 h2 1 v 1 h1 T ~ ~ nþ b¼ bþq $ rn h1 Dt h1
ð8Þ
After the displacements at tn+1 have been determined, the velocities are computed from
1 Du ð1 h2 Þv n h2 Dt
ð9Þ
The static boundary conditions are approximated in an analogous manner leading to
N T ½h1 rnþ1 þ ð1 h1 Þrn ¼ t;
on S
ð10Þ
or:
N T rnþ1 ¼ ~t;
on S
ð11Þ
where
~t ¼ 1 t 1 h1 N T rn : h1 h1
ð12Þ
The above time integration scheme is unconditionally stable provided that h1 P 12 and h2 P 12. For h1 ¼ h2 ¼ 12 it coincides with the Newmark average acceleration scheme while, for h1 = 0, a fully explicit and conditionally stable scheme is obtained. In the following, we will only consider implicit and unconditionally stable schemes. Regarding the constitutive model, we require that all relations are satisfied at time tn+1:
Fðrnþ1 Þ 6 0
$ðDuÞ ¼ Dkrr Gðrnþ1 Þ DkFðrnþ1 Þ ¼ 0; Dk P 0:
ð13Þ
4.2. Variational formulation We are now in a position to postulate that the above time discrete governing equations can be stated in terms of the following time-discrete optimization problem, or variational principle,
min Du
max
ðr;rÞnþ1
~ 1 r nþ1 iV þ hr nþ1 ; DuiV DkFðrnþ1 Þ 1=2Dt 2 hr nþ1 ; q
~ Dui h~t; Dui hrnþ1 ; $ðDuÞiV hb; V S
~ 1 rnþ1 iV þ hr nþ1 ; DuiV 1=2Dt2 hr nþ1 ; q subject to Fðrnþ1 Þ 6 0
(
$T rnþ1 þ b~ ¼ rnþ1 ; in V N T rnþ1 ¼ ~t;
on S dJ 2 ~ 1 ¼ Dt q r nþ1 þ Du ¼ 0 dr nþ1 dJ ¼ $u Dkrr Fðrnþ1 Þ ¼ 0: drnþ1
4.3. Spatial discretization The principle Eq. (14) is discretized in space by postulating finite element approximations to the state variables r, r and u. Using standard finite element notation, we have
rðxÞ N r ðxÞr^ rðxÞ N r ðxÞr^ ^ ; $u Bu ðxÞu ^ uðxÞ N u ðxÞu
hx; yiA ¼
xT ydA
ð18Þ
^ are the nodal variables, the N matrices contain the ^; ^ r and u where r shape functions, and Bu ¼ $N u . In standard finite formulations, the stress shape functions in Nr are chosen as being continuous within the elements and discontinuous between elements while functions in Nu in addition are continuous between elements and usually one polynomial degree higher than those of Nr. The shape functions for the dynamics forces are taken to be identical to those of the displacements but are discontinuous between elements. Substituting the above approximations into the variational principle Eq. (14) leads to the following discrete principle
max
^ ;^rÞnþ1 ðr
^ Du
subject to
r^ Tnþ1 BDu^ f T Du^ 1=2^rTnþ1 D^rnþ1 þ ^rTnþ1 ADu^ ^ jnþ1 Þ 6 0; Fðr
j ¼ 1; . . . ; nr ð19Þ
where we have used the notation
Z
ð17Þ
Furthermore, by addition of appropriate penalty terms to the functional, the yield and plastic consistency conditions are recovered (see [17,28,7,18] for details). From the second set of conditions, we see that the variables rn+1 are the dynamic forces ~ Du=Dt 2 . The momentum balance equations are thus verir nþ1 ¼ q fied. However, from the last set of conditions we see that the flow rule is associated, i.e. the flow potential is equal to the yield potential F. This is a consequence of the variational formulation and is not possible to circumvent directly. In this paper we use the approach suggested in [7] of replacing the original yield function by an ‘effective’ one which, when used as a flow potential, yields the correct plastic strains. This approach was shown to produce reliable results for a range of static small deformation problems and in this paper we demonstrate that it is equally reliable for dynamic problems involving very large deformations.
min ð14Þ
ð16Þ
where Dk are Lagrange multipliers. Next, taking functional derivatives, the following equations appear:
dJ ¼ dDu
q
v nþ1 ¼
~ Dui h~t; Dui J ¼ hrnþ1 ; $ðDuÞiV hb; V S
ð7Þ
where Du = un+1 un and
q~ ¼
4.2.1. Euler–Lagrange equations The equivalence between the above optimization problem and the governing equations follows by showing that the Euler–Lagrange equations associated with Eq. (14) do indeed reproduce the governing equations. To this end, define the functional:
ð15Þ
where the yield function is imposed at a finite number of points nr and
A
and r is a set of variables whose significance will become apparent shortly.
B¼
Z V
N Tr Bu dV
ð20Þ
137
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
f ¼
Z V
D¼
Z V
A¼
Z V
~ þ N Tu bdV
Z S
N Tu ~tdS
ð21Þ
~ 1 N r dV N Tr q
ð22Þ
N Tu N r dV:
ð23Þ
The minimization part of Eq. (19) may be solved first to yield
^ nþ1 þ AT ^r nþ1 f ¼ 0 BT r
ð24Þ
which is included as a constraint in the remaining maximization part of the problem:
maximize ^ ;^r Þnþ1 ðr
1=2^rTnþ1 D^r nþ1
^ nþ1 þ AT ^rnþ1 ¼ f subject to BT r ^ jnþ1 Þ 6 0; Fðr
ð25Þ
j ¼ 1; . . . ; nr :
This is the problem solved in each time step. In the course of solving the problem, the displacement increments appear as the dual variables, or Lagrange multipliers, associated with the discrete equilibrium constraints. The above problem is similar in structure to those arising from limit analysis [13,14], static elastoplasticity [17,7,18], and granular contact dynamics [19–21].
Fig. 4. Contact specification.
and those of the rigid body. The edges in this triangulation then define the potential contacts. Considering each individual particle as a free particle at a distance g0 from the boundary and with potential contact forces p and q, a variational formulation of the corresponding frictional contact problem is given by [30,20,21]:
maximize
g0 p
subject to jqj lp 6 0; p P 0
where the standard Coulomb friction law has been used, l being the friction coefficient. This principle implies the following kinematics:
DuN ¼ g 0 k
ð27Þ
DuT ¼ k 4.4. Solution algorithm
ð26Þ
The optimization problems Eq. (25) generated by the scheme in each time step are solved by means of the second-order cone programming (SOCP) solver MOSEK [29]. Such solvers have become increasingly more powerful and robust in recent years and require a computational effort similar to conventional implicit Newton– Raphson based schemes. The SOCP standard form operates with a linear objective function [as opposed to the quadratic one in Eq. (25)] and with inequality constraints that are slightly different from those of the Mohr–Coulomb criterion. However, a transformation of Eq. (25) into SOCP standard form is straightforward as has been documented previously, e.g. [14,18].
where k P 0 is a Lagrange multiplier such that k(j q j lp) = 0 and DuN and DuT are the normal and tangential displacement increments respectively. This means that the sliding rule is associated in the sense that sliding is accompanied by a dilation of relative magnitude l. As in the continuum case, this dilation is not desirable and can rarely be justified physically. However, as discussed in detail in [30], the dilation arising at frictional contacts may be viewed as an artifact of the time discretization that gradually decreases as the time step is reduced. Considering all potential contacts and incorporating the above contact problem into the original problem Eq. (25) leads to a final problem of the type
5. Treatment of contact
maximize 1=2^rTnþ1 D^r nþ1
For general problems, proper treatment of frictional contact is an important component of the overall problem. In this paper, we will consider the contact between the deforming solid and rigid boundaries. Extension to general non-stationary rigid bodies is straightforward while the extension to contact with other deformable bodies, including self-contact, is less so although it in principle can be dealt with using the methodology described in the following. The treatment of the above mentioned types of contact can be accommodated within the general variational framework by using the recently proposed granular contact dynamics method [30,20,21]. Though originally conceived for particle dynamics of the discrete element type, the basic principles are equally applicable in a finite element context. The procedure is as follows. Consider a body discretized by finite elements as shown in Fig. 4. Initially, each boundary node of the system is viewed as a particle with zero radius. Next, the particles that are likely to come into contact with the rigid body (in Fig. 4 taken to be a single segment) are identified. This can be done in various ways as described in detail in [21]. The most general approach, and the one adopted in this paper, is to discretize the boundary of the rigid body and perform a Delaunay triangulation with respect to the boundary nodes of the deformable solid
^ ;^r Þnþ1 ðr
nc X g 0;j pj j¼1
^ nþ1 þ AT ^r nþ1 þ ET q ¼ f subject to B r T
^ jnþ1 Þ 6 0; Fðr pj ¼
nTj
j ¼ 1; . . . ; nr
ð28Þ
qj ; j ¼ 1; . . . nC
^ Tj qj ; qj ¼ n jqj j lpj 6 0;
j ¼ 1; . . . nC j ¼ 1; . . . nC T
with q = (q1, q2) being the nodal forces, n the normal of the rigid ^ ¼ ðn2 ; n1 ÞT . Furthermore, E is an index matrix of boundary and n zeros and ones and nC is the number of potential contacts. 6. Examples In the following, four examples demonstrating the capabilities of the new scheme are presented. For all examples, a reasonable value of a was first determined by trial and error for a coarse discretization (in both time and space) and for a very limited range of the parameter set involved in each problem. As discussed previously, this procedure is less than ideal but at present the only feasible one. Also, in practice, it is not particularly arduous and is usually completed in a fraction of the time that it takes to solve a given problem (using more refined temporal and spatial discretizations) with a fixed value of a. By this initial calibration, a value
138
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
of a = 1.4 was adopted for the first two examples, while a = 1.6 for the two last. In all examples, the gravitational acceleration is g = 9.8 m/s2. 6.1. Cylinder–soil interaction The first example concerns the interaction between a rigid weightless cylinder and a Tresca soil. The cylinder first penetrates vertically to a depth of 0.8 m after which it is displaced horizontally for a total distance of 2 m (see Fig. 5). This type of problem involving the combined horizontal and vertical motion of a rigid cylindrical object resembles that of a pipeline subjected to either external loads due to waves and current or to internal temperature and pressure loadings. The analysis conducted uses 280 displacement increments of magnitude equal to 0.01 m (80 for the penetration and 200 for the subsequent horizontal motion). The cylinder-soil interface is assumed purely frictional with a friction coefficient l = tan 30° ’ 0.577. Other parameters are as follows: soil density q = 1.6 g/cm3, undrained shear strength c = cu = 1 kPa, and cylinder diameter D = 1 m. The analysis is fully dynamic with a time step of 1s. Under these conditions, inertial effects are negligible and the simulation may be regarded as quasi-static. The results of the simulation in terms of load–displacement curves are shown in Fig. 6. During the initial penetration phase, the resultant horizontal force remains at zero while the vertical force increases approximately linearly up to around Fy = 18 kN/m. The cylinder is then displaced horizontally while being fixed in its vertical position of uy = 0.8 m. This change in motion brings about an immediate increase in horizontal force and a decrease in the vertical force. This is followed by a more gradual change of the forces until the motion becomes steady at around ux = 1.5 m. At this point, the resultant forces are Fx ’ 9 kN/m and Fy = 7.5 kN/ m respectively. Thus, an appreciable vertical force is required to prevent the cylinder moving towards the surface in response to the purely horizontal motion. These findings echo those of [31] where a vertically constrained cylinder was plowed through a granular material resulting in a similar ratio of vertical to horizontal force. Finally, we note that the fundamental change in geometry associated with the self-contact at the back of the cylinder between uH = 1.0 m and uH = 1.5 m is handled seamlessly and without resorting to any additional rules or procedures. 6.2. Accretionary wedge The next problem concerns the deformation of an accretionary wedge as shown in Fig. 7. This is a common model problem in structural geology used to study, among other things, mountain building via folding and thrusting events [32–35]. The wedge has a height of 8 cm and a length of 140 cm and is composed of a purely frictional material with / = 30° and w = 0 and a density of 2 g/cm3. The left vertical and bottom horizontal boundaries are rigid with friction coefficients of l = 30° and l = 15° respectively. The simulation is carried out by displacing the left wall horizontally with a velocity of 0.1 cm/s. A total of 1000 steps, with a constant time step of 0.16 s, were used to yield a total wall displacement of 16 cm. Snapshots of the deformations at different times are shown in Fig. 7. Contrary to what might be expected, the deformations do not proceed by evolution of one of more dominant shear bands. Rather, the location of the shear bands changes in a somewhat unpredictable manner: the classic V-band at the initial phases first propagates in the direction of loading (Du = 8.0 cm). Further loading then gives way to a new shear band much like the initial one (Du = 10 cm). This again gives way to another family of shear
Fig. 5. Interaction of a rigid cylinder with a Tresca soil.
bands further down the direction of loading. The end result (when viewed as a continuous process) is a seemingly random appearance of shear bands. This phenomenon has recently been studied in some detail by Mary et al. [36] using a semi-analytical slip line type model. The conclusion is here that the system displays deterministic chaos of a nature similar to more common chaotic dynamical system. Without repeating this analysis in the present context,
Normalized resultant force, F (kN/m)
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
139
20 15 FV
10 FH
5
Fig. 8. Granular column problem: initial column and final deposit.
0
0
−5 0
0.4
0.5
1.0
1.5
2.0 uH
0.8 uV
Displacment, u (m) Fig. 6. Horizontal and vertical resultant forces acting on the cylinder.
we do note that the behavior of the system displays a behavior consistent with that observed and analyzed by Mary et al.
6.3. Collapse of granular columns The next example concerns the collapse of granular columns. The setup is shown in Fig. 8. A column of granular material of width d0 and height h0 is allowed to collapse under the action of gravity by quickly removing the right wall. This type of experiment, which may be seen as a model landslide, has received significant attention in recent years, both from an experimental point of view [37–39] as well as in terms of reproducing the experiments using both particle methods [40–45,30,21] and continuum models [45–47]. The two-dimensional plane strain version of the experiment considered in the following is due to Lube et al. [38]. On the basis of a rather extensive experimental program, the following fits to
the height, h1, and diameter, d1, of the final deposit were determined:
h1 ¼ d0
a;
a 6 1:15 2
1:1a5 ; a > 1:15
8 a < 1:8 > 1:6a; d1 d0 < ¼ transition region; 1:8 6 a 6 2:8 > d0 2 : 2:2a3 ; a > 2:8
ð29Þ
ð30Þ
where a = h0/d0 is the initial aspect ratio. We note that there is a transition region at 1.8 6 a 6 2.8 for which no expressions for the final width were given. To validate the proposed scheme, a series of experiments involving columns with different aspect ratios were conducted. In all simulations, the following parameters were used: density q = 2 g/cm3, cohesion c = 0, friction angle / = 31° (as quoted by Lube et al.) or / = 35°, and dilation angle w = 0. The proceed from time psimulations ffiffiffiffiffiffiffiffiffiffi t = 0 and are terminated at t ¼ t= h0 =g ¼ 4:0. Fig. 9 shows three columns with aspect ratios a = 0.5, 1.0, and 7.0 at different times instants of the collapse process using a time step of Dt ¼ 0:01. The internal deformation patterns observed, in particular the apparent ‘buckling’ of internal columns, corresponds very well with what can be observed experimentally and in DEM type simulations [30]. As expected, and predicted from Eq. (29), the two shortest
Fig. 7. Accretionary wedge problem. Colors are proportional to the equivalent shear strain. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
140
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
100 time steps
200 time steps
400 time steps
Fig. 9. Collapse of granular columns (dots represent the finite element nodes).
3.0
Normalized final height, h∞ /d0
the model of Andrade et al. [27], better fits to the experiments can be expected. To gauge the influence of the time step, the a = 1 simulation was run with different time steps. The resulting final deposits are shown in Fig. 11. As seen, the shape and dimensions of the final deposit are relatively insensitive to the time step. Indeed, the final 2
Volume change (%)
columns collapse only partially while the spread of the tallest one is much wider. A more complete comparison with the experimental fits of Lube et al. is shown in Fig. 10. We here see that the columns essentially lack resistance, even for a friction angle of 35°. That is, the final height is underpredicted while the final run-out distance is overpredicted, and with the deviation between experiment and simulation increasing with increasing initial aspect ratio. As discussed in Section 3.2, the key feature missing in the constitutive model adopted is a rate-dependence of the friction angle on the shear strain rate. With such an enhancement, for example following
Fig. 11. Convergence of final deposit as function of time step for a = 1 (dots represent the finite element nodes).
2.5
0
−2 100 time steps
−4
200 time steps 400 time steps
−6
2.0
0
1
2
3 −
Normalized time, t
1.5
Fig. 12. Volume change as function of time for a = 1.
1.0
φ = 31°
0.5 0
φ = 35° 0
2
4
6
8
10
12
Normalized final width, (d∞ -d0 )/d0
Aspect ratio, a = h0 /d0 12 100 cm
10 8 6 4
φ = 31°
2 0
35 cm
φ = 35° 0
2
4
6
8
10
12
8 cm
Aspect ratio, a = h0 /d0 50 cm Fig. 10. Normalized final height and width of granular columns as function of aspect ratio (dots represent the finite element nodes).
Fig. 13. Silo discharge problem.
4
141
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
0.0 s
3.5 s
1.0 s
4.0 s
0.0 s
3.0 s
6.0 s
7.0 s
5.0 s
2.0 s
4.5 s
8.0 s
Fig. 15. Discharge of rough walled silo (dots represent the finite element nodes).
Fig. 14. Discharge of smoothed walled silo (dots represent the finite element nodes).
7
6.4. Discharge of silo The final example concerns the discharge of a two-dimensional silo with the dimensions shown in Fig. 13. The material parameters are: density q = 1.5 g/cm3, cohesion c = 0, friction angle / = 31°, and dilation angle w = 0°. The solid mass confined by the silo wall is initially at rest after which, at time t = 0, the bottom wall is removed. Two different sets of simulations were run: assuming perfectly smooth walls and assuming perfectly rough walls (corresponding to a wall-solid friction coefficient of l = tan /). Snapshots of the discharge in the two cases are shown in Figs. 14 and 15. As can be seen, the wall-solid interface properties have a significant influence on the overall flow pattern. Again, these patterns are consistent with those observed in both experiments and particle simulations [48]. Moreover, the wall-solid interface conditions influences the total discharge time, from approximately 5 s for the smooth silo to almost twice that for the rough silo.
Δ t = 0.04s
Volume in silo (1000 cm3 )
height changes very little between the three different time steps while the run-out distance does increase somewhat between the first two simulations involving 100 and 200 time steps respectively. A similar trend is observed for other aspect ratios. Finally, we investigate the mass balance properties of the scheme as function of the time step. The results, shown in Fig. 12, reveal the loss of mass to be around 2–6% and decreasing as the time step is decreased (the number of time steps increased).
6
Δ t = 0.02s
5
Δ t = 0.01s Δ t = 0.005s
4 3 2 1 0
0
1
2
3
4
5
6
7
Time (s) Fig. 16. Silo volume versus time for different time steps.
The convergence of the results in terms the discharge versus time for different time steps is shown in Fig. 16 for the smooth walled silo. As in the previous example, some 500 time steps (corresponding to Dt = 0.01) appears to be sufficient to ensure a reasonable degree of accuracy. This trend is the same for the rough walled silo. 7. Conclusions The variant of the Particle Element Method has been developed and applied to a range of geomechanics problems involving very large deformations. A variational approach has been taken whereby
142
X. Zhang et al. / Computers and Geotechnics 54 (2013) 133–142
the discrete governing equations give rise to a second-order cone program that can be solved in an efficient and robust manner using available tools. Moreover, frictional contact is accounted for in a straightforward manner and does not give rise to any complications. A number of examples have been presented and while the results generally are encouraging, there are several issues that warrant further consideration. Firstly, regarding the alpha-shape method which is at the heart of the PFEM. As discussed, the parameter a is at present chosen on the basis of simple trial and error. Clearly, a more rigorous or at least automated procedure would be highly desirable. This could well involve a parameter that would vary throughout the simulation in contrast to the constant value currently used for the entire simulation. Secondly, while the simple rate-independent model used in this paper in many cases probably is quite adequate, it has its limitations. This is most apparent in the granular column collapse simulations where the deviation between simulation and experiment increases as the column height, and thereby the characteristic velocity, increases. A rate-dependent model, essentially implying an increase in friction angle with velocity, is here likely to improve the results.
References [1] Nazem M, Sheng DC, Carter JP, Sloan SW. Arbitrary Lagrangian–Eulerian method for large-strain consolidation problems. Int J Numer Anal Meth Geomech 2008;32:1023–50. [2] Nazem M, Carter JP, Airey D. Arbitrary Lagrangian–Eulerian method for dynamic analysis of geotechnical problems. Comput Geotech 2009;36:549–57. [3] Belytschko T, Liu WK, Moran B. Nonlinear finite elements for continua and structures. Wiley; 2000. [4] Donea J, Huerta A, Ponthot JP, Rodriguez-Ferran A. Arbitrary Lagrangian Eulerian method. In: Stein E, de Borst R, Hughes TJR, editors. Encyclopedia of Computational Mechanics, vol. 1; 2004. [5] Hu Y, Randolph MF. A practical numerical approach for large deformation problems in soil. Int J Numer Anal Meth Geomech 1998;22:327–50. [6] Yu L, Hu YX, Liu J, Randolph MF, Kong XJ. Numerical study of spudcan penetration in loose sand overlying clay. Comput Geotech 2012;46:1–12. [7] Krabbenhoft K, Karim MR, Lyamin AV, Sloan SW. Associated computational plasticity schemes for nonassociated frictional materials. Int J Numer Meth Eng 2012;89:1089–117. [8] Idelsohn SR, Onate E, Del Pin F. The particle finite element method: a powerful tool to solve incompressible flows with free-surfaces and breaking waves. Int J Numer Meth Eng 2004;61:964–89. [9] Onate E, Idelsohn SR, Del Pin F, Aubry R. The particle finite element method an overview. Int J Comput Meth 2004;1:267–307. [10] Onate E, Idelsohn SR, Celigueta MA, Rossi R. Advances int the particle finite element method for the analysis of fluid-multibody interaction and bed erosion in free surface flows. Comput Meth Appl Mech Eng 2008;197:1777–800. [11] Lyamin AV, Sloan SW. Lower bound limit analysis using non-linear programming. Int J Numer Meth Eng 2002;55:573–611. [12] Lyamin AV, Sloan SW. Upper bound limit analysis using linear finite elements and non-linear programming. Int J Numer Anal Meth Geomech 2002;26:181–216. [13] Krabbenhoft K, Damkilde L. A general nonlinear optimization algorithm for lower bound limit analysis. Int J Numer Meth Eng 2003;56:165–84. [14] Krabbenhoft K, Lyamin AV, Sloan SW. Formulation and solution of some plasticity problems as conic programs. Int J Solids Struct 2007;44:1533–49. [15] Vicente da Silva M, Antao AN. A nonlinear programming method approach for upper bound limit analysis. Int J Numer Meth Eng 2007;72:1192–218. [16] Vicente da Silva M, Antao AN. Upper bound limit analysis with a parallel mixed finite element formulation. Int J Solids Struct 2008;45:5788–804. [17] Krabbenhoft K, Lyamin AV, Sloan SW, Wriggers P. An interior-point method for elastoplasticity. Int J Numer Meth Eng 2007;69:592–626.
[18] Krabbenhoft K, Lyamin AV. Computational Cam clay plasticity using secondorder cone programming. Comput Meth Appl Mech Eng 2012;209– 212:239–49. [19] Krabbenhoft K, Lyamin AV, Huang J, Vicente da Silva M. Granular contact dynamics using mathematical programming methods. Comput Geotech 2012;43:165–76. [20] Krabbenhoft K, Huang J, Vicente da Silva M, Lyamin AV. Granular contact dynamics with particle elasticity. Granul Matter 2012;14:607–19. [21] Huang J, Vicente da Silva M, Krabbenhoft K. Three-dimensional granular contact dynamics with rolling resistance. Comput Geotech 2013;49:289–98. [22] Edelsbrunner H, Mucke EP. Three dimensional alpha shapes. ACM Trans Graph 1994;13(1):43–72. [23] Cremonesi M, Frangi A, Perego U. A lagrangian finite element approach for the analysis of fluid-structure interaction problems. Int J Numer Meth Eng 2010;84:610–30. [24] Onate E, Idelsohn SR, Celigueta MA, Rossi R, Marti J, Carbonell JM, et al. Advances in the particle finite element method (PFEM) for solving coupled problems in engineering. In: Particle-based methods fundamentals and applications computational methods in applied sciences, vol. 25. Springer; 2011. p. 47–54. [25] Ryzhakov P, Onate E, Rossi R, Idelsohn S. Improving mass conservation in simulation of incompressible flows. Int J Numer Meth Eng 2012;90:1435–51. [26] Forterre Y, Jop P, Pouliquen O. A constitutive law for dense granular flows. Nature 2006;441:727–30. [27] Andrade JE, Chen Q, Lee PH, Avila CF, Evans TM. A constitutive law for dense granular flows. J Mech Phys Solids 2012;60:1122–36. [28] Krabbenhoft K. A variational principle of elastoplasticity and its application to the modeling of frictional materials. Int J Solids Struct 2009;46:464–79. [29] Andersen ED, Roos C, Terlaky T. On implementing a primal–dual interior-point method for conic quadratic optimization. Math Program 2003;95:249–77. [30] Krabbenhoft K, Lyamin AV, Huang J, Vicente da Silva M. Granular contact dynamics using mathematical programming methods. Comput Geotech 2012;43:165–76. [31] Ding Y, Gravish N, Goldman DI. Drag induced lift in granular media. Phys Rev Lett 2011;106:028001. [32] Souloumiac P, Leroy YM, Maillot B, Krabbenhoft K. Predicting stress distributions in fold-and-thrust belts and accretionary wedges by optimization. J Geophys Res 2009;114:B09404. [33] Souloumiac P, Krabbenhoft K, Leroy YM, Maillot B. Failure in accretionary wedges with the maximum strength theorem: numerical algorithm and 2d validation. Comput Geosci 2010;14:793–811. [34] Cubas N. Prediction of thrusting sequences in accretionary wedges. J Geophys Res 2008;113:B12412. [35] Maillot B, Leroy Y. Kink-fold onset and development based on the maximum strength theorem. J Mech Phys Solids 2006;54:2030–59. [36] Mary BCL, Maillot B, Leroy YM. Deterministic chaos in frictional wedges revealed by convergence analysis. Int J Numer Anal Meth Geomech 2013. http://dx.doi.org/10.1002/nag.2177 [in press]. [37] Lube G, Huppert HE, Sparks RSJ, Hallworth MA. Axisymmetric collapses of granular columns. J Fluid Mech 2004;508:175–99. [38] Lube G, Huppert HE, Sparks RSJ, Freundt A. Collapses of two-dimensional granular columns. Phys Rev E 2005;72:041301. [39] Lajeunesse E, Mangeney-Castelneau A, Vilotte JP. Spreading of a granular mass on an horizontal plane. Phys Fluids 2004;16:2731–8. [40] Lacaze L, Phillips JC, Kerswell RR. Planar collapse of a granular column: experiments and discrete element simulations. Phys Fluids 2008;20:063302. [41] Lacaze L, Kerswell RR. Axisymmetric granular collapse: a transient 3D flow test of viscoplasticity. Phys Rev Lett 2009;102:108305. [42] Zenit R. Computer simulations of the collapse of a granular column. Phys Fluids 2005;17:031703. [43] Staron L, Hinch EJ. Study of the collapse of granular columns using twodimensional discrete-grain simulation. J Fluid Mech 2005;545:1–27. [44] Staron L, Hinch EJ. The spreading of a granular mass: role of grain properties and initial conditions. Granul Matter 2007;9:205–17. [45] Lagree PY, Staron L, Popinet S. The granular column collapse as a continuum: validity of a two-dimensional Navier–Stokes model with a l(i)-rheology. J Fluid Mech 2011;686:378–408. [46] Mangeney-Castelnau A, Bouchut F, Vilotte JP, Lajeunesse E, Aubertin A, Pirulli M. On the use of saint venant equations to simulate the spreading of a granular mass. J Geophys Res 2005;110:B09103. [47] Kerswell RR. Dam break with Coulomb friction: a model for granular slumping? Phys Fluids 2005;17:057101. [48] Gonzalez-Montellano C, Ayuga F, Ooi JY. Discrete element modelling of grain flow in a planar silo: influence of simulation parameters. Granul Matter 2011;13:149–58.