A Constitutive Model for Locally Drained Shear Bands in Globally Undrained Dense Sand

When saturated granular materials which are dilative in nature are subjected to the undrained deformation, their strength increases due to the generation of negative excess pore pressure. This phenomenon is known as dilative hardening and can be witnessed in saturated dense sand or rocks during very fast loading. However, experimental evidence of undrained biaxial compression tests of dense sand shows a limit to this dilative hardening due to the formation of shear bands. There is no consensus in the literature about the mechanism which triggers these shear bands in the dense dilative sand under isochoric constraint. The possible theoretical reasoning is the local drainage inside the specimen under the globally undrained condition, which is challenging to be monitored experimentally. Hence, both incept of localisation and post-bifurcation of the saturated undrained dense sand demand further numerical investigation. Pathological mesh dependency hinders the ability of the finite element method to represent the localisation without advanced regularisation methods. This paper attempt to provide a macroscopic constitutive behaviour of the undrained deformation of the saturated dense sand in the presence of a locally drained shear band. Discontinuation of dilatant hardening due to partial drainage between the shear band and the adjacent material is integrated into the constitutive model without changing governing equilibrium equations. Initially, a classical bifurcation analysis is conducted to detect the inception and inclination of the shear band based on the underlying drained deformation. Then a post-bifurcation analysis is carried out assuming an embedded drained or partially drained shear band at gauss points which satisfy bifurcation criterion. The smeared shear band approach is utilised to homogenise the constitutive relationship. It is observed that the dilatant hardening in the saturated undrained dense sand is reduced considerably due to the formation of shear bands.


Introduction
Dilatant materials cause a reduction in pore pressure during undrained plastic shearing. This leads to an increase in effective stress which is popularly known as dilative hardening in saturated soils and rocks. If the bifurcation criterion is regarded as the singularity of the undrained acoustic tensor, this dilative hardening phenomenon can delay the instability. Hence the strengthening can continue until the tendency to dilate ceases at the critical state. The only exception is when the suction reduces enough to reach the cavitation pressure, then the soil is desaturated, lifting the strict isochoric constraint. This leads to the rationale that under large hydrostatic pressure in sub-sea conditions, dilative silty sand can possess very high undrained strength. However, [1], [2] showed using the perturbation analysis that homogeneous undrained deformation of dilative materials can be potentially unstable when the condition for localisation is met in terms of the underlying drained response. On the contrary, [3] states that in the absolute absence of inhomogeneities, no shear band can occur even though the drained localisation criterion is met. In practical circumstances, the specimens are not homo- * e-mail: hem42@cam.ac.uk geneous. A small imperfection can trigger local drainage. This phenomenon can possibly terminate the dilative hardening leading to a localised deformation in the globally undrained dense sand.

Literature Review
The strain localisation in the undrained dilative granular material is a subject of controversy. Observations of the experimental investigations carried out over the past few decades are not sufficient to reach a unique consensus. This is mainly due to the fact that coupled localisation is objective of the size of the specimen, aspect ratio, grain size, loading rate, material properties (friction, dilation, permeability) and initial conditions. [4] reported that the onset of shear band formation was delayed in the undrained dense sand until the cavitation, which was governed by the initial back pressure. [5] reported the evidence of locally drained shear bands despite global undrained conditions.
Due to the lack of experimental evidence, most studies on the locally drained shear bands in the globally undrained sand are limited to theoretical and numerical investigations. [6], [7] and [8] presented the concept of the locally drained shear band under globally undrained conditions. Only this type of instability is possible for dilative dense sand [7,8]. A local volume change is always accompanied by a water flow associated with a pore pressure gradient. This means there should be a jump in the volume change across the shear band when the pore pressure is dissipated.
The onset of the shear band is dictated by the rate of local volume change. [7] theoretically showed it is most likely that a shear band initiates with a minimum deviation from the isochoric constraint, after the peak friction ratio. Because in this region, the jump in local volume change required to trigger a dilative shear band is smaller. Hence, a small heterogeneity in the sample can trigger localisation. On the other hand, if the isochoric constraint is lifted due to cavitation, full drainage is possible, hence the localisation can occur very early [7]. Therefore, in most experimental conditions, the internal fluid flow is suppressed and cavitation can precede that. Moreover, there are possible time and scale effects which prevent local drainage in displacement controlled tests with impermeable boundaries [4,6].
Modelling the internal drainage within the finite element framework is not very reliable. Generally, fluid migration happens within the pore volume which has a much smaller scale than the element size. Several attempts have been made to allow the internal drainage within the elements. [3] developed an embedded shear band approach which allows strong discontinuities within finite elements. [9] augmented the constitutive relations governing the macroscopic undrained response in the presence of a locally drained shear band. This method captured the exchange of fluid inside the element. However, these strategies do not recognise the time dependence of pore fluid diffusion.

Objective
The objective of this study is to develop a phenomenological constitutive relationship for the globally undrained dense sand, which encompasses a fully or partially drained shear band. It should be able to capture the internal fluid movements and volume changes between the shear band and intact material taking both the rates of loading and the thickness of the band into consideration. The macroscopic strength is decided by the interplay between the excess pore pressure generation and diffusion at the grain scale. Thus, the influence of underlying micro-kinematics can be apprehended with sufficient accuracy without changing equilibrium equations.
The hypothesis of this development is illustrated in Figure 1. Figure 1(a) shows a globally undrained biaxial compression specimen with an extra-fine mesh. A locally drained shear band has been propagated in it due to the weak material points embedded at the bottom left corner. Figure 1(b) shows a single element with the same dimensions and boundary conditions. Impermeable boundaries are presumed in both tests to maintain the globally undrained and constant volume condition. The aim is to replicate the mechanical response of the extra-fine mesh by a phenomenological constitutive relationship of the single element. The assumption here is that the finite element analysis with an extremely fine mesh can reproduce the actual local drainage in dense sand. This is the most plausible alternative due to the lack of experimental evidence.

Numerical Implementation
The details of mathematical formulations derived for the embedded shear band approach is described in this section. A shear band with a finite thickness is assumed to occur in the single element/or material point after it has reached the bifurcation criterion ( Figure 1(b)). After bifurcation point, the material response no longer can be accurately described by a continuum constitutive model. The homogeneity might be valid separately in each region (inside and outside the band) but not as a whole. The mode of deformation inside and outside the shear band take different forms regardless of the boundary conditions of the overall element. Therefore, to accurately capture the postbifurcation response, two separate stress-strain relationships are required to describe the behaviour inside and outside the localisation zone [10]. The average macroscopic response is calculated based on the volume averaging procedure following works of [9], [10]. Finally, the averaged macroscopic stress is utilised to solve the equilibrium equations according to global boundary conditions. Hence, the evolution of underlying micro-deformation dictates the macro response. Main formulations of this embedded shear band approach are summarised below.
The numerical investigation is carried out using the finite element software ABAQUS with an user-defined subroutine. The onset of bifurcation and the inclination angle of the shear band should be calibrated for the extra-fine mesh for different loading rates. The bifurcation is decided upon the first time the drained acoustic tensor becomes zero in any material point in the extra-fine mesh. The shear band angle is assumed to be a constant although it slightly changes during the deformation.
All constitutive relations are developed assuming a plane strain condition. The homogeneous deformation prior to bifurcation is described with reference to the global coordinate system x-y. For clarity, the postbifurcation analysis is attached to a local coordinates system n-t which is aligned with the shear band. All calculations are conducted in local coordinates until they are transformed back to global coordinates. The transformation matrix is defined as The stress and strain measures after the transformation are denoted as A single material point is assumed to be consist of two sub-elements made out of the shear band and intact materials. Both share the same material properties but different modes of deformation. This approach integrates both the size and orientation of the shear band into constitutive relationships beyond the localisation. During post-bifurcation analysis, macro effective stress and strain increment vectors are symbolised by σ i j and d i j respectively where i, j represent n and t. Two constitutes involved: shear band and intact material are denoted by superscripts sb and out respectively. The area factor-µ is defined as the ratio of shear band area-A sb to the total area of the element or material point-A element . This can be approximated as the ratio of shear band thicknesst sb to the average element or material point lengtha.
The weak discontinuity approach or Hill-Mandel condition assumes a velocity jump or a change in the gradients of displacement across the boundaries of the shear band. Hence a linear scaling of the strain rate is valid.
A unique feature of the undrained deformation is negligible volumetric strains. Therefore, without comprehensive calculations, it is hypothesised that the shear strain is concentrated inside the shear band while the volumetric strain is shared by both shear band and intact material. Thus Equation 5 can be rewritten as At the onset of bifurcation, the material response inside and outside the shear band diverge from the same point. Therefore, the effective stresses, void ratios and excess pore pressures inside and outside the shear band are equalised to homogeneous effective stress-σ i j , void ratioe and pore pressure-U (negative) at the onset of bifurcation.
e sb = e out = e (8) During the post-bifurcation analysis, effective stresses, void ratios and pore pressures inside and outside the shear band are updated separately from time step n to n + 1 .
During the update, they should be corrected for the pore pressure generation and dissipation as given in sections 4.1 and 4.2 before being carried forward for the next step. To simplify the calculation, the generation and dissipation of excess pore pressure is decoupled as shown in Figure 2 (a) and (b), even though both occur simultaneously. The notations for stress, strain and excess pore pressures before the dissipation are symbolised by a superscript star. In shear band material, the pore pressure change is caused by both volumetric change and the shear induced dilation. Outside the shear band, the only reason for the pore pressure generation is the volumetric change.

Calculation of local drainage
Since the material behaviour inside and outside are locally homogeneous, their stress strain relationships take the general form.
The shear induced pore pressure rise dU sb dilation is calculated as the increase in the mean effective stress inside the shear band.
The generated pore pressure increment due to the volumetric change during undrained condition is where K f is the bulk modulus of fluid. Due to the difference in the accumulated pore pressures inside and outside the shear band, an excess pore pressure gradient is created as shown in Figure 2 (a). The main hypothesis of this analysis is that a certain percentage of this gradient should be dissipated within the same time step (t n+1 − t n = dt). After pore fluid dissipation the excess pore pressure inside the band reduces and the pressure outside the band increases as shown in Figure 2 (b). Therefore, Equations 10 can be corrected as The amount of dissipated pore pressure-dU sb dissipation depends on the difference in generated pore pressures inside and outside the shear band-∆U. Depending on the duration of loading and diffusive properties of the material, the pore pressure gradient can be accumulated from previous steps.
The magnitude of pore pressure remaining inside the band after the dissipation can be calculated from the diffusivity Equation 19. The amount of dissipation is governed by the diffusivity coefficientc v as well as the time stepdt of the analysis. c v depends on the permeabilityk and volume the compressibilitym v .
dU(z, dt) is the pore pressure profile at end of the time step. The initial pore pressure profile dU(z, 0) is assumed to be a step function with a magnitude of ∆U inside the shear band and zero outside. This is shown as the initial profile (red) in Figure 3. The main assumption here is that the dissipation happens after generation. It is observed in Figure 3, that depending on the loading duration, the shear band can be fully or partially drained. The dissipated amount of pore pressure during a single time step can be calculated as dU(z, dt) depends on the duration of external loading, hence this includes the time dependent behaviour in the context of time independent plasticity. Generally, it is assumed that the fully undrained response of the saturated soil is time independent. The hypothesis here is that although the global behaviour is independent of time, the local drainage depends on the time, hence should be taken into the consideration. dU out recived depends on drainage distance as shown in Figure 3. For very small area factor, dU out recived can be assumed as zero.

Correction of stress for the local drainage
The pore pressure gradient created due to enhanced shear strain inside the band is shown in Figure 2(a) in blue colour while the homogeneous pore pressure is shown in black dashed line. During the loading of saturated soil, the dissipation of excess pore pressure also occurs simultaneously with the generation. If the loading is slow enough, the excess pore pressure inside can be dissipated making the shear band fully or partially drained which is shown in Figure 2(b). The remaining excess pore pressure after dissipation contributes to the final effective stresses in and out of the shear band. Therefore, the ultimate strength at the end of a time step results from the competition between the rate of generation and dissipation, which is governed by the permeability, compressibility of soil as well as the loading rate. Hence, updated stress inside the shear band from Equation 13 should be corrected as σ sb The outside effective stress is assumed not to be changed by the dissipation.

Calculation of volumetric strain inside the shear band
Due to the undrained global boundary conditions, the volumetric strain at a material point is very small. However, it is delineated in the above section that a portion of excess pore pressure generated due to undrained mean pressure rise inside the band is reduced due to the dissipation. This changes the boundary conditions of shear band material from undrained to partially drained. In other words, there is a net volumetric strain increment inside the shear band due to the dissipated pore pressure. The corrected strain increments inside the shear bandd sb i j can be recalculated from the corrected stress increments.

Calculation of macroscopic stress
The macroscopic stress of the material point is established by averaging the corrected shear band stress (Equation 21) and outside stress (Equation 22). The foundation for this averaging technique stems from combining the average principal of virtual work and the strain decomposition in Equation 5 [10].
Finally, the homogenised effective stress at the end of the time step is transformed back to global coordinates xy. It is utilised by ABAQUS to solve equilibrium equations at nodes. Corrected stresses, void ratios and excess pore pressures of respective material are carried forward for the next time step. The proposed constitutive method is termed as the diffusion SB model hereafter.

Calibration of diffusion SB model with undrained extra-fine mesh
Single element biaxial compression test is conducted to numerically validate the proposed constitutive model. It is compared with extra-fine mesh (element size 0.006 25 m) with same boundary conditions. Some weak material points are embedded in the bottom right corner of the extra-fine mesh to trigger the localisation. Plane strain reduced integration elements with pore pressure (CPE8RP) are utilised for all simulations. The initial and boundary conditions are specified in the first step. In both specimens, the bottom right node is pinned and other bottom nodes are roller supported. The top and side nodes are not constrained. All vertical and horizontal boundaries are assumed to be impermeable. The initial void ratio of sand is 0.55 and the permeability is 0.001 m/s. The initial pore pressure is prescribed as 10 kPa throughout the specimens. A confining pressure of 100 kPa is applied to both top and side edges during the second step. A transient consolidation analysis is conducted with ABAQUS finite element software. At the final step, a displacement of 0.1 m is applied to top nodes. The different loading durations are mentioned in Table 1. Both specimens are assumed to remain saturated throughout the deformation. Hence the phenomenon of cavitation is ruled out.
Formulations in section 4 are valid for an arbitrary critical state material model. In this study, the Nor-Sand constitutive model is adopted. A detailed implementation of this model can be found in [11]. Since it is a state parameter based model, the dilation is a function of the current void ratio and effective mean pressure. This facilitates the implementation of the proposed approach.
The extra-fine mesh is simulated with the original Nor-Sand model whereas the single element is modelled with the proposed diffusion SB model. The thickness of the shear band is assumed to be 0.0125 m, similar to the band thickness of extra-fine mesh. The root mean square average of a gauss point in the single element is 0.175 m. Hence the area factor calculated from Equation 4 is 0.071. The shear band angle is 54 • which is calculated from the extra-fine mesh. Both shear band inclination and thickness are assumed to be constants and measured after a steady state shear band is formed in the extra-fine mesh.
The onset of bifurcation and diffusive coefficient in the diffusion SB model are calibrated to match with the global response of the extra-fine mesh. The calibrated input parameters of the diffusion SB model are given in Table 1. In fact, the shear localisation in the single element starts when the acoustic tensor criterion is first time satisfied by the extra-fine mesh. Figure 4 displays the total reaction forces of extra-fine mesh simulated by original Nor-Sand model along with the single element simulations of the diffusion SB model. The response of the single element modelled by the original Nor-Sand model is also included for the comparison.   It is observed that the reaction forces of the diffusion SB model reach constant plateaus after the shear bands are formed for the loading duration 1 s,0.1 s and 0.01 s. According to the assumption 6(a), the outside shear strain is zero, thus the mean pressure increases only inside the shear band. When the loading time step is large enough to fully dissipate the dilation induced pore pressure inside the shear band, effective mean pressure becomes a constant. Hence there is no longer an increase in shear strength. However, when the loading duration is 0.001 s, initially pore pressure gradient accumulates inside the shear band which leads to a slight increase of stiffness as shown in Figure 4 (orange dotted line). Nevertheless, large pressure gradient, in turn, accelerates the dissipation and reduces the shear strength soon. Another reason for the observed plateaus is the stress averaging in Equation 26. When the area factor is small (shear band thickness is much smaller than average element/or material point length), the macroscopic strength is governed by the outside material.
In reality, the strength in material outside the shear band also increases due to two reasons. First, it also undergoes a slight dilation and produces negative pore pressure. Second, it receives dissipated negative pore pressure from the shear band. The diffusion SB model does not take this into consideration since the amount of received pressure depends on the proximity to the shear band. The proposed model is created for the ultimate motive of modelling large scale boundary value problems, for which the shear band thickness is several orders of magnitude smaller than the element size. Therefore the received pore pressure is ignored here. As a result, the macroscopic mean pressure increase is hindered after the shear band is formed.
Several factors should be taken into the account when the global response of the extremely fine mesh is compared with the material response of the diffusion SB model. Localisation in a continuum mesh is progressive in nature. Different material points reach bifurcation criteria at different stages of deformation. Generally, weak elements initiate the shear band. The change of stiffness in the global response of extra-fine mesh is likely when the acoustic tensor criterion is first time satisfied by any material point. This generally occurs after the maximum stress ratio. This heterogeneity cannot be captured by a homogeneous constitutive model. As a solution, in Figure  4 the onset of the shear band is calibrated to match with the extra-fine mesh. In a similar attempt to build a material model with a locally drained shear band, [9] utilised the experimentally observed values for the shear band initiation.
Moreover, even when the shear band is fully drained (when the loading duration is greater) the reaction forces in extra-fine mesh continue to increase. This is because different material points reach the critical state at different stages and the strength is continuously building up. The shear band develops progressively widening its thickness. Further, outside elements unload at the start of localisation and reload again. These features cannot be accounted by the constitutive relations in the diffusion SB model. Hence a small amount of shear should be included in outside material to match with the results of extra-fine mesh.
Even though the post-localised progressive failure is not apprehended very accurately, the results of diffusion SB model are superior to the undrained behaviour predicted by the original Nor-Sand model (displayed by a red dash line in Figure 4). The original model cannot perceive the underlying micro-kinematics. Hence it shows an uniform undrained response independent of the rate of loading.

Conclusion
The negative pore pressure generated during the undrained deformation of dilative sand results in an enhanced strength which is called dilative hardening. This phenomenon can occur in the saturated dense sand with a higher percentage of fine particles during fast loading rates. However, the shear localisation associated with the local drainage interrupts the dilative hardening and terminates the strength increase. In finite element method, this local drainage is mesh dependent and correctly captured only if the mesh size is smaller than the shear band thickness. A macroscopic constitutive model is developed in this study to capture the local drainage which occurs at micro-scale. Both rate and scale effects which govern the pore fluid diffusion in the shear band are taken into consideration. The proposed method is calibrated with the finite element response of a very fine mesh. For smaller loading rates, the strength increase is terminated by the proposed model indicating a fully drained shear band. For faster loading rates, the shear strength is enhanced initially due to the accumulated pore pressure gradient but it is decreased eventually as the excess pore pressure is dissipated. The ultimate strengths predicted by this method are much lower than the homogeneous undrained strength simulated by the original constitutive model.