## 4 Consistent Tangent Stiffness
The same results can be also observed in a \(p-q\) plane in Figure 6. As expected, this test confirmed that the implicit integration scheme was more robust and allowed for the use of a larger time step size, possibly reducing the global computing time of the analysis. #### 5.1.3 Test 3: Loading and Unloading

In this test, a loading and unloading evolution of the strain was imposed. First, it increased linearly (\(\overline{\Delta e}_{V}=10^{-4}\) for 10 time steps), then it remained constant for a while (\(\overline{\Delta e}_{V}=0\) for 10 time steps); finally, it decreased linearly again to zero (\(\overline{\Delta e}_{V}=-10^{-4}\) for 10 time steps). Due to the viscous behavior, in the region with constant strain, a stress relaxation took place. The other parameters were: \(\mathbf{n}_{2}=[1\,0\,0\,0\,0\,0]^{T}\), \(\beta=30^{o}\), and \(p_{0}=5\cdot 10^{5}\) Pa. In Figure 7, results in the \(p-q\) plane are shown. Figure 5: Test 2. Time evolution of \(\sigma_{11}\) stress at different time steps for explicit and implicit time integrations. Figure 6: Test 2. Stress point evolution in the \(p-q\) plane. As in the previous cases, implicit and explicit integration schemes gave very close results when a relative small time step was used (\(\Delta t=0.1\)). On the contrary, when a larger time step was used (\(\Delta t=1.0\)), only the proposed implicit scheme could provide a non-oscillating solution. ### On the Validation of the Consistent Stiffness Matrix in a Standard FE Code

To verify the implementation of the consistent tangent stiffness matrix (60) in a commercial FE code, the following procedure was followed. The implicit time integration of the constitutive model was first carried out over a time step of a regular size for a prescribed strain increment. Stresses \(\sigma^{n+1}\) and viscoplastic strains \(\epsilon^{vp\,n+1}\) at the end of the step were therefore computed as an output of the implicit time integration. Starting from this computed stress state, six very small (of the order of \(10^{-8}\)) strain increment vectors \(\delta e_{i}\), \(i=1,\ldots 6\) were assigned. Each strain increment vector \(\delta e_{i}\) consisted of a very small non-zero \(i\)-th component, while the remaining five components were zero. The corresponding vector \(\delta(\sigma_{\text{loc}})_{i}\) was then computed by solving the local backward difference problem stated through the nonlinear system (35) and the update (21). This stress vector was compared to the estimate obtained as \(\delta\big{(}\sigma_{\text{approx}}\big{)}_{i}=\left[\frac{\partial\sigma}{ \partial\epsilon}\right]^{n+1}\delta e_{i}\), where \(\left[\frac{\partial\sigma}{\partial\epsilon}\right]^{n+1}\) is the tangent stiffness matrix (60) at time \(t^{n+1}\), in correspondence with the computed stresses \(\sigma^{n+1}\). This estimate is more acceptable the smaller the strain increment \(\delta\epsilon_{i}\) is. For each of the six strain increment vectors \(\delta\epsilon_{i}\), the accuracy of the \(i\)-th column of the tangent stiffness matrix was tested. The tests described above had to be repeated computing the tangent stiffness matrix for several different initial stress states \(\sigma^{n+1}\). As an example, for an initial hydrostatic stress state \(p_{0}=5\cdot 10^{5}\) Pa, the verification of \([\partial\sigma/\partial\epsilon]^{n+1}\) was carried out at six different stress values \(\sigma^{n+1}=p_{0}\mathbf{m}+\Delta\sigma\) (with \(\Delta\sigma\) a function of \(\Delta\epsilon=\mathbf{n}_{2}\,\Delta\epsilon\), obtained by varying \(\mathbf{n}_{2}\) with one non-zero entry at a time and \(\Delta\epsilon=10^{-4}\)). The result, in terms of the \(\delta\sigma_{11}\) component, is shown in Figure 8, evidencing in all cases a small, but unavoidable difference between \(\delta(\sigma_{\text{loc}})_{11}\) and \(\delta(\sigma_{\text{approx}})_{11}\).