Memory repositioning in soil plasticity models used in contact problems

We aim to enhance the stability of finite element models of dynamic structural contact with multiphase granular soils which are described by advanced soil plasticity models that can simulate monotonic and cyclic behaviour of multiphase soils. Often, numerical oscillations cannot be avoided in these contact models and can cause advanced soil models to significantly overshoot stress, leading to unrealistic discontinuities in the stress paths. This situation can challenge the stability of the stress integration scheme and the global finite element solver and lead to the early termination of the analysis. We specifically address the issue of stress overshooting by presenting novel solutions and the corresponding stress integration schemes for a representative soil model for unsaturated granular soils. Also, several examples are provided to evaluate the integration scheme and show the advantages and limitations of the proposed overshooting solutions in solving a contact-impact problem involving unsaturated granular soils.


Introduction
Modelling dynamic structural contact with multiphase granular soils (i.e., saturated and unsaturated soils) has many applications such as pile installation, ground improvement [1] and ground assessment using dynamic indentation. However, this class of problems is computationally challenging since the accurate representations of the interface, the contributions of suction and water pressure to the elastoplastic soil behaviour, and drainage conditions generally demand the adoption of "fully coupled models" [2][3][4][5][6][7][8] which are nonlinear and cannot generally guarantee unconditional stability. Therefore, novel solutions to improve the stability of such nonlinear models are needed as geotechnical applications of B Javad Ghorbani Javad.ghorbani@monash.edu 1 ARC Smart Pavements Hub, Department of Civil Engineering, Monash University, Clayton, VIC 3800, Australia 2 School of Engineering, The University of Newcastle, Callaghan, NSW 2308, Australia 3 University of California, San Diego, La Jolla, CA 92093-0085, USA these contact models are broad and elastoplastic soil models are becoming more complex.
This paper presents a novel solution to improve the stability of numerical models of structural contact with multiphase granular soils which are mathematically described by multisurface plasticity models (e.g., [9]) that can simulate cyclic and monotonic behaviour of granular soils. These plasticity models generally use an algorithm to automatically detect the occurrence of a load reversal. To capture the experimentally observed increase in the soil stiffness in the initial stages of a load reversal during cyclic loading [10,11], the plastic modulus is often reset to infinity (i.e., a very large value in the context of numerical analyses) whereby a stiff pseudoelastic behaviour is predicted [12]. Nonetheless, if spurious numerical oscillations (instead of true loading reversals) trigger a tiny reversal event that is immediately followed by reloading, the plasticity model can unrealistically become stiff and significantly overshoot stress. The occurrence of overshooting is an undesirable situation in numerical simulations that can challenge the overall stability of the numerical solutions. In dynamic and coupled contact problems, numerical oscillations generally cannot be avoided, and therefore, early termination of the analyses due to stress overshooting can become a major problem demanding appropriate treatment.
Several solutions have been proposed to resolve the overshooting issue in advanced plasticity models for saturated soils. In [13] using a "virtual bounding surface" for solving the overshooting problem was proposed; a solution that appears to be only applicable to the multilayer model for saturated soils proposed by the authors. The authors also developed an implicit integration scheme that included their proposed strategy. An alternative approach for solving the overshooting problem is the development of measures to automatically prevent the plastic modulus from drastic changes when spurious oscillations trigger the reversal [9,[14][15][16][17][18][19]. In this approach, a criterion is defined to automatically evaluate if the reversal is a true reversal or an unrealistic reversal that is induced by spurious oscillations. Such an evaluation is often made by checking if the magnitude of the plastic deviatoric strain during the load reversal is above or below a predefined threshold [15].
To the best of the authors' knowledge, there are no solutions for the overshooting problem in advanced unsaturated soil models that can simulate cyclic behaviour. Experimental results suggest a significant impact of suction and the degree of saturation on cyclic and monotonic volume change behaviour of unsaturated granular soils (e.g., [20][21][22]). To capture this impact in unsaturated soils, different from saturated soil models, key aspects of plasticity of unsaturated soils such as the hardening law, flow rule, and plastic modulus are often defined as a function of suction and/or the degree of saturation. In analyses of dynamic structural contact with unsaturated soils, effective drainage of pore water and air pressure may not be possible during the analyses. In such conditions, the unsaturated soil model can be prone to stress overshooting since unrealistic reversal events can be triggered by spurious oscillations of suction/air pressure, pore water pressure, the degree of saturation, and the effective stress (in addition to displacement). Potential solutions to overshooting problems in unsaturated soil models can be developed by generalising existing solutions which were developed specifically for saturated soil models. Nonetheless, such potential solutions need to be rigorously tested against spurious oscillations in factors such as suction and degree of saturation, which are exclusive to the analyses of unsaturated soils. Also, robust stress integration schemes with these solutions need to be developed to facilitate the implementation of advanced unsaturated soil models to finite element packages. This study will address these knowledge gaps and propose two solutions to mitigate the overshooting problem in a representative multisurface plasticity model known as MUD (Model for Unsaturated soil Dynamics) [23]. The performance of the integration scheme and the overshooting solutions will be evaluated in a series of analyses of a contact-impact problem involving multiphase granular soils.
We will begin with a short introduction to the representative constitutive model and then discuss the issue of overshooting by giving several numerical examples. The examples highlight the limitations of an existing approach to dealing with this issue. To resolve these limitations, we formulate a novel approach to effectively prevent stress overshooting and present the associated numerical integration strategy. In a subsequent section, boundary value problems are solved to evaluate the efficiency of the proposed scheme in solving frictionless contact problems involving multiphase soils.

General definition
The soil model MUD is a multisurface model for capturing the cyclic and monotonic behaviour of saturated and unsaturated granular soils. This model establishes a direct relationship between the changes in the stress ratio and the kinematic hardening parameter by defining a linear loading surface in the space defined by the mean effective stress and deviatoric stress. Through such a definition, the model can eliminate the need to enforce the consistency conditions at the end of the integration process; a feature that simplifies the stress integration algorithm [24]. The model also eliminates the purely elastic region by defining a "loading surface" [25] that always passes through the current stress state which further simplifies the explicit stress integration algorithm as it allows the removal of the algorithms that are needed for calculating the plastic portion of the applied strain increment in each time step of the analysis. The model can simulate soil behaviour under a wide range of hydro-mechanical loads by considering the impact of stress-induced anisotropy on hardening and flow rules by employing Bishop's effective stress. It should be noted that the definition of the Bishop's effective stress parameter is not unique and some notable suggestions for this parameter are the degree of saturation, [26], an "effective degree of saturation" [27,28], and a function of suction [29]. For granular soils with a negligible residual degree of saturation in a dry state, the first two suggestions can be considered nearly the same where the Bishop's parameter is equal to the degree of saturation [26,[30][31][32]. Nonetheless, often a combination of Bishop's effective stress and an additional strain-like quantity such as suction, p c is used to more accurately describe hardening/softening induced by the changes in the degree of saturation (or suction) [23,31,33,34]. In MUD, the Bishop's parameter is taken as the degree of saturation, S w [26] and the strain-like quantity is the "bonding variable", ξ [35]. The parameter ξ is a scalar quantity and denotes the average intergranular force generated by the presence of water menisci in the soil skeleton and is defined as where p c is a normalised suction (with respect to the atmospheric pressure) and g( p c ) describes the intergranular force exerted by the water meniscus for which the hyperbolic formula given in [30] is used, i.e., g(p c ) = 1 +p c 10.7+2.4p c . The model is formulated within the critical state theory and uses the concept of the "state parameter" [36] that enables a unified description of the response of granular soils at different initial void ratios. The Soil Water Retention Curve (SWRC) relates the matric suction to the degree of saturation, and the shape is a function of the pore size distribution and voids ratio [37,38]. The critical state line is a generalisation of the model proposed in [39] by incorporating a bonding variable. In addition, MUD employs a combination of kinematic and isotropic hardening laws to simulate different stress paths. Stress ratio changes are macroscopically modelled by a kinematic hardening law. The evolution of isotropic hardening parameter, α iso is defined by using the concept of the Limiting Compression Curve (LCC) which is considered as a straight line in the bi-logarithmic space defined by ln e − ln p following [40] and is active under constant stress ratios and close to the LCC.
The incremental relationship between the effective stress, σ , strain ε, and ξ is written as [23] dσ = D ep dε + S ep dξ (2) where D ep is the elastoplastic matrix describing the relationship between stress and strain, and S ep is the elastoplastic matrix that relates the rate of change in the effective stress to the bonding variable. It is important to note that upon reaching the state of saturation, the Bishop's effective stress will smoothly converge to the Terzaghi's effective stress. In such a condition, key features of MUD including the kinematic hardening rule and dilatancy will degenerate to those used in the "SANISAND" family of models proposed by Dafalias and coworkers (e.g., [11,41]). Aspects of the model formulations such as the definitions of D ep and S ep as well as dilatancy and the plastic multiplier,λ are presented in Appendix A.

Stress overshooting
The kinematic hardening law and the definition of the plastic modulus are defined in this section to explain why stress overshooting occurs in the model. Similar to many multisurface soil models that are developed for modelling cyclic and monotonic behaviour of soils, in MUD, the plastic modulus is reset to infinity at the initiation of a new stress reversal. Such a resetting mechanism, which is based on the approach proposed in [11] (to be discussed in this section), results in stiff and pseudo-elastic behaviour immediately after the start of a new loading direction. When trivial changes in the loading direction happen (e.g., due to numerical oscillations rather than a true change in the loading), the resulting changes in the plastic modulus can be drastic and cause an unrealistic discontinuity in the stress-strain path if a tiny reversal event is followed by reloading. In this situation, significant overestimation (or underestimation) of the predicted stress values and early termination of numerical analyses can occur. Therefore, for a robust performance of the model in problems prone to oscillatory stresses, such as dynamic contact with multiphase soils, proper treatment of the issue is needed with particular attention to improving the hardening law and the definition of the plastic modulus as shown in [15].
In this section, we will discuss the approach of resetting the plastic modulus in MUD at the initiation of a load reversal and show the numerical issues that can arise as a result of using this approach. To this end, we employ Voigt notation to describe the model in multiaxial space. We define the stress ratio vector as follows where p is the mean effective stress and s is the vector of deviatoric stress defined by s = σ − p m with σ being the effective stress vector and m T = 1 1 1 0 0 0 . We consider that the critical state ratio, M can generally vary between its values obtained in triaxial compression and extension. To formulate the changes, a Lode angle, θ , is defined in the direction of η σ − α k (with α k as the kinematic hardening vector). We use the equation proposed in [42] as follows where α = M e M c (with M e and M c respectively denoting the values of M in triaxial extension and compression tests). Also, we define n as follows Based on this definition, cos3θ = √ 6 n x x 3 + n yy 3 + n zz 3 where the terms n α (α = xx, yy, zz) representing the first three components of vector n.
Also, we define α β as the multiaxial counterpart of a given triaxial quantity α β as follows α β = 2 3 α β n. Given the foregoing descriptions, and following [11], in MUD, the kinematic hardening vector, α k can evolve as follows where α b = 2 3 α b n and α b (or the triaxial counterpart of α b ) connects the hardening law to the state parameter, ψ [36], as follows where m iso is a small positive quantity (≈ 0.05M c following [43]) and n b is a material parameter and denotes the MacCauley brackets. The parameter h in Eq. (6) is defined by where p atm is the atmospheric pressure, and h 0 is a material parameter, and c h regulates the contribution of the void ratio, e in the hardening law. The value of α k upon the initiation of a load reversal will be stored in α in which hereafter is referred to as the memory of load reversal. The rule for updating is as follows. If upon the initiation of loading stepi where the superscript (i) denotes loading step i. For the cases where the stress ratio changes (for whi the kinematic hardening is used), the plastic modulus, K p can be described by Equations (8) to (10) indicate how overshooting may be triggered during numerical oscillations. From these equations, it can be inferred that if numerical oscillations are interpreted as true stress reversals, we will get (α k − α in ) T n = 0 which results in K p approaching infinity, leading to a potential for a discontinuity in the stress-strain path.
An example of stress overshooting imposed by numerical oscillations is shown in Fig. 1. An isotropically consolidated sample of saturated Toyoura sand with the initial mean effective stress of 100 kPa and void ratio of 0.735 was sheared under undrained and triaxial conditions. The experimental data was reported in [44] and the test was simulated by using strain increments of magnitude 0.00001. The parameters of the mechanical model were taken from [23] and are summarised in Appendix B. It is also notable that the adaptive explicit stress integration scheme developed in [23] was used in this simulation, and the stress integration tolerance was set to 10 −4 in this example.
The grey dotted line shows the stress-strain path if no oscillation is induced in the analysis. Note that the q axis represents the deviatoric stress in this graph. The green dashed line shows the consequence of inducing a small oscillation in the stress-strain path if no treatment is applied. This path is simulated by applying a numerical oscillation in the form of a small load reversal event that commences upon reaching an axial strain of 0.08. The reversal step is simulated by applying an axial strain increment of -0.00017. It is seen that such a small increment could result in a drastic change in the predicted deviatoric stress because a stiffer plastic modulus is adopted in the subsequent reloading step.
The curves in Fig. 2 demonstrate the evolutions of the term (α k − α in ) T n during the simulations of this test and a test with no induced oscillation. The graph shows discontinuities in the predicted values of (α k − α in ) T n near the point of reversal and during the subsequent reloading path after the reversal. When the induced reversal occurs, initially we obtain (α k − α in ) T n < 0. However, since negative values are not acceptable, (α k − α in ) T n will be set to zero following Eq. (9), leading to a very large value of plastic modulus per Eqs. (8) and (10) and the resulting discontinuities in the stress-strain path seen in Fig. 1.
On the same graph, the performance of a remedy for the stress overshooting effect is demonstrated. This remedy follows the proposed method in [15] by defining a "threshold" for the equivalent deviatoric plastic strain during a loading process, ε p q . According to this method, α in is updated as follows: where α r in represented the value of α in that is stored in the previous loading step (i.e., step (i − 1)). Also, the weighting coefficient m q is computed according to: where j is a parameter that is set to a default value of 1 [15] and e p q is the deviatoric plastic strain vector during the load reversal in step (i) of the analysis. It should be noted that, compared with the updating rule outlined in Eq. (9), the term m q may be viewed as a repositioning term where the memory of the load reversal, α in , is repositioned from the position α to mitigate the stress overshooting problem. It may also be noted that if the developed deviatoric plastic strain during the load reversal is way below the threshold,m q → 1, leading In other words, if the unloading event is trivial, it will not significantly affect the memory of the load reversal and the plastic modulus in the subsequent reloading process in step (i + 1). Nonetheless, in unloading events where m q approaches zero, the repositioning scheme will degenerate to the updating rule outlined in Eq. (9), resulting in an infinite plastic modulus at the start of the reloading process in step (i + 1). This approach requires storing α r in , which is a vector, in addition to the kinematic hardening vector at the instant of stress reversal.
A robust integration scheme is developed to accommodate all the repositioning strategies discussed in this study in an explicit integration scheme with automatic error control and substepping. To this end, we have modified the integration scheme proposed in [23] for MUD. This approach is an extension of the explicit integration scheme in [45] to fit in with the changes induced by the presence of matric suction as a strain-like quantity in the model formulations and the absence of a purely elastic region in the model.
The integration scheme divides the strain increment into several substeps. In each sub-step initiated at a pseudo-time oft, with 0 ≤t ≤ 1, the integration error is evaluated using the difference between a second-order accurate modified Euler solution and a first-order accurate Euler solution. A prescribed stress integration tolerance denoted by STOL is used to control the magnitude of the integration error. The following algorithms describe the process of the development of a stress integration scheme that includes the repositioning strategy proposed in [15]. The developed integration scheme also presents a refined error calculation to accommodate the needed calculation of the plastic deviatoric strain during reversal steps.
An important part of the algorithm is a state variable termed "Lstepcount". This state variable is updated each time the load is reversed. At the start of the analysis, it possesses a value of zero and is updated according to the following mechanism. Upon the initiation of load reversal, which is identified by checking if (α k − α in ) T n < 0, Lstepcount needs to be updated to Lstepcount + 1. If Lstepcount is an odd number, it will be viewed similar to step (i) of the load sequence discussed in Eq. (12). In such steps, it is essential to store e p q per Eq. (12). However, a better strategy is to store ε p q = 2 3 e p q T e p q instead, since ε p q is a scaler quantity requiring smaller storage, and will be directly used in Eq. (12). If Lstepcount is an even number, the algorithm will perform repositioning treatments similar to that outlined in Eq. (12) for step (i + 1). The following sub-algorithm (Algorithm 1) is proposed for the foregoing updating mechanism.
By using this sub-algorithm, the whole integration scheme using the substepping technique is presented in Algorithm 2. In this algorithm ε and z represent strain and fabric vectors, respectively and ε The repositioning scheme proposed in [15] will effectively prevent overshooting in many cases where after the update of α 1 in using Algorithm 1, we get (α k − α in ) T n > 0, a necessity to effectively prevent drastic change of the plastic modulus if oscillation occurs. Nonetheless, this algorithm cannot generally guarantee that after repositioning, we obtain (α k − α in ) T n > 0 as highlighted in [15].
More recently, in [17] this problem was addressed by defining an additional criterion to decide if a load reversal is "formal" or "informal". In the latter, α 1 in will experience no update. All reversals are initially assumed informal unless proven otherwise according to the following criterion. If the difference between the accumulated stress ratio and α 1 in is larger than a prescribed tolerance, the reversal is considered formal and repositioning will be performed according to Eq. (11). In this study, we use a different approach that, unlike the scheme proposed in [17], does not require an additional check of reversal events and can effectively address the limitation of the repositioning scheme proposed in [15].
Before presenting our approach to solving the overshooting problem, it should be noted that if after repositioning we get (α k − α in ) T n > 0, the repositioning scheme must be abandoned and the model will degenerate to Eq. (9). In this regard, Algorithm 1 needs to be modified to Algorithm 3.
Having explained the problem with the previous repositioning scheme, we propose a new repositioning approach that guarantees (α k − α in ) T n > 0 after the repositioning process, a potential cause of the ineffectiveness of the scheme outlined in Eq. (11). To this end, by defining J = (α k − α in ) T n 1 , we assume where J r is a state variable that stores the last positive value of (α k − α in ) T n before the current loading process is reversed. In Eq. (16), if a true loading event occurs (i.e., ε p q ≥ ε p q ), m q will be zero, leading to plastic modulus approaching infinity and a pseudo-elastic response similar to the original mechanism in MUD. However, for trivial load reversals (e.g., numerical oscillations), the remedy will define J (i) 1 as a continuous function of J r , preventing (α k − α in ) T n = 0 and the resulting drastic change of the plastic modulus.
After some manipulations of Eq. (16), we can rewrite the rule for updating the reversal stress memory as follows where the term −J r 1 − ε p q ε p q j n (i) can be considered as the repositioning term in this equation. Note that since J r is always positive for trivial load reversals (e.g., numerical oscillations), the proposed repositioning scheme cannot lead to (α k − α in ) T n > 0 after repositioning thus resolving the issue with the repositioning scheme proposed in [15] discussed previously. Algorithm 4 shows the repositioning strategy proposed in this paper.
The term J r can be calculated during the stress integration scheme per Algorithm 2 by modifying step 10 of the algorithm per Algorithm 5. Figure 1 also shows that the proposed repositioning scheme successfully prevented overshooting. A comparison of the number of failed and successful substeps in the two repositioning schemes in that example is given in Table 1. In this table, E q is obtained by the following equation where q re f is a reference deviatoric stress obtained from the uninterrupted stress path of the same test (i.e., the path without any oscillations) using a tight stress integration tolerance of 1E-12.
In the performed series of tests, the total number of substeps and the number of failed substeps increase as STOL increases whereas the relative error decreases. Also, while for STOL = 1E-4 the number of total substeps and the number of failed substeps are almost the same in both repositioning schemes, these numbers become consistently smaller when the proposed repositioning scheme is used with tighter values of STOL. Moreover, the relative errors are smaller in the analyses with the proposed repositioning scheme.
It may also be noted that the scheme proposed in [15] and approaches similar to this scheme, such as that suggested in [17,19], demand storing and additional vectors such as α r in . However, the proposed scheme stores J r which is a scalar quantity, therefore reducing the computational storage and the number of state variables in the analysis. Such a reduction in the state variables can be particularly advantageous in fully coupled problems and unsaturated soil dynamics that often require a large number of state variables due to factors such as suction, the degree of saturation, and hydraulic hysteresis [37].
Both repositioning schemes have shown promising potential to overcome stress overshooting. Nonetheless, a scenario should be mentioned where both schemes may fail to tackle the overshooting problem. It may be noted that to effectively overcome the overshooting issue, the load sequence must match that outlined in describing Eq. (11) where following loading at step (i − 1), we will have a reversal at step (i) that will be followed by reloading at step (i + 1) as schematically shown in Fig. 3a.
Nonetheless, the load at step (i) can be too small so that the reversal of loads in transition from step (i − 1) to (i) cannot be detected by the algorithm (i.e., (α k − α in ) T n > 0) which (among other factors) can be due to an inappropriate combination of the step size and the stress integration tolerance. However, if the algorithm detects the occurrence of reversal in transition between the load in step (i) and the subsequent reloading at step (i +1) in Fig. 3a, the algorithm will consider step (i) as the continuation of the loading process initiated at step (i − 1). The consequences of such an event are outlined in Fig. 3b. It can be seen that the loading step (i + 1) in Fig. 3a is now wrongly deemed as a reversal with the likely consequence of ε p q being much greater thanε p q , leading to a potential for ineffectiveness of both schemes in handling the overshooting issue.
To the best of the authors' knowledge, this issue has not been highlighted in any other works on repositioning schemes and it appears that there is no remedy to deal with this issue. Nonetheless, the occurrence of this issue can be considered rare and can be to some extent managed by selecting a different combination of STOL and step sizes.
It should be noted that when the soil is in an unsaturated state, in addition to oscillations in the displacement fields which result in stress overshooting (e.g., similar to that shown in Fig. 1), oscillations in the suction field can also create a potential for overshooting in the model. A key feature of MUD is the establishment of a direct relationship between the changes in the kinematic hardening tensor and stress ratio tensor through defining a linear loading surface [24]. That is the term (α k − α in ) T can be a surrogate for the total stress ratio change during a loading event (before it is reversed) η σ while the term n can be a surrogate for the direction of the stress ratio increment δη σ |η σ | in a given time step of the analysis. In modelling triaxial loading of isotopically consolidated granular soils, α in = 0 and as the loading proceeds, we obtain (α k − α in ) T n > 0, denoting that η σ > 0 and δη σ > 0. However, when in a time step of the analysis, we obtain δη σ < 0, the memory parameter α in will be updated as discussed earlier. Unlike saturated soils, δη σ can change in response to oscillations in the degree of saturation and/or suction (since the mean effective stress can be affected by these two parameters). This dependency can also create a potential for overshooting if the oscillations in suction fields lead to (α k − α in ) T n < 0 through a change δη σ as shown in the following example.
An isotopically consolidated unsaturated soil with an initial void ratio of 0.735, a suction value of 7.68 kPa, and a degree of saturation of 0.85 is first sheared in a drained condition. The sample has an initial mean net stress of 1 kPa to mimic the stress condition that can occur at integration points that are located near a surface with negligible overburden pressure. After shearing the sample by applying an axial strain of 0.01 over 1000 steps, the sample is unloaded by applying an axial strain of −0.01 over 1000 steps. During the unloading process, we imposed a trivial oscillation in the suction with a magnitude of 4% of the prescribed suction at step 1040 of the analysis. This amount was sufficient to lead to (α k − α in ) T n < 0 with the consequence that is shown in Fig. 4.
The results in Fig. 4 that the induced oscillation leads to a discontinuity in the stress-strain path when no overshooting treatment is used. However, both repositioning schemes perform reasonably well in dealing with this overshooting issue.

Finite element analysis
In this section, we perform a more rigorous comparison of the two schemes by solving highly oscillatory contact-impact problems involving multiphase soil and compare the stability and speed of numerical solutions when the two schemes are used. The researchers in [26,46] were among the pioneers who extended Biot's theory of dynamic consolidation and developed a "fully coupled model" to simulate the response of multiphase soils within the finite element framework by assuming immiscible flows and an isothermal environment. While the fundamental assumption of having immiscible flows in soils has remained the same in the majority of the subsequent modifications to this framework, features such as the ability to simulate hydraulic hysteresis [7,37], large deformations [47,48], temperature changes [49,50], the effect of volume changes on SWRC [30,51], and stressinduced anisotropy [24] have been among the improvements applied to this framework. More recently, in [52] the frictionless mortar-type discretisation [53,54] was extended to develop a fully coupled framework for simulating contact with multiphase soils. This approach is used in this study and will be briefly introduced in the following section.
In the mortar algorithm, a segment-to-segment discretisation is formulated based on the projection of the segments on the surface of one of the interacting bodies onto the segments of the other body. The contact and separation which are governed by the Kuhn-Tucker conditions states t n ≤ 0, g n ≥ 0, t n .g n = 0 (19) where t n is the contact stress and g n represents the normal gap function [55] used to obtain the minimum distance between two points on the surfaces of the contacting body. By using the penalty method to enforce these constraints, the additional penalty energy from the contribution of the contact segment, i (on boundary C i ), C i , to the proposed fully coupled equations can be written as where ε c is the penalty coefficient. After linearization and the inclusion of the additional energy arising from the contact in M wÜ + Q T wU + C wwṖw + C wcṖc + H ww P w + H wc P c = F w (22) M cÜ + R T cU + C cwṖw + C ccṖw + H T wc P w + H cc P c = F c (23) where n s is the number of contact segments. The symbols F, M, Q, H, and, C in the above equations refer to the force vectors, and the mass, coupling, flow, and damping matrices. Furthermore, the symbols U,P w and P c represent the nodal displacement, pore water pressure and suction which are taken as the unknown variables in developing the presented finite element solution. In addition, the superimposed dot denotes the time derivative of a variable. All matrices and force vectors (using the same notations) are defined in Appendix A. The definitions of K NC i and F NC i can be found in the work presented in [52]. The generalised-α method is used for time integration of the global equations of motion [4]. In all the following examples the tolerance for the iterative scheme based on the Newton-Raphson method is set to 10 −3 . Furthermore, the penalty coefficient is set to 10 7 kN/m 3 when obtaining the numerical results of the analyses presented here.

Finite element results
The problem investigated here considers the response of unsaturated soil to dynamic indentation by a rigid impervious circular platen with a diameter of 300 mm by assuming small deformations. On the top of the rigid plate, we have applied a triangular impact load as shown in Fig. 5. The load has a peak of 6.5 kN and a duration of 0.028 s. As shown in Fig. 5, we idealised this problem as an axisymmetric body where a 2 m by 2 m domain of unsaturated soil is discretised by 900 nodes and 410 6-noded quadratic elements for displacement that are coupled with two 3-noded elements for pore water pressure and suction [4]. The side boundaries of the mesh are restrained against horizontal displacements whereas vertical displacements are not allowed at the bottom boundary. Drainage is only possible from the portion of the top boundary that is not in contact with the impervious plate. Unless stated otherwise, the initial void ratio that is assigned to the soil domain is 0.7, the initial suction is 10 kPa and the at-rest earth pressure coefficient is set to 0.4. Initially, we establish geostatic stress and then, we initiate dynamic of the plate resting on the soil. The impact process is simulated by using 1333 time steps of the size of 0.00003 s. The SWRC is shown in Fig. 6a (with p * c as the normalized suction defined in Eq. (68) in Appendix B). All Displacement, pore water pressure, suction, bonding variables, and the degree of saturation at the centre of the interface of the contact between the impervious plate and the unsaturated are also shown in Fig. 6. These results are obtained from simulations with no stress overshooting treatment. The results indicate that during the loading stage (i.e., the first 0.014 s of the analysis), pore water pressure and the degree of saturation increase whereas suction decreases at the monitoring point. The bonding variable has an inverse relationship with the degree of saturation per Eq. (1), therefore, it initially decreases. Nonetheless, upon the initiation of unloading, the aforementioned trend is reversed. The pore water pressure and suction results show strong oscillations that are initiated in the final stage of the impact. Overall, given the strong oscillations observed, this problem provides a good example for studying the performance of the discussed schemes for correcting for stress overshooting.
We solved the same contact-impact problem using the two repositioning schemes discussed earlier. The predicted displacements on the contact interface are shown in Fig. 7. It may be noted that in the absence of significant oscillations, there is very good agreement between the results obtained from the analyses with and without using the repositioning schemes. Also, when true unloading occurs (due to reversal of the applied load at t = 0.014 s), all the analyses are in perfect agreement as shown in Fig. 7. However, when oscillations become more pronounced, differences between the predicted displacements can be observed.
We conducted 45 analyses to explore the performance of the proposed models for a diverse range of initial suction and void ratio. The first set of analyses was performed by selecting the initial suction values of 0, 2, 5,  with no overshooting treatment, the model with the overshooting treatment proposed in [15], and the model with the overshooting treatment in this paper. We compared the CPU time and the number of iterations in the global generalised-α solver based on Newton-Raphson's method. In the second series of tests, we selected initial void ratios of 0.6, and 0.9. It should be noted that the discussed repositioning schemes generally have competing effects on the speed of the stress integration and the speed of the generalised-alpha solver. Both repositioning schemes add complexities as well as storage and computational costs to the integration scheme as outlined in Algorithm 1 to Algorithm 5. Nonetheless, if the proposed schemes can effectively mitigate the stress overshooting problem, the storage and computational costs can be compensated by the improvement in the stability of the solution, fewer iterations and a reduced substepping frequency, and overall enhancement of the speed of the analysis. Table 2 summarises the information on the relative CPU time and the number of iterations in the two series of analyses with and without overshooting treatment. The CPU of the system that is used for these analyses is Intel(R) Core(TM) i7-3630QM and the system is equipped with 16 GB of RAM. The operating system is Windows 10, 64-bit. Table 2 and Table 3 show that in the majority of the analyses, the use of the repositioning schemes resulted in a decrease in the CPU time and iterations. However, the results show a clear dependency on the analysis condition. The unsaturated cases generally exhibit higher CPU time and iterations closer to saturation and as the initial suction  375  376  1431  1431  1435   2  641  607  638  1616  1558  1603   5  662  651  667  1861  1721  1903   8  568  544  Unsuccessful  analysis   1792  1728  -10  543  511  561  1787  1720  1795   15  574  538  Unsuccessful  analysis   1812  1728  -20  510  520  515  1747  1765  1821   40  518  493  512  1848  1807  1826   60  541  501  534  1792  1784  1960   100  Unsuccessful  analysis   457  Unsuccessful  analysis   -1820  -200  522  500  549  1943  1928  1960   350  192  192  192  1344  1344  1344   500  179  181  178  1333  1333  1333 increases, the analyses were performed faster and with fewer iterations. Figure 8 highlights how some of the key observations in these coupled problems when the initial suction changes. It should be noted that all the results presented in this figure were obtained from the analysis with no overshooting treatment. We have defined r u = p w σ yy0 (i.e., the excess pore water pressure normalised by the initial vertical effective stress, σ yy0 ) and showed its evolutions at 2 cm below the centre of impact in Fig. 8.a in the analyses with the initial suction values of 2, 5, 60, and 500 kPa. The graph shows that r u can significantly reduce when the initial suction increases. Such a reduction in r u is due to smaller excess pore water pressure when the soil is in dryer states and higher initial effective stress (≈ − 2.0, − 4.6, − 29.0 and − 171.5 kPa corresponding to suction values of 2, 5, 60, and 500 kPa) when suction increases. The changes in volumetric strain, ε v with depth beneath the centre of impact when the peak of the load is applied (at 0.014 s) is shown in Fig. 8.b. The Table 3 Information on the efficiency of the analyses with and without the overshooting treatments in analyses with initial suction of 10 kPa and void ratios of 0. 6  versus time b volumetric strain, ε v with depth c ξ ξ0 with depth observed trend in the graph resembles the regime of volume changes below the centre of a rigid cylinder in contact with unsaturated soils analysed in [52] where far from saturation, the volume change beneath the centre of impact is purely compressive. Nonetheless, as the initial suction decreases, the formation of a dilative zone can be identified in analysis with initial suction of 2 kPa. Figure 8.c shows the changes in the bonding variable normalised by its initial value, ξ 0 (≈ 0.014, 0.074, 0.543, 0.8 corresponding to suction values of 2, 5, 60, and 500 kPa) with depth. It is seen that at the chosen time step, overall, the bonding variable decreases because the saturation degree increases and suction decreases. Nonetheless, in the analyses with initial suctions of 5 and 2 kPa, a zone can be identified where the bonding variable increases. In four instances (e = 0.7 with suction values of 8, 15, and 100 kPa in Table 2, and e = 0.6 with suction values of 10 in Table 3), the use of the repositioning scheme resulted in the successful completion of the analyses that failed when no repositioning scheme was used. In one instance (initial suction = 100 kPa in Table 2), the use of the proposed repositioning scheme resulted in the successful completion of the analysis, a result that could not be achieved when no overshooting treatment or the repositioning scheme proposed in [15] was used. Figure 9 compares the evolution of the vertical effective stress at 2 cm below the centre of impact obtained from the analyses with the initial suction value of 100 kPa in Table 2 with the two repositioning schemes. The point where the analysis with the repositioning scheme proposed in [15] failed is shown in the graph. These results indicate that both analyses are in agreement before the oscillatory phase of the analyses. However, a drastic jump in the predicted vertical effective stress is seen before the termination of the analysis with the repositioning scheme proposed in [15]. Nonetheless, relatively better performance is seen in the analysis with the proposed repositioning scheme and furthermore, the analysis could be completed.
In the final series of analyses, we investigated the performance of the proposed schemes when the penalty coefficient changes by conducting 12 additional analyses with the penalty coefficient being set to 1E5, 1E6, 5E6, and 8E6 k N /m 3 and keeping other aspects of the analysis the same. The penalty coefficient represents the stiffness of a continuous spring on the contact surface where higher penalty coefficients generally yield smaller penetrations. Compared with methods such as the "Augmented Lagrangian" [56], the penalty approach has the advantage of removing the need for additional degrees of freedom in the analysis [57] which improves the speed of fully coupled analyses. However, the results can exhibit sensitivity to the selected penalty parameter [52,57]. It should be noted that small penalty coefficients can generally lead to violation of the contact constraints in Eq. (19) as penetration cannot be effectively prevented. To fully enforce these constraints, the penalty coefficient should ideally approach infinity, however, due to numerical issues such as strong oscillations and ill-conditioning, the penalty coefficient cannot be too large. Therefore, often an optimum penalty coefficient can be selected that is large enough to effectively prevent penetration of the surfaces in contact but not too large that it can cause early termination of the analysis. As noted in [57] a general applicable mathematical approach for finding the optimum penalty coefficient cannot be developed, as finding such a penalty coefficient generally depends on the characteristics of the problem such as the material model, parameters, and boundary constraints. Therefore, trial and error may be needed in nonlinear contact problems to obtain the most suitable penalty coefficient. Figure 10.a shows the predicted displacement at the surface in analyses with various penalty coefficients chosen for study. Under a constant dynamic force, smaller penalty coefficients lead to higher peaks of displacement during the impact. Also, Fig. 11 indicates that penetration is more pronounced when smaller penalty coefficients are selected and for the case with the highest penalty coefficient (1E7 k N /m 3 ) no noticeable penetration can be observed. The pore water pressure changes at the centre of impact in the analyses with the lowest and highest penalty coefficients are compared in Fig. 10.b.
The results indicate that when the penalty coefficient is low, numerical oscillations are small, however, numerical oscillations can significantly increase when a penalty coefficient of 1E7 k N /m 3 is selected. Overall, effective prevention of penetration and increasing the accuracy of the predictions by selecting a high penalty coefficient will add numerical oscillations and challenge the stress integration scheme and the constitutive model (e.g., because of the potential for stress overshooting). Table 4 summarises the information on the efficiency of the analyses with and without the repositioning schemes when different penalty coefficients are used. As explained earlier, due to stronger oscillations, overall, CPU time and iterations increase when larger penalty coefficients are used. In the analysis with the lowest penalty coefficient, no changes in the iterations were observed. Nonetheless, CPU time and iterations in the analyses are reduced in both repositioning schemes when higher penalty coefficients are adopted.

Conclusions
We discussed two memory repositioning options for solving the problem of stress overshooting in a constitutive model for multiphase soils known as MUD. The first scheme was proposed in [15] and the second was proposed in this paper to resolve a highlighted numerical issue with the first repositioning scheme. Also, we presented adaptive explicit integration schemes for MUD with the two repositioning schemes. The advantages and limitations of both algorithms  were highlighted by several examples. Moreover, an oscillatory dynamic contact problem was solved to compare the performance of the two schemes. It was found that both schemes can generally improve the speed of the analyses and the number of iterations, hence helping the stability of the solutions. The results showed the dependency of the performance of these schemes on the initial analysis conditions. Nonetheless, an instance was highlighted where the proposed repositioning scheme could outperform the scheme proposed in [15]. Furthermore, we showed that the hardening law has a pivotal role in solving the overshooting problem. As stated in Sect. 2.1, upon saturation, the hardening law of MUD will degenerate to that in the "SANISAND" family of models. Therefore, the proposed solutions and the integration schemes can be readily applied to saturated soil models that employ a similar hardening law. It has also been demonstrated that to achieve more accurate results, the penalty coefficient should be large enough to effectively prevent penetration of the surfaces in contact. However, larger penalty coefficients led to stronger oscillations. In such conditions, both repositioning schemes could effectively reduce the iterations and CPU time.  The finite element matrices and vectors in Sect. 3 are defined by C cc = N T pc C * 2 N pc d (47) R c = ∫ B T C 6 mN pc d The load and the flow vectors appearing in the governing finite element equations are defined by: where b represents the vector of body force.
In the above equations, the quantities ρ g , ρ w , and ρ s are the mass densities of air, water and solid. Based on these definitions, the overall density of the mixture, ρ, is defined as: ρ = (1 − n)ρ s + nS w ρ w + nS g ρ g (54) with n as porosity and S w and S g as the degrees of saturation of the water and air phases.
The permeability for each phase, k β (β = w, g) (with the subscripts w and g representing the water and gas phases, respectively), is defined as where k int is the intrinsic permeability. The parameters η β and k r β are respectively the viscosity and the relative permeability (the ratio between the permeability at an unsaturated state and its counterpart at the saturated state) of the water and air phases. Also, by taking K w and K g as the bulk moduli of the water and air, we define.
C * 3 = S g n K g (58) The imposed Dirichlet boundary conditions of the primary variables are: u = u on u (60) p c = p c on p c (62) whereas the Neumann boundary conditions for the prescribed tractions and fluxes are: .n * =ẇ w on q w (64) k g −∇ p g + ρ g (b −ü) .n * =ẇ g on q g (65) whereẇ β (β = w, g) are the prescribed values of the outflow rate at the permeable boundaries q β (β = w, g). In addition, we define  with the terms n α (α = x, y, z) representing components of the outward unit normal vector to the boundaries in the x, y and z directions. Table 5 shows the parameters of the constitutive model and Table 6 demonstrates the parameters of the SWRC model.

Appendix B: Material parameters
The degree of saturation is assumed to be a function of suction and void ratio following the approach by taking is a material parameter, we define a modified suction, p * c as follows The equation of SWRC is defined in [37] and [38,58] S w = S r min + (S r max − S r min ) * ln exp(1) + p * c P a n x −m x where P a denotes the air-entry value. n x and m x are the model parameters that control the slope of the SWRC. Also, S rmax and S r min denote the maximum and minimum values of the residual degree of saturations.