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.

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

2. Nonassociated Drucker–Prager Model

The nonassociated Drucker–Prager model is a classic three-parameter model: the cohesion , 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 and , where , , , and , where is the second-order unit tensor. In the elastic state, the stress variable belongs to the constitutive stress cones , as shown in Figure 1(a). The expression of these stress cones is as follows.where stands for the Euclidean norm and the 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 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, , is differentiable, and the flow direction can be expressed as

At the apex point, namely, , is nondifferentiable, and the flow direction cannot be determined. A convex set is used here. There exists a cone such that

So, associating (2) with (3), the plastic rate cone is unified as follows:

The implicit flow laws of the nonassociated Drucker–Prager model can be written as follows:where and 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 , and 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 direction by a distance , the cone translates from the location to , as shown in Figure 1(b). A new stress cone can be defined aswith being the dual cone of and is defined as

According to the geometric features of , , , and and combining them with (6), while setting , 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 , 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 is convex with respect to when remains constant and convex with respect to when remains constant [26]. So, there exists

The 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 incremental bipotential function of the nonassociated Drucker–Prager model , a deriving strategy from the rate form to an incremental one is applied, as shown in Figure 2.

3.1. The Rate Form

In the elastic state, the bipotential consists of two convex potentials: which is the strain energy density and 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

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 into a hydrostatic and a deviatoric part, there iswhere in the hydrostatic part,

Considering (1) and (4), (16) becomes

The Cauchy–Schwarz inequality is used to deal with the deviatoric part. By combining it with (1), the following inequality arises:

Then, from (14), (15), (17), and (18), the rate form of plastic bipotential can be expressed as

Taking and 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, . The incremental elastic bipotential is written as follows.

Due to the definitions of and ,

The elastic incremental bipotential of the Drucker–Prager model is

The incremental stress is written as , 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

The incremental plastic bipotential is written as

Then,

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:

The symbol in (29) is called the inf-convolution of and with respect to the space of . It is defined by the following.

Considering the strain decomposition , the incremental elastoplastic bipotential can be expressed as

We could also state that

And finally, from (22), (28), and (32), the expression of the incremental elastoplastic bipotential is

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 Figure 3. At the Gauss point level, and represent the incremental stress and strain tensors. is the local tangent matrix. At the element level, is the element internal force and is the element tangent stiffness. stands for the strain matrix, depending on the element shape function. In the global level, , , and stand for the incremental residual, external, and internal force vector. and represent the total and incremental displacement vector. is the total tangent stiffness matrix. The symbol denotes the assembling operation, and and 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.

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 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 , Poisson’s ratio , cohesion , internal friction angle , and dilatancy angle .

4.1.1. Constitutive Curves

The constitutive cone of the nonassociated Drucker–Prager model has been introduced in Section 2, and the theoretical cone 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 , 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 and that the plastic parts are on the boundary of . This shows that the simulation satisfies the constitutive law of the nonassociated Drucker–Prager model.

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. Figure 6 shows the volume expansion under compression/tension with different dilatancy angles. The constitutive simulation reveals the requirements of soil dilatancy.

Dilatancy angle influences the constitutive curves. Figure 7 shows the results for different dilatancy angles at the equivalent stress and the equivalent strain . 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].

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 and the plastic potential . When , the plastic potential equals . 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 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 Figure 9(a), where , , , and . The mesh contains 1289 nodes and 608 6-node triangle elements, as shown in Figure 9(b). A uniform pressure is applied on the surface of the foundation trench to equivalently represent the loading of the deep footing. Soil parameters are Young’s modulus , Poisson’s ratio , soil cohesion , internal friction angle , and a dilatancy angle of .

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 , the slip bands can be expressed by the contour plot of the plastic strain. Under the same value of pressure , 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].

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 Figure 11. In Figure 11, , , and . The parameters of the ground are , , , and . A uniform pressure is given on the surface of the ground. The von Mises stress distribution is plotted in Figure 12.

Comparing for different dilatancy angles, Figure 13 gives the vertical displacement curves for the top of the tunnel, and the horizontal displacement curves for the bottom edge of the tunnel.

It can be observed that when the dilatancy angle increases, increases slightly as well, and 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.

Conflicts of Interest

The authors declare that all funding programmes do not lead to any conflict of interests regarding the publication of the paper.

Acknowledgments

The authors gratefully acknowledge the financial support of the National Natural Science Foundation of China (Grant no. 11772274) and the National Key R&D Program of China (Grant no. 2017YFB0703200).