Two-Step Exact Solution for Modified Mohr-Coulomb Viscoplastic Consistency Model with Linear Hardening/Softening: A Proof with Maple Software

It has been observed numerically that the viscoplastic consistency model by Wang (1997) with a linear yield surface and a linear hardening/softening rule converges, using the standard stress return mapping, with two steps. In this paper this numerical observation is proved analytically using the Maple software. The proof is carried out for the Modified Mohr-Coulomb viscoplastic consistency model in the corner plasticity situation, i.e. when both the Rankine and Mohr-Coulomb criteria are violated.


Introduction
The classical viscoplasticity formulations, namely the Perzyna and Duvaut-Lions models, do not utilize the consistency condition [1]. Consequently, trial stress states violating the yield criterion are not returned to the yield surface defined by the yield criterion. In those models, the viscoplastic strain rate depends on the amount of violation [1]. In contrast, the consistency condition is imposed in the more recent viscoplastic consistency model by Wang [2]. Thereby, the trial stresses violating the yield function are returned to the yield surface as in the classical rate-independent plasticity. The major difference to the rate-independent plasticity is the dependence of the yield function on the rate of the internal variable(s). The basic components of the bi-surface viscoplastic consistency model are: where , i i κ κ  denote the internal variable and its rate, respectively. Moreover, g vp,i is the viscoplastic potential, vp ε  is the viscoplastic strain tensor (vector), and i λ  is the viscoplastic increment, i.e. the rate of viscoplastic multiplier.
The yield function f vp,i can be called a dynamic yield function since it can expand and shrink in relation to its static position depending on the loading rate. Moreover, the loadingunloading conditions of Kuhn-Tucker form must hold.
As a consequence, the robust stress return algorithms of rate-independent computational plasticity can be employed with minor modifications concerning the presence of the rate of the internal variable(s). It is shown in this paper that a two-step exact solution exists for a viscoplasticity problem with a linear bi-surface yield function and linear hardening/softening laws. This fact has, of course, been observed in numerical simulations but a formal proof is given here. The proof is performed for the Modified Mohr-Coulomb (MMC) viscoplastic consistency model in the corner plasticity situation, i.e. when both the Rankine and Mohr-Coulomb criteria are violated. Moreover, the proof is carried out with the Maple technical computing software due to the complexity of the equations involved. The MMC model is chosen due its wide usage in geotechnical engineering analyses. The MC criterion is usually modified with the Rankine, or the principal stress, criterion to predict the correct tensile strength of geomaterials such as rock and concrete.

Modified Mohr-Coulomb Viscoplastic Consistency Model
The MMC yield criterion in the principal stress space, assuming the ordering 1 2 3 σ σ σ ≥ ≥ , consists of the MC and Rankine criteria: where symbol x ∂ denotes the derivative with respect to a variable (scalar or vector) x. The combined yield surface is illustrated in Figure  1. The rate-dependent softening/hardening rules for cohesion and tensile strength are assumed to be linear by where c 0 is the initial (intact) value of cohesion while h MC , h R and s MC , s R are the softening/hardening and viscosity moduli in compression and tension, respectively.
In order to account for correct dilatation behaviour of the material in compression, the plastic potential is chosen as where ψ is the dilatation angle. Due to its linearity, the gradients of the MC yield surface and plastic potential are constant vectors (having the angles φ and ψ as parameters). With this choice for the plastic potential, the expression for . Therefore, since the relation between the internal variables and plastic increments is a constant relation, it is more convenient to use the plastic multipliers and increments in the following derivations.

Return Mapping Scheme for Bisurface Viscoplastic Consistency Model
The intricacy involved in the multisurface plasticity is the stress return to an intersection of yield criteria. At such a point the gradient of the combined yield surface is not unique as illustrated in Figure 2 depicting a corner point in a bi-criteria case.
The set of admissible stresses is denoted by E σ in Figure 1.
, and in . The stress return algorithm should be able to distinguish between these regions and restore the consistency accordingly. This means, e.g., that if the trial stress lies in Γ 1 the single surface integration should be performed in order to return the trial stress to the surface f 1 = 0. The determination of the right return can be performed by trial and error only because the condition trial 0 i f > does not guarantee that eventually Here the procedure presented in [3] is employed for determining the active set of yield surfaces.
Stress integration, i.e. return mapping, algorithms for classical single and multi-surface plasticity and viscoplasticity are considered, e.g. in [1,3], while algorithms for consistent viscoplasticity are developed in [2,4,5]. For present purposes, the stress integration algorithm by Wang et al. [2] is utilized as modified for the linear bi-surface case. It is assumed that both criteria are violated.

Return mapping algorithm for bi-surface viscoplastic consistency model:
Compute trial state: ∆ = ε 0 and go to 1. Local iteration: The purpose of step 3 in the algorithm is to ensure that a genuine corner plasticity situation is realized. This can be performed on trial and error only, as noted above. If one of the cumulative plastic increments is found non-positive, the trial stress was located in one of the regions denoted by Γ 1 and Γ 2 in Figure 2 and, thus, no genuine corner plasticity has occurred (i.e. the trial stress was not located in region Γ 12 ).
This means that the bi-surface iteration must be terminated and a single surface integration should be performed to return the trial stress onto the correct yield surface as explained above.
In fact, in this form the algorithm above is similar to the generalized cutting plane algorithm [3]. Therefore, it is not limited to linear yield surfaces but can be used for stress integration of nonlinear models as well.

A Two-Step Solution for the MMC Viscoplasticity Problem
It is convenient for present purposes to write the MC criterion with its viscoplastic potential in the following simpler form: where f c is the uniaxial compressive strength which is substituted for c in Equation (3). In addition, the constant k MC becomes now 2 MC Next, it is shown by Maple software that a two-step exact solution, using the algorithm above, exists for the MMC viscoplasticity problem in the corner plasticity situation with the linear hardening/softening rules ( , is treated as shown in Figure 3.  Therefore, at the onset of viscoplasticity, a single step solution exists. Before continuing, the above Maple script is elaborated in some detail. After loading the linear algebra package in line 2, the linear elasticity matrix is defined in line 3 with E 1 = E(1−ν)/(1+ν)(1−2ν) and E 2 = Eν/(1+ν)(1−2ν) Lines 4-8 define the trial values of the MC and Rankine criteria, the plastic potential (for MC criterion), and their gradients determined with the diff command. In lines 9 and 10 the initial values are set for the viscoplastic multipliers and strain. Vector F and matrix G according to Equation (5)  . This is a general time step involving softening/hardening and the strain rate effects as well. The Maple code is shown in Figure 4.
Thus, in the end of the first iteration the values of the yield functions are given by (7). Accordingly, the exact solution is obtained with the first step if the viscosity moduli s MC and s R are set to zero which is the rate-independent case. In the rate-dependent (viscous) case another step is needed: