# Non-associated constitutive law of soils and its simulation based on the bi-potential theory

Y.J. Zhou, Z.Q. Feng, W.Y. Xu

► To cite this version:

Y.J. Zhou, Z.Q. Feng, W.Y. Xu. Non-associated constitutive law of soils and its simulation based on the bi-potential theory. International Journal of Structural Analysis and Design, 2014, 1 (4), pp.1-6. hal-02398148

HAL Id: hal-02398148

https://univ-evry.hal.science/hal-02398148

Submitted on 4 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Non-associated Constitutive Law of Soils and its Simulation Based on the Bi-potential Theory*

Y-J. Zhou<sup>a</sup>, Z-Q. Feng<sup>a,b</sup>, W-Y. Xu<sup>a</sup>

**Abstract**—This paper is devoted to study the non-associated constitutive law of soil materials and to construct a bi-potential for the classical Drucker-Prager (D-P) plastic model. It can be proved that the bi-potential not only satisfies the constitutive relationship of the D-P model, but also conforms to the non-associated flow rule. We present the relevant algorithms for the integration of soil constitutive equations. Two simple numerical examples illustrate the rationality of the bi-potential theory and the feasibility of its algorithm.

**Keywords**—Bi-potential, D-P model, Algorithm, Simulation

## I. Introduction

Material nonlinearity plays an important role in solid mechanics. The rationality and accuracy of material nonlinearity theory are the basis of theoretical analysis and simulation of materials and structures. For classical materials such as metal, the yield potential is associated with the flow potential. Its theory is relatively well established. But for the non-associated materials such as soils, the yield potential is different from the flow one. Conventional model uses two independent functions to express the potentials [1]. Nowadays, the theory on associated materials is relatively well established and its simulation becomes mature. The return-mapping algorithm is a widely accepted and applied method to integrate constitutive laws of associated materials [2]. For non-associated materials, the yield potential is different from the flow one. De Saxcé uses Legendre transformation to deal with this kind of laws [3]. The bi-potential theory is put forward in [4] and the theory is applied to the non-associated material successfully. In [5] and [6], the bi-potential function has been proved to be able to express the flow rule of soils more clearly. The bi-potential theory is suitable to solve the set problems and can be used to build the constitutive laws of non-associate materials. Material like soils is a kind of non-associated materials, whose mechanical behavior is very complex. This paper summarizes previous relevant works, choosing the classical three-parameter soil D-P model to build its constitutive relationship under within the bi-potential framework. This theory leads to a simple and efficient algorithm for the numerical integration of non-associated constitutive laws.

## II. D-P model

The D-P model is given in [7] and illustrated in Figure 1. In order to simplify the problem, the ideal plastic D-P model is used in this paper. This model is represented by a pair of conjugate tensor $\sigma = (s, s_m)$ and $\dot{\epsilon}^p = (\dot{e}^p, \dot{e}_m^p)$, in which $s_m = 1/3 \cdot Tr(\sigma)$ is the mean stress, $s = \sigma - s_m \mathbf{1}$ the stress deviator tensor, $\dot{e}_m^p = Tr(\dot{\epsilon}_m^p)$ the trace of the plastic strain rate tensor, $\dot{e}^p = \dot{e}^p - 1/3 \cdot \dot{e}_m^p \mathbf{1}$ the deviator tensor and $\mathbf{1}$ the second order unity tensor.

During all the process, the stress must be in the stress cone $K_a$ defined by

$$ K_a = \{ (s_m, s) / \|s\| / k_d + s_m \tan \phi \le c \} \quad (1) $$

where $c$ is the cohesion, $\phi$ the friction angle, $k_d$ a constant, $\|\|$ the Euclidean norm.

The flow rule of D-P model is limited in the plastic strain rate cone $K_c$, which can be represented by two parts: the regular set $K_r$ and the apex set $K_a$.

$$ K_r = \{ (\dot{e}_m^p, \dot{e}_m^r) / \dot{e}_m^r \ge k_d \tan \theta \|\dot{e}^p\| \} \quad (2) $$

$$ K_a = \{ (\dot{e}_m^p, \dot{e}_m^r) / \dot{e}_m^r \ge k_d \tan \phi \|\dot{e}^p\| \} \quad (3) $$

<figure><img src="image_2.png" /><figcaption>Figure 1 D-P model cones</figcaption></figure>

As is known in [7], the D-P model is the implicit standard material. It is clear that the plastic strain rate is not orthogonal to the yield surface, so the conventional orthogonal law is not applicable. Form above, the theory of bi-potential gives the strain rate component of projection,

$$ \dot{\epsilon}^p = (\dot{e}^p, \dot{e}_m^p + k_d(\tan \phi - \tan \theta) \|\dot{e}^p\|) \quad (4) $$

<footnote>* This work is supported by and performed within the National Natural Science Foundation of China (Grant N° 11372260).</footnote>

<footnote><sup>a</sup> Southwest Jiaotong University, China</footnote>

<footnote><sup>b</sup> Université d'Evry-Val d'Essonne, France</footnote>
Equation (4) can be written as a differential collection set form to express the flow rules of D-P model. The detailed proofs are given in [9].

$$ (\dot{\mathbf{e}}^p, \dot{\mathbf{e}}_m^p + kd(\tan \phi - \tan \theta)) \|\dot{\mathbf{e}}^p\| \in \partial \Psi_{\sigma, \mathbf{c}}(\sigma) \quad (5) $$

### III. Bi-potential for soils

In [3], the bi-potential for both explicit and implicit standard material is introduced in detail. So, the rate form and increment form bi-potential for D-P model can be obtained.

#### A. Rate form

##### 1) Elasticity

When D-P model is in elastic state, it can be considered as an explicit standard material. At this time

$$ b_c(\dot{\mathbf{\epsilon}}^c, \mathbf{\sigma}) = \varphi_c(\dot{\mathbf{\epsilon}}^c) + \varphi_c^*(\mathbf{\sigma}) \geq \dot{\mathbf{\epsilon}}^c \cdot \mathbf{\sigma} \quad (6) $$

Moreover, the model satisfies the elastic constitutive laws

$$ \varphi_c(\dot{\mathbf{\epsilon}}^c) = \dot{\mathbf{\epsilon}}^c [\mathbf{D}^c] \dot{\mathbf{\epsilon}}^c / 2 = K_c (\dot{e}_m^c)^2 / 2 + \mu \|\dot{\mathbf{e}}^c\|^2 \quad (7) $$

$$ \varphi_c^*(\sigma) = \sigma [\mathbf{D}^c]^{-1} \sigma / 2 = (s_m)^2 / 2 K_c + \|\mathbf{s}\|^2 / 4 \mu \quad (8) $$

The bi-potential function can be obtained under the elastic state,

$$ b_c(\dot{\mathbf{\epsilon}}^c, \mathbf{\sigma}) = K_c (\dot{e}_m^c)^2 / 2 + \mu \|\dot{\mathbf{e}}^c\|^2 + (s_m)^2 / 2 K_c + \|\mathbf{s}\|^2 / 4 \mu \quad (9) $$

Where $K_c$ and $\mu$ are respectively the bulk and shear modulus. Equation (9) can be easily verified to meet the orthogonal law

$$ \dot{\mathbf{\epsilon}}^c = \partial_c \varphi_c^*(\sigma) = \partial_c b_c(\dot{\mathbf{\epsilon}}^c, \sigma) \quad \text{and} \quad \sigma = \partial_c \varphi_c(\dot{\mathbf{\epsilon}}^c) = \partial_c b_c(\dot{\mathbf{\epsilon}}^c, \sigma) \quad (10) $$

##### 2) Plasticity

When D-P model is in plastic state, it can be considered as an implicit standard material. At this time,

$$ b_p(\dot{\mathbf{\epsilon}}^p, \mathbf{\sigma}) \geq \dot{\mathbf{\epsilon}}^p \cdot \mathbf{\sigma} \quad (11) $$

$\dot{\mathbf{\epsilon}}^p \cdot \mathbf{\sigma}$ can be decomposed as,

$$ \mathbf{\sigma} \cdot \dot{\mathbf{\epsilon}}^p = s_m \cdot \dot{e}_m^p + \mathbf{s} \cdot \dot{\mathbf{e}}^p \quad (12) $$

The first item of (12) can be rewritten as follow,

$$ s_m \cdot \dot{e}_m^p = c \dot{e}_m^p / \tan \phi + (s_m - c / \tan \phi) \dot{e}_m^p \quad (13) $$

Equation (1) and (2) give,

$$ s_m \cdot \dot{e}_m^p \leq c \dot{e}_m^p / \tan \phi + k_d \tan \theta (s_m - c / \tan \phi) \|\dot{e}_m^p\| \quad (14) $$

For the second item of (12), considering the Cauchy-Schwarz inequality and the relationship (1), another inequality is given by

$$ s \cdot \dot{e}^p \leq \|\mathbf{s}\| \cdot \|\dot{\mathbf{e}}^p\| \leq -k_d \tan \phi (s_m - c / \tan \phi) \|\dot{\mathbf{e}}^p\| \quad (15) $$

So, in the case of plastic flow, the plastic bi-potential can be obtained by combining (11), (12), (13) and (14).

$$ b_p(\dot{\mathbf{\epsilon}}^p, \mathbf{\sigma}) = c \dot{e}_m^p / \tan \phi + k_d (\tan \theta - \tan \phi) (s_m - c / \tan \phi) \|\dot{\mathbf{e}}^p\| \quad (16) $$

The flow rule of D-P model under plastic stage could also be verified.

$$ \dot{\mathbf{\epsilon}}^p \in \partial_c b_p(\dot{\mathbf{\epsilon}}^p, \mathbf{\sigma}) \quad (17) $$

#### B. Incremental form

According to the incremental relationships

$$ \Delta \mathbf{\epsilon} = \Delta t \cdot \dot{\mathbf{\epsilon}} \quad , \quad \Delta \mathbf{\sigma} = \mathbf{\sigma}_1 - \mathbf{\sigma}_0 \quad (18) $$

The subscript 0 indicates the initial iteration and the subscript 1 stands for the final iteration. According to relationships above, the incremental bi-potential function of the D-P model can be built.

##### 1) Elasticity

According to (6) and (18), $\Delta b_c(\Delta \mathbf{\epsilon}^c, \Delta \mathbf{\sigma})$ is defined as

$$ \Delta b_c(\Delta \mathbf{\epsilon}^c, \Delta \mathbf{\sigma}) = \Delta \varphi_c(\Delta \mathbf{\epsilon}^c) + \Delta \varphi_c^*(\Delta \mathbf{\sigma}) \geq \Delta \mathbf{\sigma} \cdot \Delta \mathbf{\epsilon}^c \quad (19) $$

So, under the condition of elastic constitutive, the incremental elastic bi-potential of D-P model is

$$ \Delta b_c(\Delta \mathbf{\epsilon}^c, \Delta \mathbf{\sigma}) = \frac{K_c}{2} (\Delta \mathbf{e}_m^c)^2 + \mu \|\Delta \mathbf{e}^c\|^2 + \frac{1}{2 K_c} (\Delta s_m)^2 + \frac{1}{4 \mu} \|\Delta \mathbf{s}\|^2 \quad (20) $$

##### 2) Plasticity

In the period of time $\Delta t$, the plastic bi-potential must satisfy

$$ \Delta b_p(\Delta \mathbf{\epsilon}^p, \Delta \mathbf{\sigma}) \geq \Delta \mathbf{\epsilon}^p \cdot \Delta \mathbf{\sigma} \quad (21) $$

Furthermore,

$$ \Delta \mathbf{\epsilon}^p \in \partial_c b_p(\Delta \mathbf{\epsilon}^p, \mathbf{\sigma}_0 + \Delta \mathbf{\sigma}) \quad (22) $$

$$ \Delta \mathbf{\sigma} \in \partial_c b_p(\Delta \mathbf{\epsilon}^p, \mathbf{\sigma}_0 + \Delta \mathbf{\sigma}) - \mathbf{\sigma}_0 \quad (23) $$

Thus the $b_p(\Delta \mathbf{\epsilon}^p, \Delta \mathbf{\sigma})$ can be expressed as

$$ \Delta b_p(\Delta \mathbf{\epsilon}^p, \Delta \mathbf{\sigma}) \geq b_p(\Delta \mathbf{\epsilon}^p, \mathbf{\sigma}_0 + \Delta \mathbf{\sigma}) - \mathbf{\sigma}_0 \cdot \Delta \mathbf{\epsilon}^p \quad (24) $$

So, under the condition of plastic flow, the incremental plastic bi-potential of the D-P model is given by

$$ b_p(\Delta \mathbf{\epsilon}^p, \Delta \mathbf{\sigma}) = c \Delta \mathbf{e}_m^p / \tan \phi + k_d (\tan \theta - \tan \phi) (s_m + \Delta s_m - c / \tan \phi) \|\Delta \mathbf{e}^p\| - s_m \Delta \mathbf{e}_m^p - s_0 \Delta \mathbf{e}^p \quad (25) $$

##### 3) Elastoplasticity
There exists an incremental elastoplastic relationship

$$ \Delta \varepsilon^r = \Delta \varepsilon - \Delta \varepsilon^p \quad (26) $$

It can be implied from (20) and (21) the incremental elastoplastic bi-potential,

$$ \Delta b(\Delta \varepsilon, \Delta \sigma) = \begin{aligned} & \inf_{s_{n+1}^r; \varepsilon_n^r \in \Gamma_n^r} \left\{ \begin{aligned} & K_c(\Delta e_m - \Delta e_m^r)^2 / 2 + \mu \left\| \Delta \mathbf{e} - \Delta \mathbf{e}^p \right\| - s_{n+1} \Delta e_m^p \\ & - s_n \Delta \mathbf{e}^p + c \Delta e_m^p / \tan \varphi + k_d(\tan \theta - \tan \varphi)(s_{n+1} + \Delta s_m - c / \tan \varphi) \end{aligned} \right\} \end{aligned} \quad (27) $$

# IV. Algorithm of the D-P model

## A. Finite element method

According to the bi-potential theory and the equilibrium equation we have

$$ \int_V \mathbf{B}^T \frac{\partial \Delta b[\Delta \varepsilon(\Delta \mathbf{u}), \Delta \sigma]}{\partial \Delta \varepsilon} dV - \int_V \mathbf{N}^T \Delta \mathbf{f} dV - \int_S \mathbf{N}^T \Delta \mathbf{T} dS = \mathbf{0} \quad (28) $$

In this equation, the variable is the node incremental displacement, which can be solved by the finite element method, and then the incremental strain field can be obtained. Furthermore, the incremental stress field can be obtained by integrating the constitutive laws.

For the structure in each loading step, the incremental residual force vector should converge to zero, namely

$$ \Delta \mathbf{F}_r = \int_{\Omega} \mathbf{B}^T \Delta \sigma - \Delta \mathbf{F} \to \mathbf{0} \quad (29) $$

For a single Gauss point, the incremental residual stress tensor should converges to zero too

$$ \Delta \sigma_r = \partial \Delta b(\Delta \varepsilon, \Delta \sigma) / \partial \Delta \varepsilon - \Delta \sigma \to 0 \quad (30) $$

<figure><img src="image_3.png" /><figcaption>Figure 2 (ΔU, Δσ) updates</figcaption></figure>

Figure 2 shows the solution procedure. $A_d$ is the $(\Delta U, \Delta \sigma)$ in the trial conditions, $\Gamma$ represents $(\Delta U, \Delta \sigma)$ of the real value. $\mathbf{D}^+$ realizes the update of incremental stress on a single Gauss point, explaining (30). $\mathbf{D}^-$ realizes the update of incremental displacement of the global structure, explaining (29). Finally, after numbers of iterations, the incremental stress on every Gauss point should satisfy the D-P model, and the incremental displacement should make the equilibrium

equation established. This way of updating $(\Delta U, \Delta \sigma)$ is called LATIN method [8].

## B. Extremism conditions

Dealing (27) with the Lagrange multiplier method and considering the condition (2), a Lagrange expression is given as

$$ \begin{aligned} L(\Delta \varepsilon, \Delta \sigma, \lambda) = & \lambda (k_d \tan \theta \| \Delta \mathbf{e}^p \| - \Delta e_m^p) + K_c (\Delta e_m - \Delta e_m^p)^2 / 2 \\ & + \mu \| \Delta \mathbf{e} - \Delta \mathbf{e}^p \|^2 + c \Delta e_m^p / \tan \varphi + k_d (\tan \theta - \tan \varphi) \\ & \cdot (s_{n+1} + \Delta s_m - c / \tan \varphi) \| \Delta \mathbf{e}^p \| - s_{n+1} \Delta e_m^p - s_n \Delta e_m^p \end{aligned} \quad (31) $$

So, there are three extremities, which can be shown as follows

$$ \partial L / \partial \Delta \mathbf{e}^p = 0, \quad \partial L / \partial \Delta e_m^p = 0, \quad \lambda \cdot \partial L / \partial \lambda = 0 \quad (32) $$

According to the extremisms above and the consistency condition, the third item of (32) can be expressed as following relationships [9]

$$ k_d \tan \theta \| \Delta \mathbf{e}^p \| - \Delta e_m^p = 0 \quad \text{and} \quad \lambda > 0 \quad (33) $$

$$ \text{or} \quad k_d \tan \theta \| \Delta \mathbf{e}^p \| - \Delta e_m^p < 0 \quad \text{and} \quad \lambda = 0 \quad (34) $$

## C. Update $(\Delta \varepsilon^p, \Delta \sigma)$

According to (33) and (34), the $(\Delta \varepsilon^p, \Delta \sigma)$ updating method is obtained. The update method can be discussed in two situations: regular point and apex point.

### 1) Regular point

At elastic stage, plastic strain disappears, $\Delta s$ and $\Delta s_m$ can be explicitly calculated by

$$ \Delta s = 2\mu\Delta\mathbf{e} \quad \text{and} \quad \Delta s_m = K_c \Delta e_m \quad (35) $$

Under the plastic phase, (33) is satisfied, then

$$ \| \Delta \mathbf{e}^p \| = \frac{\| \mathbf{n}_o + \Delta \mathbf{e} \| - \varepsilon_d \{ c - \tan \varphi (s_{n+1} + \Delta s_m) + \tan \theta (\Delta s_m - K_c \Delta e_m) \}}{1 + \varepsilon_d k_d^2 \tan^2 \theta} \quad (36) $$

$$ \Delta \mathbf{e}^p = \| \Delta \mathbf{e}^p \| \mathbf{n} \quad \text{and} \quad \Delta e_m^p = k_d \tan \theta \| \Delta \mathbf{e}^p \| \quad (37) $$

$$ \text{where} \quad \varepsilon_d = \frac{k_d}{2\mu}, \quad \varepsilon_c = \frac{K_c}{2\mu}, \quad \mathbf{n}_o = \frac{\mathbf{s}_o}{2\mu}, \quad \mathbf{n} = \frac{\mathbf{n}_o + \Delta \mathbf{e}}{\| \mathbf{n}_o + \Delta \mathbf{e} \|} $$

As we known, $\Delta s$ and $\Delta s_m$ are expressed as

$$ \Delta s = 2\mu(\Delta \mathbf{e} - \Delta \mathbf{e}^p) \quad \text{and} \quad \Delta s_m = K_c(\Delta e_m - \Delta e_m^p) \quad (38) $$

According to (36), (38) and the implicit expressions of $\Delta s$ and $\Delta s_m$, the residual stress $\Delta s$, and $\Delta s_m$ should finally tend to 0 and 0 by implicit iterations

$$ \Delta s_r = 2\mu (\| \Delta \mathbf{e} + \mathbf{n}_o \| - \| \Delta \mathbf{e}^p \|) \mathbf{n} - (\mathbf{s}_o + \Delta s) \to \mathbf{0} \quad (39) $$
$$ \Delta s_m = K_c (\Delta e_m - \Delta e_m^p) - \Delta s_m \to 0 \quad (40) $$

Because $||\Delta e^p||$ must be nonnegative, so we can get the plastic yield criterion at regular point of the D-P model

$$ ||\mathbf{n}_0 + \Delta \mathbf{e}|| - \varepsilon_d \{c - \tan \varphi (s_{m0} + \Delta s_m) + \tan \theta (\Delta s_m - K_c \Delta e_m)\} \ge 0 \quad (41) $$

## 2) Apex point

An apex point has the following characteristics

$$ s_m = c / \tan \varphi \quad \text{and} \quad ||s|| = 0 \quad (42) $$

Updating $(\Delta\varepsilon^p, \Delta\sigma)$ at an apex point in elastic stage is the same as that at a regular point.

When the apex point under the plastic phase, (42) is satisfied, then

$$ ||\Delta \mathbf{e}^p|| = ||\mathbf{n}_0 + \Delta \mathbf{e}^p|| = ||\mathbf{s}_0 + 2 \mu \Delta \mathbf{e}^p|| / 2 \mu \quad (43) $$

So, $(\Delta\varepsilon^p, \Delta\sigma)$ can be explicitly expressed,

$$ \Delta \mathbf{e}^p = ||\Delta \mathbf{e}^p|| \mathbf{n} \quad \text{and} \quad \Delta e_m^p = (s_{m0} + K_c \Delta e_m - c / \tan \varphi) / K_c \quad (44) $$

$$ \Delta \mathbf{s} = 2\mu(\Delta \mathbf{e} - \Delta \mathbf{e}^p) \quad \text{and} \quad \Delta s_m = K_c (\Delta e_m - \Delta e_m^p) \quad (45) $$

Because $\partial L / \partial \Delta e_m^p = 0$, the plastic yield criterion at apex point of the D-P model is given by

$$ k_d \tan \theta ||\mathbf{n}_0 + \Delta \mathbf{e}^p|| - (s_{m0} + K_c \Delta e_m - c / \tan \varphi) / K_c \le 0 \quad (46) $$

# D. Elastoplastic matrix

According to the definition of elastoplastic matrix

$$ \mathbf{D}_{ep} = \partial \Delta \boldsymbol{\sigma} / \partial \Delta \varepsilon \quad (47) $$

Because of the implicit expression of the incremental stress and incremental strain, the elastoplastic matrix at regular point and apex point can be represented as a unity

$$ \mathbf{D}_{ep} = \frac{\partial \Delta \boldsymbol{\sigma}}{\partial \Delta \boldsymbol{\sigma}_r} \frac{\partial \Delta \boldsymbol{\sigma}_r}{\partial \Delta \varepsilon} = - \left( \frac{\partial \Delta \boldsymbol{\sigma}_r}{\partial \Delta \boldsymbol{\sigma}} \right)^{-1} \frac{\partial \Delta \boldsymbol{\sigma}_r}{\partial \Delta \varepsilon} = (\mathbf{I} - \mathbf{D}_{ep})^{-1} \mathbf{D}_{ep} \quad (48) $$

## 1) Regular point

Under elastic stage, the elastoplastic matrix is equal to the elastic matrix.

$$ \mathbf{D}_{ep} = \mathbf{D}_{ep} \quad (49) $$

Under plastic stage, the elastoplastic matrix can be calculated by (48) and

$$ \mathbf{D}_i = \beta \mathbf{1} \otimes \mathbf{1} / k_d + \alpha (\mathbf{I} - \mathbf{1} \otimes \mathbf{1} / 3) + \delta \mathbf{n} \otimes \mathbf{n} \\ - \beta \tan \theta (\mathbf{1} \otimes \mathbf{n} + \mathbf{n} \otimes \mathbf{1}) \quad (50) $$

$$ \mathbf{D}_c = (\beta(\tan \theta - \tan \varphi)\mathbf{n} \otimes \mathbf{1}) / 3K_c \quad (51) $$

$$ \text{with } \alpha = 2\mu \left( 1 - \frac{||\Delta \mathbf{e}^p||}{||\Delta \mathbf{e} + \mathbf{n}_0||} \right), \quad \beta = \frac{1}{1 + \varepsilon_d k_d^2 \tan^2 \theta}, \quad \delta = k_d \beta \tan^2 \theta. $$

## 2) Apex point

Due to the constraints (42), the elastoplastic matrix at the apex point can be inferred as

$$ \mathbf{D}_{ep} = \mathbf{0} \quad (52) $$

# E. Algorithms

## 1) Constitutive integration

The constitutive integration located at one Gauss point is illustrated in Figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3 Constitutive integration algorithm</figcaption></figure>

## 2) Structure

The algorithm for a structure in one iteration is shown in Figure 4.

<figure><img src="image_5.png" /><figcaption>Figure 4 Algorithm for structure in one step</figcaption></figure>
# V. Simulation examples

## A. Compression/traction test

### 1) Numerical solution

Let us take a 4-node unit as an example under a uniform compression and traction condition. The aim of this example is to check the correctness of integrating constitutive equations on one stress point, i.e., Gauss point. The characteristics are:
$E = 0.5 \times 10^5$ MPa, $\nu = 0.33$, $c = 30$ MPa, $k_d = 1.01566$, $\phi = 40^\circ$, $\theta = 40^\circ$. Because the D-P model is an ideal plastic model, displacement boundary conditions are applied to make the results converge. The model is shown in Figure 5.

<figure><img src="image_6.png" /><figcaption>Figure 5 Compression/traction unit</figcaption></figure>

On the coordinate of $||s|| - s_m$, the D-P model constitutive curve is obtained and plotted in Figure 6.

<figure><img src="image_7.png" /><figcaption>Figure 6 D-P model $||s|| - s_m$ curve</figcaption></figure>

As compared with Figure 1, an excellent match of D-P model constitutive curve can be seen by simulation and analysis. With different dilatancy angles, i.e. $\theta = \phi = 40^\circ$, $\theta = \phi / 2 = 20^\circ$, $\theta = \phi / 4 = 10^\circ$ and $\theta = 0^\circ$, we can simulate four curves under compression and traction respectively as follows,

<figure><img src="image_8.png" /><figcaption>Figure 7 $u_z - \sigma$ curve of compression/traction</figcaption></figure>

We observe from Figure 7 that the limit stress increases with the increase of dilatancy angle. The $u_z - \sigma$ curve

conforms to the ideal plastic. When the dilatancy angle does not equal to the internal friction, softening phenomenon appears at the beginning of the loading in traction case. The greater gap between dilatancy angle and internal friction is, the more obvious softening phenomenon will be.

### 2) Analytical solution

For this simple example, an analytical limit stress is available and given by [10]

$$ \sigma = \frac{c}{\tan \varphi} \frac{2}{1 \pm \left(2 - \tau \xi / \sqrt{3\xi(2 - \tau^2\xi)}\right)} \quad (53) $$

$$ \text{with } \tau = \tan \theta / \tan \varphi \quad \text{and} \quad \xi = 1 / (3 \cdot k_d^2 \tan^2 \varphi) \quad (54) $$

where “-” stands for compression, “+” stands for traction. The accuracy of the simulation is verified by comparing numerical solution and analytical solution in Table 1.

<table><caption>Table 1 Numerical solution and analytical solution</caption><thead><tr><th rowspan="2">θ (°)</th><th colspan="2">Analytical solution<br>(MPa)</th><th colspan="2">Numerical solution<br>(MPa)</th></tr><tr><th>compression</th><th>traction</th><th>compression</th><th>traction</th></tr></thead><tbody><tr><td>40</td><td>-128.669</td><td>27.97837</td><td>-128.669</td><td>27.9784</td></tr><tr><td>20</td><td>-121.094</td><td>27.60294</td><td>-121.094</td><td>27.6029</td></tr><tr><td>10</td><td>-115.046</td><td>27.27608</td><td>-115.046</td><td>27.2761</td></tr><tr><td>0</td><td>-108.439</td><td>26.88765</td><td>-108.439</td><td>26.8876</td></tr></tbody></table>

## B. Plane strain problem

The second example is devoted to apply the D-P model to a structure under plane strain conditions. The material characteristics are the same as in the above example. The geometry parameter and the finite element mesh are shown in Figure 8. The mesh involves 128 rectangle elements.

<figure><img src="image_9.png" /><figcaption>Figure 8 Geometry parameter and mesh</figcaption></figure>

Assuming that the soil is weightless, given a displacement of $u_z = -0.001m$ on the right upper corner, we can simulate $u_z$, $\sigma_z$ and $\varepsilon_z$ for $\theta = \phi = 40^\circ$. The result is displayed by the in-house post-processing software FERView.

<figure><img src="image_10.png" /><figcaption>Figure 9 Displacement field</figcaption></figure>
The displacement fields (Figure 9) show the continuity of the deformation. Due to the compressibility of the soil, an obvious phenomenon of compression would appear when enlarging the magnification factor.

<figure><img src="image_11.png" /><figcaption>Figure 10 Stress fields</figcaption></figure>

The maximum stress in Figure 10 is greater than the limit stress in the first example. The reason can be explained by the ferruling phenomenon under the three-side fixed constraints. The compressive stress is mainly concentrated in the areas of displacement constraints.

<figure><img src="image_12.png" /><figcaption>Figure 11 Plastic strain fields</figcaption></figure>

As shown in Figure 11, the plastic zones can be clearly illustrated. In addition, it is only in the places where large shape changes that exists the zones in plastic stage. Other zones are still in elasticity.

Besides, more results under different dilatancy angles can be obtained as indicated in Table 2.

<table><caption>Table 2 Results of different dilatancy angles</caption><thead><tr><th rowspan="2">θ (°)</th><th colspan="2">Displacement (m)</th><th colspan="2">Stress (MPa)</th><th colspan="2">Plastic strain</th></tr><tr><th>max</th><th>min</th><th>max</th><th>min</th><th>max</th><th>min</th></tr></thead><tbody><tr><td>40</td><td>5.294e-5</td><td>-0.001</td><td>1.847</td><td>-157.4</td><td>4.533e-4</td><td>≈0</td></tr><tr><td>20</td><td>4.936e-5</td><td>-0.001</td><td>1.902</td><td>-155.5</td><td>6.131e-4</td><td>≈0</td></tr><tr><td>0</td><td>4.438e-5</td><td>-0.001</td><td>1.981</td><td>-153.5</td><td>8.363e-4</td><td>≈0</td></tr></tbody></table>

According to Table 2, a rule can be summarized: the maximum stress at the loading direction enlarges with the increase of the dilatancy angle. The plastic zones are also expanding with the dilatancy angle. The simulation results almost match with that of ANSYS code. The rationality of the bi-potential theory and the feasibility of its algorithm can thus be proved.

## VI. Conclusion

The bi-potential of the D-P model for non-associated soil materials is established. It has been shown that the bi-potential theory allows expressing better the flow rule of non-associated materials. By the finite element discretization, a constitutive integration method for the D-P model is obtained. Numerical solutions are validated as compared with theoretical ones. Through the simulation of the plane strain problem, the

rationality of the bi-potential theory and the feasibility of its algorithm are well verified. To conclude, the bi-potential theory in dealing with set problems shows its great potential in future research. Considering large deformation and frictional contact interfaces, this approach can be extended to deal with a rigid wall standing on the ground. The fields of the soils under the influence of the rigid wall could be analyzed. What's more, the simulation of the pile driving, a still complicated problem, could be analyzed when considering dynamic effects.

## References

[1] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, vol. 2, 4th Edition, McGraw-Hill, London, 1991.

[2] J.C. Simo, T.J.R. Hughes, Computational Inelasticity. Springer (1998).

[3] G. de Saxcé, A generalization of Fenchel's inequality and its applications to the constitutive law. C.R.Acad, Sci, Paris 314, Serie 125-II, 129(1992)

[4] A. El Hanafi, J. Chaoufi, C. Vallée, A. Germaneau, K. Atchonouglo, H. Fatmaoui, A. Ghafiri. Construction of a bipotential representing a linear non-associated constitutive law. Comptes Rendus Mecanique 341 (2013) 667-671

[5] A. Berga. Mathematical and numerical modeling of the non-associated plasticity of soils—Part 1: The boundary value problem. International Journal of Non-Linear Mechanics,47(2012) 26-35

[6] C. Bouby, G. de Saxcé, J.B. Tritsch. A comparison between analytical calculations of the shakedown load by the bipotential approach and step-by-step computations for elastoplastic materials with nonlinear kinematic hardening. International Journal of Solids and Structures 43 (2006) 2670-2692

[7] D.C. Drucker, W. Prager, Soil mechanics and plasticity analysis or limit design, Quarterly Applied Mathematics. 10(2) (1953) 157-165.

[8] M. Hjiaj, J. Fortin, G. de Saxcé. A complete stress update algorithm for the non-associated Drucker-Prager model including treatment of the apex. International Journal of Engineering Science 41 (2003) 1109-1143

[9] M. Hjiaj, D.-L. Dao, G. de Saxcé. A family of bi-potentials describing the non-associated flow rule of pressure-dependent plastic models. Acta Mech 220, 237-246 (2011)

[10] A. Berga. Mathematical and numerical modeling of the non-associated plasticity of soils—Part 2: Finite element analysis. International Journal of Non-Linear Mechanics, 47(2012) 36-45