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

/e 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.


Introduction
e 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.e 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.e 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].
e Uzawa algorithm was used to replace the Newton-Raphson algorithm to solve the nonlinear equations.is 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].e 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].e 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.is 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.
is paper is structured as follows.e cones of the nonassociated Drucker-Prager model are introduced in detail in Section 2. 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.

Nonassociated Drucker-Prager Model
e nonassociated Drucker-Prager model is a classic threeparameter model: the cohesion c, the internal friction angle ϕ, and the dilatancy angle θ [29].e 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.ey are σ � (s m , s) and _ ε p � ( _ e p m , _ e p ), where s m � Tr(σ)/3, s � σ − s m 1, _ e p m � Tr(_ ε p ), and _ e p � _ ε p − _ e p m 1/3, where 1 is the second-order unit tensor.In the elastic state, the stress variable σ � (s m , s) belongs to the constitutive stress cones K σ , as shown in Figure 1(a).e expression of these stress cones is as follows.Advances in Civil Engineering where ‖ • ‖ stands for the Euclidean norm and the k d parameter is a function of the internal friction angle ϕ as 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).e regular and the apex points on the stress cone have different flow expressions.At a regular point, namely, (s m ≠ c/tan ϕ and ‖s‖ ≠ 0), K σ is differentiable, and the flow direction K εr can be expressed as At the apex point, namely, (s m � c/ tan ϕ and ‖s‖ � 0), K σ is nondifferentiable, and the flow direction cannot be determined.A convex set is used here.ere exists a cone K εa such that So, associating ( 2) with ( 3), the plastic rate cone K ε is unified as follows: e implicit flow laws of the nonassociated Drucker-Prager model can be written as follows: where ΨK σ (σ) and ΨK ε (_ ε p ) are indicator functions.Equation ( 5) can also be written with the hydrostatic part and the deviatoric part: From Figure 1(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.ere 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 s m direction by a distance c/tan ϕ, the cone translates from the location (s m , ‖s‖) to (0, 0), as shown in Figure 1(b).A new stress cone K * σ can be defined as with K * ε being the dual cone of K * σ and is defined as According to the geometric features of K ε (_ ε 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 Equation ( 9) is divided into a hydrostatic part and a deviatoric part: 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.

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.e 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.e bipotential b(σ, _ ε) is convex with respect to σ when _ ε remains constant and convex with respect to _ ε when σ remains constant [26].So, there exists e flow rules are therefore expressed by subdifferential mappings: 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 Advances in Civil Engineering incremental bipotential function of the nonassociated Drucker-Prager model Δb ep (Δσ, Δε), a deriving strategy from the rate form to an incremental one is applied, as shown in Figure 2.

3.1.
e Rate Form.In the elastic state, the bipotential consists of two convex potentials: V e (_ ε e ) which is the strain energy density and W e (σ) which is the complementary energy density.
e 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 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 Separating σ • _ ε p into a hydrostatic and a deviatoric part, there is where in the hydrostatic part, Considering ( 1) and ( 4), (16) becomes e Cauchy-Schwarz inequality is used to deal with the deviatoric part.By combining it with (1), the following inequality arises: en, from ( 14), ( 15), (17), and ( 18), the rate form of plastic bipotential can be expressed as Taking σ � (s m , s) and _ ε p � ( _ e p m , _ e p ) as partial derivatives to (19), it is proved that the implicit flow laws can be expressed.

e Incremental
e 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 e incremental plastic bipotential is written as en, Due to (11), it can be also obtained that From ( 24) to (25), the incremental plastic bipotential can be expressed as the following inequality: In addition, combining ( 19) and ( 27), the bipotential is In order to combine the incremental elastic bipotential and the incremental plastic one, a new function is proposed as follows: e symbol © y in ( 29) is called the inf-convolution of b 1 (x, y − y ′ ) and b 2 (x, y ′ ) with respect to the space of y.It is defined by the following.

Numerical Examples
With the incremental elastoplastic bipotential of the nonassociated Drucker-Prager model, a nite element method is applied to conduct the numerical implementation.We apply the Newton-Raphson method to solve nonlinear equations.e process of the strategy is demonstrated in Figure 3.At the Gauss point level, Δσ and Δε p represent the incremental stress and strain tensors.D ep is the local tangent matrix.At the element level, Δinter f is the element internal force and k T is the element tangent sti ness.B stands for the strain matrix, depending on the element shape function.In the global level, ΔR, ΔF, and Δinter F stand for the incremental residual, external, and internal force vector.U and ΔU represent the total and incremental displacement vector.K T is the total tangent sti ness matrix.e symbol A()  Advances in Civil Engineering denotes the assembling operation, and δ 1 and δ 2 are tolerances of the convergence.
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.

Constitutive Curves.
e constitutive cone of the nonassociated Drucker-Prager model has been introduced in Section 2, 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 /k d − s m , 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 σ . is shows that the simulation satis es the constitutive law of the nonassociated Drucker-Prager model.
Furthermore, dilatancy is a special property in rock-soils.Associated ow rule leads to a nonnegligible volume expansion after the model enters a plastic stage.However, using a nonassociated ow rule is suitable to control the dilated value.Figure 6 shows the volume expansion under compression/tension with di erent dilatancy angles.e  Dilatancy angle in uences the constitutive curves.Figure 7 shows the results for di erent 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.e results satisfy the consequence in [24].

Limit Stress Analysis.
In numerical elds, the explicit algorithm introduced in [30] is a typical method for simulating the nonassociated models.e 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 de ned by the yield potential and the plastic potential When ϕ θ, the plastic potential Q equals F. erefore, the model is associated.What is more, the literature [31] gives us the analytical equivalent stress under uniaxial compression and tension loadings.erefore, Figure 8 gives the comparison between the bipotential, explicit, Simo-Taylor, and the theoretical value.
Figure 8 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 satis ed 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.e error between the return mapping algorithm and the analytical result is smaller than 3.04%.e 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 satis es the theoretical value when the dilatancy angle decreases.e peak value that the return mapping algorithm returns is 5.03% smaller than the analytical one.e 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.

Example 2: e E ect 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.e geometry of soils under deep footing is shown in Figure 9(a), where l 100 m, h 50 m, a 12.5 m, and b 12.5 m. e mesh contains 1289 nodes and 608 6-node triangle elements, as shown in Figure 9(b).A uniform pressure q 16.0 kN/m 2 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 × 10 6 kPa, Poisson's ratio ] 0.33, soil cohesion c 2 kPa, internal friction angle ϕ 30 °, and a dilatancy angle of θ 20 °.
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 in uence of di erent 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.24 kN, 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.is satis es the results from the typical landslide theory of soils mechanics [5].

Example 3: Modeling of the Ground Surrounding a Tunnel.
is example is concerned by modeling the ground surrounding a tunnel in its cross section.e mesh contains 498 nodes and 225 6-node triangle elements.Both the geometry and the mesh are shown in Figure 11.In Figure 11, l 14 m, h 20 m, and r 6 m. e parameters of the ground are E 30, 000 kN/m 2 , ] 0.3, c 10 kN/m 2 , and ϕ 20 °.A uniform pressure q 10kN/m 2 is given on the surface of the ground.
e von Mises stress distribution is plotted in Figure 12.
Comparing for di erent dilatancy angles, Figure 13 gives the vertical displacement U y curves for the top of the tunnel, and the horizontal displacement U x curves for the bottom edge of the tunnel.
It can be observed that when the dilatancy angle increases, U x increases slightly as well, and U y decreases accordingly.In other words, the larger the θ is, the sti er the ground is. is is in accordance with what happens in reality for tunnel engineering [32].erefore, it is possible to apply the proposed numerical strategy in civil engineering.Advances in Civil Engineering

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) e volume expansion of the material after plastication can be controlled by dilatancy angles.(iii) From the comparisons of stress-displacement curves calculated with di erent 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.e practicability of the bipotential algorithm in civil engineering was therefore veri ed.(v) e 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.

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 displacement u z on the body to avoid nonconvergence, as shown in Figure 4.A kind of dolomite is chosen in the example.e material parameters are Young's modulus E 0.5 × 10 5 MPa, Poisson's ratio ] 0.33, cohesion c 30 MPa, internal friction angle ϕ 40 °, and dilatancy angle θ 20 °.

Figure 9 :Figure 10 :
Figure 9: Geometry (a) and mesh (b) of the soil under a deep rigid strip footing.

Figure 11 :
Figure 11: Geometry and mesh for the tunnel example.

Figure 13 :Figure 12 :
Figure 13: Displacements of the tunnel with di erent dilatancy angles: (a) vertical displacement U y at the top of the tunnel; (b) horizontal displacement U x at the bottom edge of the tunnel.