# Application of the Bipotential Theory to a Nonassociated Drucker–Prager Model

Authors: Yang-Jing Zhou; Zhi-Qiang Feng; Juan Antonio Rojas Quintero; Jian-Bo Zhou; Hui-Jian Chen
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

The bipotential theory allows us to describe nonassociated material laws. In this paper, we propose its application to the Drucker–Prager model. With a new description of the implicit flow rules, we propose dual constitutive cones as well as five forms of the bipotential function: the elastic stage in rate form, the plastic stage in rate form, the elastic stage in incremental form, the plastic stage in incremental form, and the elastoplastic stage in incremental form. By combining these with the finite element method, a numerical strategy that deals with the nonassociated Drucker–Prager model is obtained. Two examples are simulated to verify the accuracy, the stability, and the practicability of the algorithm in civil engineering.


---

## Body

## 1. Introduction

The exploration of the material constitutive laws plays an important role in its numerical modeling. Much work has been done in dealing with associated plastic materials such as metals. A return mapping algorithm was proposed by Simo and Taylor to simulate the behavior of associated elastoplastic materials [1]. On the basis of the return mapping algorithm, many other algorithms were proposed to study different associated constitutive laws of materials [2, 3, 4]. However, as we know, most rock-soil materials belong to the nonassociated materials category, and their behavior is more complex [5] because the plastic flow direction is not orthogonal to the yield surface. Typical methods dealing with nonassociated rock-soils consist in generating two potential functions: a yield function which delimits the zone of plasticity, and the flow potential which defines the plastic flow direction. The use of these two potentials to describe nonassociated rock-soils breaks the frame of orthogonality between stress and plastic strain, which is a classical property in solid mechanics. However, it is desirable to obtain a unique potential which can express both the yield and the flow rules. The bipotential method put forward by De Saxcé and Feng [6] achieves this goal. More recently, Sun et al. proposed the fractional plasticity which brings some new ideas to the study of nonassociated plasticity [7, 8].On the basis of the Legendre transformations, the Fenchel inequality with respect to two conjugate variables was established [9]. According to whether the free energy could be separated into two convex potential functions or not, materials were divided into the generalized standard materials (GSMs) and implicit standard materials (ISMs) [10]. De Saxcé et al. have shown that the bipotential algorithm could not only satisfy the constitutive laws of ISMs but also meet the implicit flow requirements of material conjugated variables [11, 12]. Nowadays, the bipotential method can be used to solve constitutive laws of ISMs, including frictional contact, nonassociated metal materials, or even nonassociated rock-soils among other examples.Feng proposed a bipotential method combining both the Signorini and Coulomb conditions to solve the frictional contact [13]. De Saxcé and Feng summarized the main process of the bipotential algorithm and proved that the algorithm had a high efficiency in dealing with contact issues [14]. The Uzawa algorithm was used to replace the Newton–Raphson algorithm to solve the nonlinear equations. This improved the speed of the solving strategy [15]. Others have also developed the algorithm in frictional contact considering large deformations [16], anisotropy [17], and impact dynamics [18].The bipotential method is also applied to solve the constitutive laws of materials. On the basis of De Saxcé’s work, Bodovillé succeeded in applying the bipotential method in kinematic hardening models [19]. Bouby et al. conducted a shakedown analysis, by comparing different models with the linear unlimited, the linear limited, and the nonlinear kinematic hardening [20]. Chaaba et al. stated the superiority of combining the bipotential method with sequential limit analysis in disposing the nonassociated nonlinear hardening models [21]. Zhou et al. compared the accuracy and stability of the return mapping algorithm against the bipotential algorithm in two typical nonassociated models [22]. Cheng et al. simulated the behaviors of a hollow sphere modeled with a nonassociated material under an isotropic loading [23]. Aiming at ductile porous materials, the bipotential method was applied into limit analysis and a stress-based variational model was developed. In addition, some homogenization problems were also studied. Bousshine et al. derived the bipotential function of nonassociated rock-soils in detail [24]. Hjiaj et al. discussed the singular point on the nonassociated model, keeping the integrity of the bipotential algorithm [25]. Berga aimed at the nonassociated plasticity of rock-soils, giving a detailed process in establishing its bipotential function [26]. Furthermore, he proposed a numerical process and simulated some basic examples [27]. On the basis of Berga’s work, Zhou et al. applied the bipotential algorithm in the geotechnical context [28].The abovementioned investigations show some advantages provided by the bipotential method in dealing with ISMs. However, regarding the expression of the dual cone of the Drucker–Prager model, further detail is needed to explain how the new transmission meets the implicit flow requirements. This paper focuses on the implicit flow laws of the perfect elastoplastic Drucker–Prager model, detailing the derivation process of the incremental elastoplastic bipotential of the Drucker–Prager model. When applying the bipotential algorithm into the finite element method, a large time increment is used to conduct the numerical implementation of the nonassociated Drucker–Prager model of the whole structure.This paper is structured as follows. The cones of the nonassociated Drucker–Prager model are introduced in detail in Section2. Combining them with the bipotential theory, in Section 3, we establish the incremental bipotential function of the nonassociated Drucker–Prager model by a deriving strategy. Two numerical examples are proposed in Section 4, verifying the accuracy, the stability, and the practicability of the bipotential Drucker–Prager model implementation. Finally, Section 5 gives some conclusions of the presented work.


## 2. Nonassociated Drucker–Prager Model

The nonassociated Drucker–Prager model is a classic three-parameter model: the cohesionc, the internal friction angle ϕ, and the dilatancy angle θ [29]. The relations of stress and plastic strain rate can be expressed by a Drucker–Prager constitutive cone, which consists of the hydrostatic part and the deviatoric part. They are σ=sm,s and ε˙p=e˙mp,e˙p, where sm=Trσ/3, s=σ−sm1, e˙mp=Trε˙p, and e˙p=ε˙p−e˙mp1/3, where 1 is the second-order unit tensor. In the elastic state, the stress variable σ=sm,s belongs to the constitutive stress cones Kσ, as shown in Figure 1(a). The expression of these stress cones is as follows.(1)Kσ=sm,sskd+smtanϕ≤c,where ⋅ stands for the Euclidean norm and the kd parameter is a function of the internal friction angle ϕ as kd=32/9+12tan2ϕ. In the principal stress space, the Drucker–Prager cone is the inscribed compressive meridian of the Mohr–Coulomb pyramid. Due to the nonassociated flow rules of the Drucker–Prager model, the dual plastic rate cone Kε is obtained as shown in Figure 1(a). The regular and the apex points on the stress cone have different flow expressions. At a regular point, namely, sm≠c/tanϕands≠0, Kσ is differentiable, and the flow direction Kεr can be expressed as(2)Kεr∈e˙mp,e˙pe˙mp=kdtanθe˙p.Figure 1
Cones of the nonassociated Drucker–Prager model. (a) Constitutive cones. (b) Dual cones.
(a)
(b)At the apex point, namely,sm=c/tanϕands=0, Kσ is nondifferentiable, and the flow direction cannot be determined. A convex set is used here. There exists a cone Kεa such that(3)Kεa∈e˙mp,e˙pe˙mp>kdtanθe˙p.So, associating (2) with (3), the plastic rate cone Kε is unified as follows:(4)Kε∈e˙mp,e˙pe˙mp≥kdtanθe˙p.The implicit flow laws of the nonassociated Drucker–Prager model can be written as follows:(5)σ∈∂ΨKεε˙p,ε˙p∈∂ΨKσσ,where ΨKσσ and ΨKεε˙p are indicator functions.Equation (5) can also be written with the hydrostatic part and the deviatoric part:(6)sm∈∂e˙mpΨKεε˙ps∈∂e˙pΨKεε˙p,e˙mp∈∂smΨKσσe˙p∈∂sΨKσσ.From Figure1(a), it can be noted that when θ=ϕ, Kε and Kσ are dual cones, the material is associated, and the flow direction is normal to its yield surface. If θ≠ϕ, the material is nonassociated, and there exists an angle ϕ−θ between the plastic flow direction and the yield surface, as shown in Figure 1(a). So, the implicit flow laws between stress and plastic strain are therefore not met. There is a new expression which not only satisfies the nonassociated flow rules but also meets the constitutive requirements. Moving the stress cone along the negative sm direction by a distance c/tanϕ, the cone translates from the location sm,s to 0,0, as shown in Figure 1(b). A new stress cone Kσ∗ can be defined as(7)Kσ∗=sm,sskd+smtanϕ≤0,with Kε∗ being the dual cone of Kσ∗ and is defined as(8)Kε∗=e˙mp,e˙pe˙mp≥kdtanϕe˙p.According to the geometric features ofKεε˙p, Kε∗ε˙p, Kσσ, and Kσ∗σ and combining them with (6), while setting cϕ=c/tanϕ, the implicit flow laws of the nonassociated Drucker–Prager model under new constitutive cones can be stated as(9)σ∈∂ΨKεε˙p=∂ΨKε∗ε˙p+cϕ,s,ε˙p∈∂ΨKσσ=∂ΨKσ∗σ+kdtanθ−tanϕe˙p,s.Equation (9) is divided into a hydrostatic part and a deviatoric part:(10)sm∈∂e˙mpΨKε∗ε˙p+cϕs∈∂e˙pΨKε∗ε˙p,e˙mp∈∂smΨKσ∗σ+kdtanθ−tanϕe˙pe˙p∈∂sΨKσ∗σ.Equation (10) is the new expression of the implicit flow laws based on the dual cone of the Drucker–Prager model. On the basis of the bipotential theory, we try to find a bipotential function bσ,ε˙p, which satisfies the constitutive laws of the nonassociated Drucker–Prager model and meets the requirements of implicit flow laws.


## 3. Bipotential of the Nonassociated Drucker–Prager Model

In solid mechanics, many materials can be represented by a relationship of a set of coupled variablesσ,ε˙. According to their free energy, materials are divided into generalized standard materials (GSMs) and implicit standard materials (ISMs) [12], where GSMs are similar to associated materials, whose plastic flow directions are normal to their yield surfaces. ISMs have the same characteristics of nonassociated materials, whose plastic flow directions and yield surfaces are not orthogonal. Under the traditional plastic mechanics, the expressions of their free energy are different. However, in the framework of the bipotential theory, the materials potential can be unified. The energy potential of the ISM cannot be separated. A bipotential was promoted by De Saxcé and Feng in [6], which describes the materials free energy by dual potentials. The bipotential bσ,ε˙ is convex with respect to σ when ε˙ remains constant and convex with respect to ε˙ when σ remains constant [26]. So, there exists(11)GSM:∀σ,ε˙Wσ+Vε˙≥σ:ε˙,ISM:∀σ,ε˙bσ,ε˙≥σ:ε˙.The flow rules are therefore expressed by subdifferential mappings:(12)GSM:σ∈∂Vε˙,ε˙∈∂Wσ,ISM:σ∈∂ε˙bσ,ε˙,ε˙∈∂bσσ,ε˙.As shown in the section above, in the elastic state, the bipotential of the Drucker–Prager model conforms to the characteristics of GSMs. And, in the plastic state, the model satisfies the features of ISMs. So, in order to establish the incremental bipotential function of the nonassociated Drucker–Prager modelΔbepΔσ,Δε, a deriving strategy from the rate form to an incremental one is applied, as shown in Figure 2.Figure 2
Incremental form of the elastoplastic bipotential deriving strategy.


### 3.1. The Rate Form

In the elastic state, the bipotential consists of two convex potentials:Veε˙e which is the strain energy density and Weσ which is the complementary energy density. The superscript “e” stands for the elastic state. Due to (11), the elastic bipotential of the Drucker–Prager model in the rate form can be expressed as(13)beσ,ε˙e=Veε˙e+Weσ≥σ:ε˙.What is more, under the plastic state, the Drucker–Prager model falls into the category of ISMs. According to (11), the rate form of the plastic bipotential is(14)bpσ,ε˙p≥σ:ε˙+ΨKσ∗σ+ΨKε∗ε˙p.Separatingσ⋅ε˙p into a hydrostatic and a deviatoric part, there is(15)σ:ε˙=sm⋅e˙mp+s⋅e˙p+ΨKσ∗σ+ΨKε∗ε˙p,where in the hydrostatic part,(16)sm⋅e˙mp=cϕe˙mp+sm−cϕe˙mp.Considering (1) and (4), (16) becomes(17)sm⋅e˙mp≤cϕe˙mp+kdtanθsm−cϕe˙p.The Cauchy–Schwarz inequality is used to deal with the deviatoric part. By combining it with (1), the following inequality arises:(18)s⋅e˙p≤s⋅e˙p≤−kdtanϕsm−cϕe˙p.Then, from (14), (15), (17), and (18), the rate form of plastic bipotential can be expressed as(19)bpσ,ε˙p=cϕe˙mp+kdtanθ−tanϕsm−cϕe˙p+ΨKσ∗σ+ΨKε∗ε˙p.Takingσ=sm,s and ε˙p=e˙mp,e˙p as partial derivatives to (19), it is proved that the implicit flow laws can be expressed.


### 3.2. The Incremental Form

Assume the plastic strain rate is constant, namely,Δε=Δt⋅ε˙. The incremental elastic bipotential is written as follows.(20)ΔbeΔσ,Δεe=ΔVeΔεe+ΔWeΔσ≥σ:ε˙.Due to the definitions ofΔVeΔεe and ΔWeΔσ,(21)ΔVeΔεe=Δσ:Δεe2=Δsm⋅Δeme+Δs:Δee2=KcΔeme22+μΔee2,ΔWeΔσ=Δσ:Δεe2=Δsm⋅Δeme+Δs:Δee2=Δsm22Kc+4μΔs22Kc.The elastic incremental bipotential of the Drucker–Prager model is(22)ΔbeΔσ,Δεe=KcΔeme22+μΔee2+Δsm22Kc+Δs4μ2.The incremental stress is written asΔσ=σ1−σ0, where the subscript 0 means the initial value of an increment step and the subscript 1 stands for the final value. From the implicit flow laws of (12), there exists(23)Δεp∈∂σbpσ0+Δσ,Δεp,Δσ+σ0∈∂Δεpbpσ0+Δσ,Δεp.The incremental plastic bipotential is written as(24)σ0+Δσ:Δεp≤bpσ0+Δσ,Δεp.Then,(25)Δσ:Δεp≤bpσ0+Δσ,Δεp−σ0:Δεp.Due to (11), it can be also obtained that(26)Δσ:Δεp≤bpΔσ,Δεp.From (24) to (25), the incremental plastic bipotential can be expressed as the following inequality:(27)ΔbpΔσ,Δεp≥bpσ0+Δσ,Δεp−σ0:Δεp.In addition, combining (19) and (27), the bipotential is(28)ΔbpΔσ,Δεp=cϕΔemp+kdtanθ−tanϕsm0+Δsm−cϕΔep−sm0Δemp−s0:Δep+ΨKσ∗σ+ΨKε∗ε˙p.In order to combine the incremental elastic bipotential and the incremental plastic one, a new function is proposed as follows:(29)bx,y=b1x,y−y′©yb2x,y′.The symbol©y in (29) is called the inf-convolution of b1x,y−y′ and b2x,y′ with respect to the space of y. It is defined by the following.(30)b1x,y−y′©yb2x,y′=infy′b1x,y−y′+b2x,y′.Considering the strain decompositionΔεe=Δε−Δεp, the incremental elastoplastic bipotential can be expressed as(31)ΔbepΔσ,Δε=ΔbeΔσ,Δε−Δεp©ΔεpΔbpΔσ,Δεp.We could also state that(32)ΔbepΔσ,Δε=infΔεpΔbeΔσ,Δε−Δεp+ΔbpΔσ,Δεp.And finally, from (22), (28), and (32), the expression of the incremental elastoplastic bipotential is(33)ΔbepΔσ,Δε=infΔemp≥kdtanθΔepKcΔem−Δemp22+μΔe−Δep+cϕΔemp+kdtanθ−tanϕsm0+Δsm−cϕΔep−sm0:Δemp−s0:Δep+Δs4μ2+Δsm22Kc+ΨKσ∗σ+ΨKε∗ε˙p.


## 3.1. The Rate Form

In the elastic state, the bipotential consists of two convex potentials:Veε˙e which is the strain energy density and Weσ which is the complementary energy density. The superscript “e” stands for the elastic state. Due to (11), the elastic bipotential of the Drucker–Prager model in the rate form can be expressed as(13)beσ,ε˙e=Veε˙e+Weσ≥σ:ε˙.What is more, under the plastic state, the Drucker–Prager model falls into the category of ISMs. According to (11), the rate form of the plastic bipotential is(14)bpσ,ε˙p≥σ:ε˙+ΨKσ∗σ+ΨKε∗ε˙p.Separatingσ⋅ε˙p into a hydrostatic and a deviatoric part, there is(15)σ:ε˙=sm⋅e˙mp+s⋅e˙p+ΨKσ∗σ+ΨKε∗ε˙p,where in the hydrostatic part,(16)sm⋅e˙mp=cϕe˙mp+sm−cϕe˙mp.Considering (1) and (4), (16) becomes(17)sm⋅e˙mp≤cϕe˙mp+kdtanθsm−cϕe˙p.The Cauchy–Schwarz inequality is used to deal with the deviatoric part. By combining it with (1), the following inequality arises:(18)s⋅e˙p≤s⋅e˙p≤−kdtanϕsm−cϕe˙p.Then, from (14), (15), (17), and (18), the rate form of plastic bipotential can be expressed as(19)bpσ,ε˙p=cϕe˙mp+kdtanθ−tanϕsm−cϕe˙p+ΨKσ∗σ+ΨKε∗ε˙p.Takingσ=sm,s and ε˙p=e˙mp,e˙p as partial derivatives to (19), it is proved that the implicit flow laws can be expressed.


## 3.2. The Incremental Form

Assume the plastic strain rate is constant, namely,Δε=Δt⋅ε˙. The incremental elastic bipotential is written as follows.(20)ΔbeΔσ,Δεe=ΔVeΔεe+ΔWeΔσ≥σ:ε˙.Due to the definitions ofΔVeΔεe and ΔWeΔσ,(21)ΔVeΔεe=Δσ:Δεe2=Δsm⋅Δeme+Δs:Δee2=KcΔeme22+μΔee2,ΔWeΔσ=Δσ:Δεe2=Δsm⋅Δeme+Δs:Δee2=Δsm22Kc+4μΔs22Kc.The elastic incremental bipotential of the Drucker–Prager model is(22)ΔbeΔσ,Δεe=KcΔeme22+μΔee2+Δsm22Kc+Δs4μ2.The incremental stress is written asΔσ=σ1−σ0, where the subscript 0 means the initial value of an increment step and the subscript 1 stands for the final value. From the implicit flow laws of (12), there exists(23)Δεp∈∂σbpσ0+Δσ,Δεp,Δσ+σ0∈∂Δεpbpσ0+Δσ,Δεp.The incremental plastic bipotential is written as(24)σ0+Δσ:Δεp≤bpσ0+Δσ,Δεp.Then,(25)Δσ:Δεp≤bpσ0+Δσ,Δεp−σ0:Δεp.Due to (11), it can be also obtained that(26)Δσ:Δεp≤bpΔσ,Δεp.From (24) to (25), the incremental plastic bipotential can be expressed as the following inequality:(27)ΔbpΔσ,Δεp≥bpσ0+Δσ,Δεp−σ0:Δεp.In addition, combining (19) and (27), the bipotential is(28)ΔbpΔσ,Δεp=cϕΔemp+kdtanθ−tanϕsm0+Δsm−cϕΔep−sm0Δemp−s0:Δep+ΨKσ∗σ+ΨKε∗ε˙p.In order to combine the incremental elastic bipotential and the incremental plastic one, a new function is proposed as follows:(29)bx,y=b1x,y−y′©yb2x,y′.The symbol©y in (29) is called the inf-convolution of b1x,y−y′ and b2x,y′ with respect to the space of y. It is defined by the following.(30)b1x,y−y′©yb2x,y′=infy′b1x,y−y′+b2x,y′.Considering the strain decompositionΔεe=Δε−Δεp, the incremental elastoplastic bipotential can be expressed as(31)ΔbepΔσ,Δε=ΔbeΔσ,Δε−Δεp©ΔεpΔbpΔσ,Δεp.We could also state that(32)ΔbepΔσ,Δε=infΔεpΔbeΔσ,Δε−Δεp+ΔbpΔσ,Δεp.And finally, from (22), (28), and (32), the expression of the incremental elastoplastic bipotential is(33)ΔbepΔσ,Δε=infΔemp≥kdtanθΔepKcΔem−Δemp22+μΔe−Δep+cϕΔemp+kdtanθ−tanϕsm0+Δsm−cϕΔep−sm0:Δemp−s0:Δep+Δs4μ2+Δsm22Kc+ΨKσ∗σ+ΨKε∗ε˙p.


## 4. Numerical Examples

With the incremental elastoplastic bipotential of the nonassociated Drucker–Prager model, a finite element method is applied to conduct the numerical implementation. We apply the Newton–Raphson method to solve nonlinear equations. The process of the strategy is demonstrated in Figure3. At the Gauss point level, Δσ and Δεp represent the incremental stress and strain tensors. Dep is the local tangent matrix. At the element level, Δinterf is the element internal force and kT is the element tangent stiffness. B stands for the strain matrix, depending on the element shape function. In the global level, ΔR, ΔF, and ΔinterF stand for the incremental residual, external, and internal force vector. U and ΔU represent the total and incremental displacement vector. KT is the total tangent stiffness matrix. The symbol A denotes the assembling operation, and δ1 and δ2 are tolerances of the convergence.Figure 3
Process of the Newton–Raphson strategy.For the numerical implementation, a program was developed in C++. In this section, two numerical examples are simulated to verify the accuracy, stability, and practicability of the proposed approach.


### 4.1. Example 1: Uniaxial Tension and Compression

Simple uniaxial tension and compression examples are simulated to verify the accuracy of the constitutive laws. Because of the perfect elastoplasticity of the nonassociated Drucker–Prager model, we apply a displacementuz on the body to avoid nonconvergence, as shown in Figure 4. A kind of dolomite is chosen in the example. The material parameters are Young’s modulus E=0.5×105MPa, Poisson’s ratio ν=0.33, cohesion c=30MPa, internal friction angle ϕ=40°, and dilatancy angle θ=20°.Figure 4
Uniaxial compression/tension of the Drucker–Prager model example.


#### 4.1.1. Constitutive Curves

The constitutive cone of the nonassociated Drucker–Prager model has been introduced in Section2, and the theoretical cone Kσ is shown in Figure 1(a). In Figure 5, the solid line stands for the simulated results and the dashed line represents the boundary of the stress cone. At the coordinate s/kd−sm, the simulation curves of compression and tension are shown in Figure 5. By comparing the numerical curve with the theoretical one, it appears obvious that the elastic parts are both in the cone Kσ and that the plastic parts are on the boundary of Kσ. This shows that the simulation satisfies the constitutive law of the nonassociated Drucker–Prager model.Figure 5
Comparison of the simulated constitutive law versus the theoretical one: (a) compression; (b) tension.
(a)
(b)Furthermore, dilatancy is a special property in rock-soils. Associated flow rule leads to a nonnegligible volume expansion after the model enters a plastic stage. However, using a nonassociated flow rule is suitable to control the dilated value. Figure6 shows the volume expansion under compression/tension with different dilatancy angles. The constitutive simulation reveals the requirements of soil dilatancy.Figure 6
Volume expansion under different dilatancy angles: (a) compression; (b) tension.
(a)
(b)Dilatancy angle influences the constitutive curves. Figure7 shows the results for different dilatancy angles at the equivalent stress σeq and the equivalent strain εeq. From Figure 7(a), the trends of the curves conform to the characteristic of perfect elastoplasticity. Maximum stress value increases as the dilatancy angle escalates. From Figure 7(b), we can see that the limit stress also increases slightly with a greater dilatancy angle. What is more, a softening phenomenon can be observed when the loading begins. The results satisfy the consequence in [24].Figure 7
Equivalent stress-strain curves for various dilatancy angles: (a) compression; (b) tension.
(a)
(b)


#### 4.1.2. Limit Stress Analysis

In numerical fields, the explicit algorithm introduced in [30] is a typical method for simulating the nonassociated models. The return mapping algorithm is widely used in the simulation of materials constitutive laws [1]. Both methods need the detailed expressions of the yield and the plastic potential functions. In this work, the two potential functions are defined by the yield potential F=sinϕ/33+sin2ϕ1I1+J2−3cosϕ/3+sin2ϕ and the plastic potential Q=sinϕ/33+sin2ϕ2I1+J2−3cosθ/3+sin2ϕ. When ϕ=θ, the plastic potential Q equals F. Therefore, the model is associated. What is more, the literature [31] gives us the analytical equivalent stress under uniaxial compression and tension loadings. Therefore, Figure 8 gives the comparison between the bipotential, explicit, Simo–Taylor, and the theoretical value.Figure 8
Comparison of different algorithms associated with various dilatancy angles based on stress-displacement curves. Compression at (a)θ = 40°, (b) θ = 20°, and (c) θ = 20°. Tension at (d) θ = 40°, (e) θ = 20°, and (f) θ = 0°.
(a)
(b)
(c)
(d)
(e)
(f)Figure8 displays the comparison of 4 algorithms on the compression and tension example for various dilatancy angles for stress-displacement. We limit such curves to 100 steps. As shown in Figure 8(a), if the material is associated and if the transition parts from elasticity to plasticity are neglected, theoretical results are properly satisfied by three methods. As illustrated in Figure 8(b) and Figure 8(c), when the dilatancy angle decreases to 0 degrees, the bipotential results are still satisfying. The error between the return mapping algorithm and the analytical result is smaller than 3.04%. The explicit simulation has a slight oscillation. In Figure 8(d), we note that the peak stress for the three methods is nearly the same as those of the analytical results when the material is associated. From Figures 8(e) and 8(f), we remark that the peak stress calculated by the bipotential method satisfies the theoretical value when the dilatancy angle decreases. The peak value that the return mapping algorithm returns is 5.03% smaller than the analytical one. The result of the explicit algorithm displays oscillations and nonconvergence. From the analysis above, the bipotential algorithm is more accurate and stable in dealing with the nonassociated Drucker–Prager model.


### 4.2. Example 2: The Effect of Deep Footing

Deep footing is a kind of structure foundation that is widely used in civil engineering. A plane strain model is used in this example. The geometry of soils under deep footing is shown in Figure9(a), where l=100m, h=50m, a=12.5m, and b=12.5m. The mesh contains 1289 nodes and 608 6-node triangle elements, as shown in Figure 9(b). A uniform pressure q=16.0kN/m2 is applied on the surface of the foundation trench to equivalently represent the loading of the deep footing. Soil parameters are Young’s modulus E=2.5×106kPa, Poisson’s ratio ν=0.33, soil cohesion c=2kPa, internal friction angle ϕ=30°, and a dilatancy angle of θ=20°.Figure 9
Geometry (a) and mesh (b) of the soil under a deep rigid strip footing.
(a)
(b)Soil mass failure happens when the plastic band generates. Although, in this example, only a small deformation case is considered, the slip bands can still be observed. Here, we study the influence of different dilatancy angles for the slip bands. By locking the maximum plastic strain atεp=0.01%, the slip bands can be expressed by the contour plot of the plastic strain. Under the same value of pressure q=10.24kN, the slip band zones for varying dilatancy angles are illustrated in Figure 10. We note that with the increase of the nonassociated property, the slip band becomes more distinct, which means that the phenomenon of soil mass failure becomes more evident. In other words, when the dilatancy angle decreases, the soil around the deep footing slides easier. This satisfies the results from the typical landslide theory of soils mechanics [5].Figure 10
Comparing the slip bands for various dilatancy angles while locking the maximum plastic strain: (a)θ=ϕ=30°; (b) θ=20°; (c) θ=10°; (d) θ=0°.
(a)
(b)
(c)
(d)


### 4.3. Example 3: Modeling of the Ground Surrounding a Tunnel

This example is concerned by modeling the ground surrounding a tunnel in its cross section. The mesh contains 498 nodes and 225 6-node triangle elements. Both the geometry and the mesh are shown in Figure11. In Figure 11, l=14m, h=20m, and r=6m. The parameters of the ground are E=30,000kN/m2, ν=0.3, c=10kN/m2, and ϕ=20°. A uniform pressure q=10kN/m2 is given on the surface of the ground. The von Mises stress distribution is plotted in Figure 12.Figure 11
Geometry and mesh for the tunnel example.Figure 12
von Mises stress.Comparing for different dilatancy angles, Figure13 gives the vertical displacement Uy curves for the top of the tunnel, and the horizontal displacement Ux curves for the bottom edge of the tunnel.Figure 13
Displacements of the tunnel with different dilatancy angles: (a) vertical displacementUy at the top of the tunnel; (b) horizontal displacement Ux at the bottom edge of the tunnel.
(a)
(b)It can be observed that when the dilatancy angle increases,Ux increases slightly as well, and Uy decreases accordingly. In other words, the larger the θ is, the stiffer the ground is. This is in accordance with what happens in reality for tunnel engineering [32]. Therefore, it is possible to apply the proposed numerical strategy in civil engineering.


## 4.1. Example 1: Uniaxial Tension and Compression

Simple uniaxial tension and compression examples are simulated to verify the accuracy of the constitutive laws. Because of the perfect elastoplasticity of the nonassociated Drucker–Prager model, we apply a displacementuz on the body to avoid nonconvergence, as shown in Figure 4. A kind of dolomite is chosen in the example. The material parameters are Young’s modulus E=0.5×105MPa, Poisson’s ratio ν=0.33, cohesion c=30MPa, internal friction angle ϕ=40°, and dilatancy angle θ=20°.Figure 4
Uniaxial compression/tension of the Drucker–Prager model example.


### 4.1.1. Constitutive Curves

The constitutive cone of the nonassociated Drucker–Prager model has been introduced in Section2, and the theoretical cone Kσ is shown in Figure 1(a). In Figure 5, the solid line stands for the simulated results and the dashed line represents the boundary of the stress cone. At the coordinate s/kd−sm, the simulation curves of compression and tension are shown in Figure 5. By comparing the numerical curve with the theoretical one, it appears obvious that the elastic parts are both in the cone Kσ and that the plastic parts are on the boundary of Kσ. This shows that the simulation satisfies the constitutive law of the nonassociated Drucker–Prager model.Figure 5
Comparison of the simulated constitutive law versus the theoretical one: (a) compression; (b) tension.
(a)
(b)Furthermore, dilatancy is a special property in rock-soils. Associated flow rule leads to a nonnegligible volume expansion after the model enters a plastic stage. However, using a nonassociated flow rule is suitable to control the dilated value. Figure6 shows the volume expansion under compression/tension with different dilatancy angles. The constitutive simulation reveals the requirements of soil dilatancy.Figure 6
Volume expansion under different dilatancy angles: (a) compression; (b) tension.
(a)
(b)Dilatancy angle influences the constitutive curves. Figure7 shows the results for different dilatancy angles at the equivalent stress σeq and the equivalent strain εeq. From Figure 7(a), the trends of the curves conform to the characteristic of perfect elastoplasticity. Maximum stress value increases as the dilatancy angle escalates. From Figure 7(b), we can see that the limit stress also increases slightly with a greater dilatancy angle. What is more, a softening phenomenon can be observed when the loading begins. The results satisfy the consequence in [24].Figure 7
Equivalent stress-strain curves for various dilatancy angles: (a) compression; (b) tension.
(a)
(b)


### 4.1.2. Limit Stress Analysis

In numerical fields, the explicit algorithm introduced in [30] is a typical method for simulating the nonassociated models. The return mapping algorithm is widely used in the simulation of materials constitutive laws [1]. Both methods need the detailed expressions of the yield and the plastic potential functions. In this work, the two potential functions are defined by the yield potential F=sinϕ/33+sin2ϕ1I1+J2−3cosϕ/3+sin2ϕ and the plastic potential Q=sinϕ/33+sin2ϕ2I1+J2−3cosθ/3+sin2ϕ. When ϕ=θ, the plastic potential Q equals F. Therefore, the model is associated. What is more, the literature [31] gives us the analytical equivalent stress under uniaxial compression and tension loadings. Therefore, Figure 8 gives the comparison between the bipotential, explicit, Simo–Taylor, and the theoretical value.Figure 8
Comparison of different algorithms associated with various dilatancy angles based on stress-displacement curves. Compression at (a)θ = 40°, (b) θ = 20°, and (c) θ = 20°. Tension at (d) θ = 40°, (e) θ = 20°, and (f) θ = 0°.
(a)
(b)
(c)
(d)
(e)
(f)Figure8 displays the comparison of 4 algorithms on the compression and tension example for various dilatancy angles for stress-displacement. We limit such curves to 100 steps. As shown in Figure 8(a), if the material is associated and if the transition parts from elasticity to plasticity are neglected, theoretical results are properly satisfied by three methods. As illustrated in Figure 8(b) and Figure 8(c), when the dilatancy angle decreases to 0 degrees, the bipotential results are still satisfying. The error between the return mapping algorithm and the analytical result is smaller than 3.04%. The explicit simulation has a slight oscillation. In Figure 8(d), we note that the peak stress for the three methods is nearly the same as those of the analytical results when the material is associated. From Figures 8(e) and 8(f), we remark that the peak stress calculated by the bipotential method satisfies the theoretical value when the dilatancy angle decreases. The peak value that the return mapping algorithm returns is 5.03% smaller than the analytical one. The result of the explicit algorithm displays oscillations and nonconvergence. From the analysis above, the bipotential algorithm is more accurate and stable in dealing with the nonassociated Drucker–Prager model.


## 4.1.1. Constitutive Curves

The constitutive cone of the nonassociated Drucker–Prager model has been introduced in Section2, and the theoretical cone Kσ is shown in Figure 1(a). In Figure 5, the solid line stands for the simulated results and the dashed line represents the boundary of the stress cone. At the coordinate s/kd−sm, the simulation curves of compression and tension are shown in Figure 5. By comparing the numerical curve with the theoretical one, it appears obvious that the elastic parts are both in the cone Kσ and that the plastic parts are on the boundary of Kσ. This shows that the simulation satisfies the constitutive law of the nonassociated Drucker–Prager model.Figure 5
Comparison of the simulated constitutive law versus the theoretical one: (a) compression; (b) tension.
(a)
(b)Furthermore, dilatancy is a special property in rock-soils. Associated flow rule leads to a nonnegligible volume expansion after the model enters a plastic stage. However, using a nonassociated flow rule is suitable to control the dilated value. Figure6 shows the volume expansion under compression/tension with different dilatancy angles. The constitutive simulation reveals the requirements of soil dilatancy.Figure 6
Volume expansion under different dilatancy angles: (a) compression; (b) tension.
(a)
(b)Dilatancy angle influences the constitutive curves. Figure7 shows the results for different dilatancy angles at the equivalent stress σeq and the equivalent strain εeq. From Figure 7(a), the trends of the curves conform to the characteristic of perfect elastoplasticity. Maximum stress value increases as the dilatancy angle escalates. From Figure 7(b), we can see that the limit stress also increases slightly with a greater dilatancy angle. What is more, a softening phenomenon can be observed when the loading begins. The results satisfy the consequence in [24].Figure 7
Equivalent stress-strain curves for various dilatancy angles: (a) compression; (b) tension.
(a)
(b)


## 4.1.2. Limit Stress Analysis

In numerical fields, the explicit algorithm introduced in [30] is a typical method for simulating the nonassociated models. The return mapping algorithm is widely used in the simulation of materials constitutive laws [1]. Both methods need the detailed expressions of the yield and the plastic potential functions. In this work, the two potential functions are defined by the yield potential F=sinϕ/33+sin2ϕ1I1+J2−3cosϕ/3+sin2ϕ and the plastic potential Q=sinϕ/33+sin2ϕ2I1+J2−3cosθ/3+sin2ϕ. When ϕ=θ, the plastic potential Q equals F. Therefore, the model is associated. What is more, the literature [31] gives us the analytical equivalent stress under uniaxial compression and tension loadings. Therefore, Figure 8 gives the comparison between the bipotential, explicit, Simo–Taylor, and the theoretical value.Figure 8
Comparison of different algorithms associated with various dilatancy angles based on stress-displacement curves. Compression at (a)θ = 40°, (b) θ = 20°, and (c) θ = 20°. Tension at (d) θ = 40°, (e) θ = 20°, and (f) θ = 0°.
(a)
(b)
(c)
(d)
(e)
(f)Figure8 displays the comparison of 4 algorithms on the compression and tension example for various dilatancy angles for stress-displacement. We limit such curves to 100 steps. As shown in Figure 8(a), if the material is associated and if the transition parts from elasticity to plasticity are neglected, theoretical results are properly satisfied by three methods. As illustrated in Figure 8(b) and Figure 8(c), when the dilatancy angle decreases to 0 degrees, the bipotential results are still satisfying. The error between the return mapping algorithm and the analytical result is smaller than 3.04%. The explicit simulation has a slight oscillation. In Figure 8(d), we note that the peak stress for the three methods is nearly the same as those of the analytical results when the material is associated. From Figures 8(e) and 8(f), we remark that the peak stress calculated by the bipotential method satisfies the theoretical value when the dilatancy angle decreases. The peak value that the return mapping algorithm returns is 5.03% smaller than the analytical one. The result of the explicit algorithm displays oscillations and nonconvergence. From the analysis above, the bipotential algorithm is more accurate and stable in dealing with the nonassociated Drucker–Prager model.


## 4.2. Example 2: The Effect of Deep Footing

Deep footing is a kind of structure foundation that is widely used in civil engineering. A plane strain model is used in this example. The geometry of soils under deep footing is shown in Figure9(a), where l=100m, h=50m, a=12.5m, and b=12.5m. The mesh contains 1289 nodes and 608 6-node triangle elements, as shown in Figure 9(b). A uniform pressure q=16.0kN/m2 is applied on the surface of the foundation trench to equivalently represent the loading of the deep footing. Soil parameters are Young’s modulus E=2.5×106kPa, Poisson’s ratio ν=0.33, soil cohesion c=2kPa, internal friction angle ϕ=30°, and a dilatancy angle of θ=20°.Figure 9
Geometry (a) and mesh (b) of the soil under a deep rigid strip footing.
(a)
(b)Soil mass failure happens when the plastic band generates. Although, in this example, only a small deformation case is considered, the slip bands can still be observed. Here, we study the influence of different dilatancy angles for the slip bands. By locking the maximum plastic strain atεp=0.01%, the slip bands can be expressed by the contour plot of the plastic strain. Under the same value of pressure q=10.24kN, the slip band zones for varying dilatancy angles are illustrated in Figure 10. We note that with the increase of the nonassociated property, the slip band becomes more distinct, which means that the phenomenon of soil mass failure becomes more evident. In other words, when the dilatancy angle decreases, the soil around the deep footing slides easier. This satisfies the results from the typical landslide theory of soils mechanics [5].Figure 10
Comparing the slip bands for various dilatancy angles while locking the maximum plastic strain: (a)θ=ϕ=30°; (b) θ=20°; (c) θ=10°; (d) θ=0°.
(a)
(b)
(c)
(d)


## 4.3. Example 3: Modeling of the Ground Surrounding a Tunnel

This example is concerned by modeling the ground surrounding a tunnel in its cross section. The mesh contains 498 nodes and 225 6-node triangle elements. Both the geometry and the mesh are shown in Figure11. In Figure 11, l=14m, h=20m, and r=6m. The parameters of the ground are E=30,000kN/m2, ν=0.3, c=10kN/m2, and ϕ=20°. A uniform pressure q=10kN/m2 is given on the surface of the ground. The von Mises stress distribution is plotted in Figure 12.Figure 11
Geometry and mesh for the tunnel example.Figure 12
von Mises stress.Comparing for different dilatancy angles, Figure13 gives the vertical displacement Uy curves for the top of the tunnel, and the horizontal displacement Ux curves for the bottom edge of the tunnel.Figure 13
Displacements of the tunnel with different dilatancy angles: (a) vertical displacementUy at the top of the tunnel; (b) horizontal displacement Ux at the bottom edge of the tunnel.
(a)
(b)It can be observed that when the dilatancy angle increases,Ux increases slightly as well, and Uy decreases accordingly. In other words, the larger the θ is, the stiffer the ground is. This is in accordance with what happens in reality for tunnel engineering [32]. Therefore, it is possible to apply the proposed numerical strategy in civil engineering.


## 5. Conclusion

In this paper, we have proposed a numerical strategy to simulate the behavior of nonassociated materials. From numerical experiments, we have found the following:(i)
In the simple uniaxial tensile/compressive loading example, the simulated constitutive cones coincide with the theoretical ones as well as the constitutive curves.(ii)
The volume expansion of the material after plastification can be controlled by dilatancy angles.(iii)
From the comparisons of stress-displacement curves calculated with different algorithms, the bipotential method proved to be more stable and more accurate.(iv)
In the deep footing example, slip bands with varying dilatancy angles are compared. When the dilatancy angle decreases, the soil around the deep footing slides easier. The practicability of the bipotential algorithm in civil engineering was therefore verified.(v)
The von Mises stress and displacements of the ground surrounding a tunnel were properly predicted.As a conclusion, the proposed strategy is well suited to study the nonassociated Drucker–Prager model, and the developed program is possible to be applied in civil engineering.



---

*Source: 4534750-2018-05-21.xml*