## 4 Consistent Tangent Stiffness
In the following examples, the strain history will be imposed by defining the strain increment vector \(\Delta\mathbf{\epsilon}\) as:

\[\Delta\mathbf{\epsilon}=\overline{\Delta\epsilon}_{V}[\mathbf{n}_{1}+c_{1}\mathbf{ n}_{2}\ ] \tag{61}\]

where \(\overline{\Delta\epsilon}_{V}\) is the desired volumetric strain (to be imposed), \(\mathbf{n}_{1}=[1\,1\,1\,0\,0\,0]^{T}\) (using Voigt's notation for the stress tensor), \(c_{1}=\frac{K\tan\beta}{2G}\), and \(\beta\) defines the relationship between \(p\) and \(q\) (\(\Delta p=\tan\beta\,\Delta q\)) in the first increment. The vector \(\mathbf{n}_{2}\) defines which stress directions were activated, and it varies from case to case. All the analyses were initialized with a non-zero hydrostatic stress through the parameter \(p_{0}\). \begin{table}
\begin{tabular}{c c c c c c c} \hline \hline \(\mathbf{\kappa}^{*}\) & \(\mathbf{\lambda}^{*}\) & \(\mathbf{\mu}^{*}\) & \(\mathbf{M}\) & \(\mathbf{\
u}\) & **OCR** & \(\mathbf{\tau}\) **(days)** \\ \hline
0.0084 & 0.061 & 0.0011 & 1.33 & 0.3 & 1.39 & 1 \\ \hline \hline \end{tabular}
\end{table} 
Table 1: Material parameters for the Vermeer–Nev (V-N) law. OCR, overconsolidation ratio. #### 5.1.1 Test 1: Monotonic Loading

In the first test, an increasing strain was imposed at a single Gauss point. Because of the imposed strain path (61), the obtained stress had both a hydrostatic and a deviatoric component. The imposed deformation increment was chosen to activate the stress component in Direction 11, i.e., aligned with the vector \(\mathbf{n}_{2}=[1\,0\,0\,0\,0\,0]^{T}\). The other parameters were: \(\overline{\Delta\epsilon}_{V}=10^{-4}\), \(\beta=30^{o}\), and \(p_{0}=10^{5}\) Pa. The evolution in time of the stress component \(\sigma_{11}\) is depicted in Figure 3. To be precise, the stress is plotted as a function of the time step number, where time steps of constant amplitude were used. The curves show an excellent agreement between implicit and explicit integration using a time step \(\Delta t=0.1\), where \(\Delta t\) is the non-dimensional time step size \(\Delta t=(t^{n+1}-t^{n})/\tau\) and the time units of \(\tau\) are years. It can be observed that the implicit scheme was able to provide an accurate solution even when using a time step 10 times (\(\Delta t=1.0\)) larger than the explicit case. When using the same larger time step (\(\Delta t=1.0\)), the implementation with explicit integration was not able to provide a solution. To confirm the accuracy of the proposed algorithm, an explicit solution with a very small time step (\(\Delta t=0.01\)) is plotted. No significant difference can be appreciated between the curves. The same results are subsequently plotted in a \(p-q\) plane in Figure 4. In this plane as well, it is possible to observe a very good agreement between the implicit solution with large time step and the explicit one obtained with a much smaller time step. A slight difference is visible in the implicit solution with the larger time step, due to the less frequent evaluation of the stress along the (nonlinear) loading history. Figure 4: Test 1. Stress point evolution in the \(p-q\) plane. Figure 3: Test 1. Time evolution of the stress \(\sigma_{11}\) at different time steps for explicit and implicit time integrations. #### 5.1.2 Test 2: Monotonic Loading

The second test was very similar to the previous one, the only difference being the activated stress component. In this second case, the increasing strain activated the stress in Direction 12, i.e., aligned with the vector \(\mathbf{n}_{2}=[0\,0\,0\,1\,0\,0]^{T}\)). Figure 5 shows the time evolution of the stress component for two different time step sizes. For the smaller time step (\(\Delta t=0.1\)), the implicit and explicit results were very close. However, when the time step size increased (\(\Delta t=1.0\)), the explicit solution tended to oscillate.