A two‐scale constitutive framework for modelling localised deformation in saturated dilative hardening material

Undrained deformation of dilative sand generates negative excess pore pressure. It enhances the strength, which is called dilative hardening. This increased suction is not permanent. The heterogeneity at the grain scale triggers localisations causing local volume changes. The negative hydraulic gradient drives fluid into dilating shear zones. It loosens the soil and diminishes the shear strength. It is essential to understand the mechanism behind this internal drainage and to capture it numerically. The purpose of this paper is to develop a macroscopic constitutive relationship for the undrained deformation of saturated dense sand in the presence of a locally fully or partially drained shear band. Separate constitutive relations are generated for the band and intact material. Both time and scale dependence during pore fluid diffusion in saturated sand are captured, eliminating the mesh dependency for finite element implementations. The model is applied to the Gauss points that satisfy the bifurcation criterion. The proposed method is calibrated to recreate the undrained macroscopic response bestowed by an extra‐small mesh. The microscopic behaviours inside and outside shear band predicted by this model are qualitatively in good agreement with individual material point behaviours inside and outside the shear band in the extra‐small mesh. Depending on the loading rate and the shear band thickness, the response inside the band can be fully or partially drained, which governs the ultimate global strength. The calibrated model is exploited to simulate an upscaled biaxial compression test with semipermeable boundaries.

dilation ceases at the critical state unless it is interrupted by shear localisation. If suction increases sufficiently to reach the cavitation pressure, the sand may desaturate releasing the isochoric constraint. There are several experimental and numerical studies that state that the localisation in globally undrained dense sand is delayed until the cavitation takes place. [1][2][3] However, in deep-sea condition, the hydrostatic pressure can suppress the cavitation. 4,5 In these circumstances, localisation induced by local volume changes precedes.
Vardoulakis, 6 Guo 7 and Guo and Stolle 8 theoretically presented the concept of the locally drained shear bands under globally undrained conditions. It is signalled by the nullity of the drained acoustic tensor when the undrained acoustic tensor is still positive. The onset of these shear bands depends on permeability, compressibility, the grain size of soil and the rate of loading. Han and Vardoulakis 9 provided experimental evidence on the local drainage in saturated dense sand during globally undrained biaxial compression. Tests were conducted under sufficient back pressure to avoid cavitation. The failure occurred after the Coulomb state of the background drained behaviour. X-ray radiographs showed conjugate periodic shear band arrays corresponding to alternative dilatant and contracting behaviour, which are linked together by the internal fluid flow. However, their observation was not substantiated by subsequent experiments of Roger et al. 1 and Mokni and Desrues. 2 They ascribed it to possible time and scale effects that prohibit internal fluid flow when the fluid flux across the boundaries of the specimen is zero. Nevertheless, previous numerical studies using finite element analysis have demonstrated the possibility of local drainage during globally undrained deformation. [10][11][12] However, they did not adequately address the mesh sensitivity during internal pore fluid diffusion.
In a conventional multiphase finite element method, pore pressure values are calculated at nodes similar to displacement. Hence, the generation of pore pressure gradient and its dissipation depend on the relative positions of nodes. In practical geotechnical applications, viable element sizes used in simulations are many orders of magnitude larger than the actual shear band thickness. These elements overlook the internal fluid flow, which happens at a finer length scale. In a drained analysis, the onset of localisation is mesh independent and can be detected solely based on prebifurcation constitutive relationship. Hence, regularisations or adaptive meshing techniques can be used for the post-localised deformation. In undrained dilative hardening materials, the initiation of shear band triggered by local drainage is also mesh sensitive. Because the mesh refinement is extremely expensive for large-scale problems, a strategy to bridge the micro and macro response is necessary.
To simulate local drainage inside undrained elements, Larsson et al. 13,14 and Larsson and Larsson 14 implemented strong discontinuity method. The regularised internal discontinuities of pore pressure and displacement are implemented in mixed variational formulations from 'enhanced strain' method. Two additional degrees of freedom, enhanced strain and enhanced pore pressure, are included to regularise discontinuities. The shear band is activated in elements that satisfy the lowest of either drained or undrained localisation criteria. Equilibrium equations are solved both at the element and the global levels. Numerical results indicated that embedded localisation minimised the mesh sensitivity during pore fluid diffusion. However, this approach needs tedious changes of equilibrium equations and involves extra degrees of freedom. Instead, Pietruszczak 15,16 opted to augment the constitutive relations governing the macroscopic undrained response in the presence of a fully drained or undrained shear band. This method captured the exchange of fluid inside the element (between the shear band and intact material) under globally undrained constraint. The macroscopic response is formulated by homogenisation based on volume averaging of both constitutes. As the thickness of the shear band is included in constitutive relations, this method is inherently mesh independent. To determine the point of transition from homogenised to localised deformation under undrained constraint, Pietruszczak 15 decided to use a physically based criterion rather than the conventional bifurcation analysis. So he selected the onset of the shear band based on experimental data. His results delineated that internal drainage at micro scale reduced the undrained strength of dilative hardening soil. This method is quite simple and promising. The drawback of this approach is that it cannot take account of the rate effect involved with pore fluid diffusion. Hence, it can only model fully drained or undrained shear bands. It will be interesting to explore the diffusion rate effect of partially drained shear bands inside undrained elements.

| OBJECTIVES
The primary objective of this study is to develop a phenomenological constitutive relationship for the undrained fully saturated dense sand, which is embedded with 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 size of the band into consideration. Therefore, both time and scale effects can be captured at the constitutive level. The macroscopic strength of undrained sand is decided by the excess pore pressure generation and diffusion at the shear band 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 1A displays a globally undrained biaxial compression specimen with an extra-small mesh. A locally drained shear band is propagated in it due to the weak material points embedded at the bottom right corner. Figure 1B shows a single element (or material point) 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-small mesh by a constitutive relationship of a single material point and ultimately use it for large-scale boundary value problems.
The objective here is to reproduce the local fluid flow of dense sand simulated with an extra-fine mesh. This is the most plausible alternative due to the lack of experimental evidence on local drainage in globally undrained dense sand. Another important assumption here is that isochoric constraint is not broken at the localisation. In other words, desaturation or phase change does not occur. This assumption is valid for the deep-sea condition where the hydrostatic pressure is large enough to delay the cavitation.

| NUMERICAL PROCEDURE FOR THE EMBEDDED LOCALISATION APPROACH
The details of mathematical formulations derived for the embedded shear band approach are described in this section. The numerical investigation is carried out using the finite element software ABAQUS with a user-defined subroutine. It defines the constitutive relationships and returns the updated stress to the finite element code.
The hierarchical formulations of the overall process are encapsulated in the flow diagram in Figure 2. Conventional constitutive formulations are carried out until the bifurcation criterion is met using macro variables (left of Figure 2). A shear band with a finite thickness t sb and a constant inclination angle θ sb is assumed to initiate in material points that have reached the bifurcation deviatoric strain ϵ sb q . A detailed discussion about the onset of bifurcation in the saturated undrained dense sand given in Section 3.1.
After the bifurcation point, the material response cannot 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 modes 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 post-bifurcation response, two separate stress-strain relationships are required to describe the behaviour inside and outside the localisation zone. 17,18 Hence, a single material point is assumed to consist of two subelements 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 the orientation of the shear band into constitutive relationships beyond the localisation. The right side of Figure 2 surmises microscopic constitutive relations inside and outside the band contemplating the pore fluid migration between two constitutes. The average macroscopic response is calculated based on volume averaging procedure following other works. [16][17][18][19] Finally, the averaged macroscopic stress is utilised to solve the equilibrium equations according to global boundary conditions. The evolution of underlying micro deformation dictates the macro response. Main formulations of this embedded shear band approach are summarised below.

| Determination of onset, inclination and thickness of the embedded shear band
The proposed embedded shear band method has the same school of thought as smeared shear band approach by Pietruszczak and Mroz 20 and Pietruszczak and Stolle 21,22 and two-scale constitutive framework by Nguyen 23 and F I G U R E 2 Flow chart for the embedded shear band method [Colour figure can be viewed at wileyonlinelibrary.com] Nguyen et al. 17,18,24 They have detected the onset and inclination of a shear band from the instability of the macroscopic constitutive equation. At the onset of bifurcation, the determinant of acoustic tensor ceases to be positive. The acoustic tensor is calculated from the current stiffness matrix, which depends on the current stress and stress path. Hence, this bifurcation criterion is applicable for any arbitrary loading path and boundary conditions. However, all these literature have only considered single-phase strain-softening plastic models or damage models. Hence, the acoustic tensor becomes singular at peak (for associative plasticity) or slightly before the peak (for nonassociative plasticity) of the stress-strain relationship. In these cases, localisation is triggered as a result of material softening.
However, the mathematical criterion for the onset of bifurcation in undrained dense sand (dilative hardening material) is not conclusive over the past few decades. According to Rice 25 and Rudnicki and Rice, 26 undrained deformation of dilative material becomes unstable only in the softening regime of underlying drained behaviour (signalled by the singularity of the drained acoustic tensor). By analysing acceleration waves under dynamic conditions, Loret and Harireche, 27 Loret and Prevost 28 and Vardoulakis 29 concluded the same; the singularity of the drained acoustic tensor marks the onset of bifurcation in dilative hardening materials.
On the contrary, Larsson and Larsson 30,31 stated that in the absolute absence of inhomogeneities, no shear band could occur even though the drained localisation criterion is met. Guo 7 and Guo and Stolle 8 also mathematically proved that the strict isochoric constraint prevents shear bands in dense dilative sand. Andrade and Borja 11 also showed that 'ideally homogenous' undrained dense sand does not show any localisation. The first author further verified this notion. For undrained deformation of dense sand, both drained and undrained acoustic tensors never become singular for two loading paths (simple shear and biaxial compression). 32 This is because second-order work is strictly positive due to the generation of negative excess pore pressure. All these evidence support Pietruszczak's statement that 'the classical bifurcation analysis is not very reliable in the context of a quasi-static undrained deformation'. 15 Due to the reasons mentioned above, the authors concluded that the conventional bifurcation analysis fails to detect the onset of shear bands in globally undrained dilative hardening materials, because it is not a material behaviour. In reality, heterogeneity of soil triggers shear zones associated with local drainage. Therefore, the authors decided to calibrate the onset and inclination of the locally drained shear band from the average behaviour of a numerical biaxial compression test. Local drainage is triggered by prescribing some weak material points. The onset of shear band coincides with the singularity of drained acoustic tensor in weak Gauss points. Hence, localisation is more likely to depend on the loading path of the weak material point.
The physical shear band thickness of sand is experimentally found to be a function of mean grain size. 33 Therefore, in this study, it is assumed as a material constant and hence calibrated with the numerical biaxial compression test. In the following sections, ϵ sb q , θ sb and t sb are input parameters that should be calibrated for the selected extra-fine mesh for different loading rates.

| Transformation of stress and strain to the local coordinate system
All constitutive relations in this study 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 post-bifurcation analysis is attached to a local coordinates system n − t, which is aligned with the shear band. The directions parallel and perpendicular to the shear band are denoted as t and n, respectively ( Figure 1B). When the criterion for the bifurcation is met, macro stress and strains are transformed from global coordinates x − y to local coordinates n − t. All postbifurcation formulations are conducted in local coordinates n − t until they are transformed back to global coordinates x − y at the end of the time step. The transformation matrix is defined as T ½ = cos 2 θ sb sin 2 θ sb 0 2 sin θ sb cos θ sb sin 2 θ sb cos 2 θ sb 0 −2 sin θ sb cos θ sb 0 0 1 0 −sin θ sb cos θ sb sin θ sb cos θ sb 0 cos 2 θ sb −sin 2 θ sb : During post-bifurcation analysis, macro effective stress and strain increment vectors are symbolised by σ 0 ij and dϵ ij , respectively, where i, j represent n and t. Two constitutes involved, shear band and intact material, are denoted by superscripts sb and out, respectively.
3.3 | Averaging of stress and strains according to area ratio 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 . It can be approximated as the ratio of shear band thickness t sb to the average element or material point length a. For reduced integration elements, the length of a material point is half the element size.
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 as shown in Figure 3. Hence, a linear scaling of the strain rate is valid according to Equation (5).

| Decomposition of strain rates into the shear band and outside material
Previous researchers used a two-scale approach to predict the post-localisation response of dry or drained geomaterials. At a very early stage, Pietruszczak and Mroz 20 assumed purely elastic strain outside and plastic strain inside the shear band. Later on, Pietruszczak and Niu, 19 Pietruszczak and Stolle 22 and Nguyen 23 conducted more accurate and rigorous calculations on the strain rate inside the shear band using the traction continuity. Thus, the internal equilibrium at shear band boundary is maintained. In this study, a two-scale method is used to explore post-localisation of undrained geomaterials. A unique feature of undrained deformation is negligible volumetric strains. Therefore, without comprehensive calculations, simple assumptions are made on strain rates inside and outside the shear band considering two scenarios. They are derived from the undrained simple shear deformation shown in Figure 4.
• Assumption Case 1. The shear strain is concentrated inside the shear band, whereas the volumetric strain is shared by both shear band and intact material.
• Assumption Case 2. The volumetric strain is shared by both shear band and intact material equally. A small shear strain is allowed outside the shear band. Satisfying linear averaging of strain, Equation (6) can be modified as

| Divergence from the bifurcation point
At the onset of bifurcation, the material responses inside and outside the shear band diverge from the same point in the stress-strain and void ratio curves as shown in Figure 5. Therefore, effective stresses, void ratios and excess pore pressures inside and outside the shear band are equalised to homogeneous effective stress σ 0 ij , void ratio e and pore pressure U at the onset of bifurcation.
To simplify the calculation, the generation and dissipation of excess pore pressure are decoupled as shown in Figure 6A,B even though both occur simultaneously. A superscript star denotes stress and pore pressure values after generation but before the dissipation of excess pore pressure.

| Calculation of stresses and pore pressures before dissipation
During the post-bifurcation analysis, effective stresses and void ratios inside and outside the shear band are updated from time step n to n + 1. This is done in two steps: before and after the dissipation. The update of variables after the generation of pore pressure is given in this section.
The stiffness matrices are calculated from stress values remembered from the end of previous time step.
The stiffness matrix of the shear band D sb ijkl is calculated from the shear band stress σ 0sb ratio e sb ½ n , whereas the stiffness matrix of intact material D out ijkl is calculated from outside stress σ 0out ij h i n and outside void ratio e out ½ n .
During post-localisation analysis, excess pore pressures become nonhomogeneous as shown in Figure 6A. In shear band material, the pore pressure change is caused by both volumetric change and the shear-induced dilation. This exceeds the homogeneous pressure calculated at nodes from global equilibrium equations and continuity equations of fluid (shown in black dash line). Outside the shear band, the only reason for the pore pressure generation is the volumetric change for Case 1. However, for Case 2, a small amount of dilation-induced pore pressure occurs in outside material as well. Hence, Equation (13b) should be modified accordingly.
Shear-induced pore pressure rise is calculated as the increase in mean effective stresses in respective materials.
The pore pressure generated due to volumetric change during undrained condition can be assumed as where K f is the bulk modulus of fluid.

| Calculation of dissipated pore pressure
Due to the difference in dilation-induced pore pressures inside and outside the shear band, an excess pore pressure gradient ΔU is created as shown in Figure 6A, which is calculated as Depending on the duration of loading and diffusive properties of the material, the pore pressure gradient can be accumulated from previous steps. 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). The dissipated amount dU sb dissipation is governed by the pressure gradient ΔU, time step dt of the analysis as well as the diffusivity coefficient c v , which is governed by permeability k and volume compressibility m v . The magnitude of pore pressure remaining after dissipation can be calculated simply from the 1D diffusivity (Equation 18). dU(z, dt) is the pore pressure profile at the 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. Figure 7 displays dUðz,dtÞ ΔU for different loading durations. The initial profile is shown in red.
= 0 otherwise ð17bÞ The dissipated amount of pore pressure in the middle of shear band (z = 0) is dU(z, dt) depends on the duration or rate 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 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 and hence should be taken into account.

| Correction of stresses and pore pressures after dissipation
The enhanced shear strain inside the shear band produces dilation-induced pore pressure inside. However, during the loading of saturated soil, dissipation of excess pore pressure co-occurs with generation. If the loading is slow enough, excess pore pressure inside can be dissipated, making the shear band fully or partially drained. After dissipation, the excess pore pressure inside the band reduces, and the pressure outside the band increases as shown in Figure 6B. Therefore, the pore pressures are corrected as The remaining excess pore pressure after dissipation contributes to the final effective stress in and out of the shear band. Hence, stress inside the shear band should be corrected as dU out recived depends on drainage distance as shown in Figure 7. If the shear band thickness is at least 8 times smaller compared with the length of the material point, dU recived is negligible. Hence, the outside effective stress is assumed not to be changed by the dissipation. 3.9 | Calculation of local volume changes inside the shear band Due to undrained global boundary conditions, the volumetric strain at a material point is negligible. However, it is delineated in the above section that a portion of excess pore pressure generated due to undrained mean pressure increase inside the shear 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 band dϵ sb − local ij can be recalculated from the corrected stress increments. Since intact material acts as almost undrained, dϵ out ii is negligible.
The void ratio inside and outside should be updated according to the respective volumetric strains.

| Calculation of macroscopic stress
According to the principle of virtual work equation, the work produced by macroscopic stress and strain rate should equilibrate the volume average of work due to stress and strain rates inside and outside the shear band.
Equation (26) can be satisfied if the volume averaging of stress is assumed along with the traction continuity across the shear band. 18 Hence, the macroscopic stress of the material point is established by averaging corrected shear band stress and outside stress.
Finally, the homogenised effective stress at the end of the time step is transformed back to global coordinates. It is utilised by ABAQUS to solve equilibrium and flow continuity equations at nodes. Corrected stresses, void ratios and excess pore pressures of respective material are carried forward for the next time step. The aforementioned method is termed as diffusion SB model hereafter.

| CALIBRATION OF DIFFUSION SB MODEL WITH UNDRAINED EXTRA-SMALL MESH
Two biaxial compression tests are conducted to validate the proposed constitutive model numerically. Figure 8 specifies the dimensions and boundary conditions of each specimen. Figure 8A contains 3200 elements (mesh size is 0.00625 m), whereas Figure 8B consists of a single element. 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. A transient consolidation analysis is conducted with ABAQUS finite element software.
The initial and boundary conditions are specified in the first step. In both specimens, the bottom left 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. At the final step, a displacement of 0.1 m is applied to top nodes under varying loading durations. Both specimens are assumed to remain saturated throughout the deformation. Hence, the phenomenon of cavitation is ruled out.
The proposed diffusion SB approach outlined in Section 3 is applied to Nor-Sand (NS) constitutive model, which is built upon critical state soil mechanics. 34 Main formulations of the model are given in Appendix A1, and material parameters are presented in Table 1. The extra-fine mesh in Figure 8A is simulated with the original NS model, whereas the single element in Figure 8B is modelled with the proposed diffusion SB model within the NS model framework.
The thickness, inclination angle and onset of the shear band are evaluated from the extra-fine mesh for different loading rates. 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 area factor is calculated as the ratio of this shear band thickness to the average length of a Gauss point in Figure 8B. The shear localisation in the single element starts when the acoustic tensor at any material point of extra-small mesh first becomes zero. The validity of both assumptions mentioned in Section 3.4 is examined. The shear factor is calibrated to match with the global response of the extra-small mesh. All calibrated parameters for each strain rate are mentioned in Table 2.

| Global response
It is observed that the single element reaction forces of the diffusion SB model (Case 1) reach constant plateaus after the shear band is formed for the loading strain rates 0.2, 2 and 20 s −1 . In Case 1, 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 strain rate is 200 s −1 , initially, the pore pressure gradient accumulates inside the shear band, which leads to a slight increase of stiffness as shown in Figure 9 (orange dotted line). Nevertheless, a 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 (27). When the area factor is small (shear band thickness is much smaller compared with the average material point length), the macroscopic strength is governed by the outside material. In Case 2, material outside the shear band is allowed to shear a little. Hence, a certain amount of dilation-induced negative pore pressure occurs outside the band as well. Therefore, a close resemblance to the extra-small mesh can be witnessed from results of diffusion SB model (Case 2). In reality, 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 consider this because the amount of received pressure depends on the proximity to the shear band. Here, the band thickness is minimal compared with material point size. Hence, strength enhancement due to received pressure is ignored.
The red line in Figure 9 displays the single element undrained response of the original NS model, and it is independent of loading rates. It is noticed that the results of the diffusion SB model are closer to those of extra-small mesh than the uniform undrained deformation predicted by the original NS model. Several factors should be taken into the account when the global reaction forces of the extra-small mesh are compared with constitutive relations of the diffusion SB model. The localisation in a continuum mesh is progressive in nature. Different material points reach bifurcation criteria at various stages of deformation. Generally, weak elements initiate the shear band. Therefore, the change of stiffness in the global response of extra-small mesh appears when acoustic tensors of weak elements start to become null. Thus, the drained localisation of the extra-small mesh is used as an input value for the proposed diffusion SB model. It happens after the Coulomb state of the homogeneous deformation. This is experimentally justifiable since Vardoulakis 6 mentioned that the onset of drained shear band occurred after the Coulomb state. Guo 7 also theoretically proved that local volume change in an undrained deformation is most likely after the maximum stress ratio.
Moreover, when the shear band is fully drained (when the loading duration is greater), the reaction forces in the extra-small mesh continue to increase. This observation is not dramatic for the partially drained cases. This is because material points reach the critical state at different stages. The shear band develops progressively, widening its thickness and sucking the pore fluid. Further, outside elements unload at the start of localisation and reload again. These features cannot be apprehended by constitutive relations of the diffusion SB model. Hence, a small amount of shear should be included in outside material to match with the results of the extra-small mesh.

| Material response inside and outside the shear band
As the global response is close enough, it is worthwhile to investigate the constitutive behaviour inside and outside the shear band predicted by the diffusion SB model. Figure 11 illustrates the stress, strain as well as volumetric relationships inside and outside the shear band for loading strain rate 200 s −1 predicted by diffusion SB model (Case 2).  Figure 10A that after bifurcation, the deviatoric strain is concentrated inside the band. Local volume changes also occur inside the band despite global isochoric constraint in Figure 10B. Both of these characteristics are captured by the diffusion SB model in Figure 11A,B. Figure 10C shows that the deviatoric stress inside the band reaches a peak after bifurcation and softens until the critical state, whereas deviatoric stress outside the band keeps on increasing albeit with a reduced stiffness. This is due to the previously mentioned reason that material outside the band rehardens by receiving the negative pore pressure dissipated from the band. Physically, dilative shear band is sucking the pore fluid. This phenomenon is replicated in Figure 11C, albeit with a comparatively higher peak deviatoric stress inside the band. Similar to Figure 10D, stress ratio inside the band reaches the critical stress ratio in Figure 11D.
In Figure 10E,G, stress path and void ratio path inside the band reverse and reach the critical state line (CSL) after localisation. Hence, the state parameter inside the band in Figure 10F approaches to zero. It indicates that, after bifurcation, the material inside the band approaches the critical state, whereas deformation outside the band is hindered. These results are consistent with experimental observations inside the localised zones from stereophotogrammetry, digital imaging correlation, X-ray tomography 35 as well as discrete element simulations. 36 In Figure 11E-G, the diffusion SB model also successfully recognises these distinct mechanical responses inside the band.
Most importantly, both Figures 10H and 11H display the difference between the accumulated pore pressures in and out of the band. Because the loading strain rate is too fast to dissipate the dilation-induced excess pore pressure inside the band, it accumulates initially. This is because the rate of generation is greater than the rate of dissipation. However, enhanced pore pressure gradient, in turn, accelerates the pore pressure dissipation. Therefore, both accumulated pore pressure and mean effective pressure start to decrease. This net pore pressure gradient demonstrates the partially drained behaviour of the shear band. It should be mentioned here that in the extra-small mesh, pore pressure values are taken from the node closest to the material point considered.
For comparison, the local volume and excess pore pressure profiles of loading strain rate 0.2 s −1 are displayed in Figure 12. Figure 12A,C depicts the response of arbitrary material points inside and outside the shear band in the extrasmall mesh simulated by the original NS model, whereas Figure 12B,D displays subelemental response of the diffusion SB model (Case 2). The main differences between loading strain rates 0.2 and 200 s −1 are (i) the early initiation of the shear band, (ii) greater local volume changes inside the band and (iii) zero pore pressure gradient across the band. For loading strain rate 0.2 s −1 , the generated pore pressure dissipates within the same time step, indicating a fully drained band. Still, the negative excess pore pressure uniformly accumulates throughout the specimen as shown in Figure 12C. The diffusion SB model sufficiently replicates all these features. However, in the extra-fine mesh, the material outside the shear band shows slight compaction that cannot be captured by the diffusion SB model. Apparently, for this loading rate, the global volume conservation mentioned by Vardoulakis 6 is satisfied. The volume increase inside the band is compensated by the contraction outside, such that isochoric constraint is maintained globally.
The main achievement of this diffusion SB approach is that constitutive relations alone capture both macroscopic and microscopic (inside and outside the shear band) behaviours of the extra-small mesh. Furthermore, these constitutive responses are sensitive to both time and scale effects observed in the extra-small mesh. In this two-scale framework, it is assumed that the material behaviour inside and outside the shear band is homogeneous. However, neither in real soil nor in a numerically modelled biaxial test the material behaviour inside the shear band is uniform throughout. In particular, the mechanical response inside the band is heterogeneous, and the band thickness evolves with deformation. Hence, hierarchical multiscale models such as FE2-method 37 and FEM-DEM methods 36 may provide a better representation of grain-scale micro mechanisms. Nonetheless, these methods are computationally demanding and require parallel computational solution strategies. Instead, this method consists of simple domain decomposition with an embedded scale refinement. Hence, it is computationally more affordable and efficient.

| APPLICATION OF THE DIFFUSION SB MODEL
As an initial step to demonstrate the applicability of the proposed diffusion SB model for large-scale problems, upscaled biaxial compression tests are conducted with drained top boundary. The specimen dimensions are 1 × 2 m. Quadratic elements with reduced integration (CPE8RP) are utilised throughout. At top nodes 0.4-m displacement is applied to maintain the same ultimate strain as in Figure 8. Upscaled biaxial compression simulations with two mesh sizes are conducted for six loading rates.
In the extra-large mesh ( Figure 13A), the element dimensions are 0.5 × 0.25 m, the same as the calibrated element in Section 4. It is simulated with both the original NS model and the calibrated diffusion SB model (Case 2). Six different sets of parameters are used for the different loading rates in the diffusion SB model. When loading path of an individual material point becomes nonproportional, the shear band angle calibrated in Section 4 for plane strain biaxial compression should be corrected as −α + θ sb to take account of the angle of principle axis rotation α = 0:5 tan − 1 2 σ 0 xy =ðσ 0 yy −σ 0 xx Þ. As a benchmark case, the upscaled specimen is discretised with extra-small mesh (0.00625 m) and simulated with the original NS model ( Figure 13C). The mechanical responses of the extra-large mesh with both the original NS model and the diffusion SB model are compared with those of extra-small mesh with the original NS model.
As graphically illustrated in Figure 13, the material response of single bifurcated element in Figure 13A can be represented by an undrained extra-small mesh with a shear band ( Figure 13B). Nonbifurcated Gauss points should follow the original NS model. Figure 13 elucidates three scales: the shear band is in micro scale, the element size is in mesoscale and dimensions of the upscaled biaxial test are in macro scale. Although this upscaled biaxial test does not represent the actual field conditions, it can merely demonstrate the validity of the proposed approach for arbitrary remote boundary conditions.  Figure 15.

| Global response
Due to numerical singularities, simulations of the extra-small mesh are terminated prematurely. The termination occurred when two shear bands intersect, leading to extreme mesh distortions at the intersection as displayed in Figure 15A,D,G. The shear strain contours in Figure 15 show localisations only in the extra-small mesh. Deformations in the extra-large mesh are more or less homogeneous except for elements near the top surface. The diffusion SB model is activated only for Gauss points that have reached the bifurcation strains mentioned in Table 2 for respective rates. It is noted that the concentrated deformation is encouraged by the diffusion SB model although it does not fully capture the localised shear strains similar to the extrasmall mesh. Figure 14 depicts that the mechanical responses of the extra-large mesh modelled with the diffusion SB model closely replicate the behaviours of the extra-small mesh when the loading strain rates are 40 and 20 s −1 . The deviatoric strain contours of the extra-small mesh corresponding to those loading rates show a single set of crossshaped shear bands.
When the loading durations are smaller (or for faster loading rates), results of these static simulations are doubtful. The deviatoric strain contours of extra-small mesh in Figure 15A display pore pressure shocks (alternative dilative and contractive shear bands). Moreover, the total reaction forces in Figure 14A show uniform undrained deformations without mesh sensitivity. Therefore, a dynamic simulation with inertia terms will be more appropriate for these rates.
On the other hand, when the loading rates are slower, there is enough time for the fluid movements at internal nodes. In these cases, local volume changes may occur in the extra-large mesh contradicting the hypothesis of the undrained diffusion SB model. The global force-displacement relationships of loading strain rates 4, 2 and 0.2 s −1 have softening branches. It ascertains that the drainage takes place both internally and externally. Moreover, when the loading rate is 0.2 s −1 , deviatoric strain contour of extra-small mesh shows arrays of tributary localised bands near the top surface as displayed in Figure 15G. These tributary shear bands facilitate further drainage at the top boundary. Therefore, for these loading rates, the onset of bifurcation predicted by the diffusion SB model (which is calibrated for an undrained element) is delayed than the predictions of the extra-small mesh. In other words, the diffusion SB model overestimates the peaks of extra-small mesh for slower loading rates. It is evident that the applicability of the diffusion SB model depends on the loading rate, permeability, drainage length as well as the remote boundary conditions. Nevertheless, as compared with simulations of extra-large mesh with the original NS model, the error from the proposed method is marginal.

| Average elemental response
Apart from the material behaviour, the global response of the biaxial compression test is also dominated by remote boundary constraints. Hence, global performance itself is not adequate to conclude about the validity of the proposed method. Therefore, it is worthwhile to compare the average constitutive behaviour of each element (four Gauss points) in the extra-large mesh with the average response of a respective subset of elements in the extra-small mesh. The purpose is to check whether the proposed constitutive model captures the average mechanical response of the extra-small mesh locally.
Only results of loading strain rate 20 s −1 are considered for the discussion below. The extra-small mesh is partitioned such that each partition is equivalent to a respective element numbers in the extra-large mesh as shown in Figure 13. The average mechanical responses of partitions in the extra-small mesh are compared with those of equivalent elements in the extra-large mesh.
The progressive development of two crossed shear bands in the extra-small mesh ( Figure 15D) can be described using Figure 16B, which encloses the stress-strain responses of each partition. In the extra-small mesh, there are only eight partitions that contain shear bands (1, 4, 6, 7, 13, 10, 11 and 16). Out of them, partitions 1 and 13 have drained top boundaries and, hence, shear bands initiate in them. Therefore, these two partitions inherit the highest deformation as shown in Figure 16B. Also, their ultimate strength is smaller because there is no volume constraint. These shear bands then progress to partitions 6 and 10; thus, they have second highest deformation. With further axial displacement on the top, these shear bands approach to partitions 7/11 and then 4/16 successively. Their deformations are also in the successive order. The average strength increases when partitions are further away from the drained boundary.
The partitions without shear bands (9,5,3,15,12,8,14 and 2) display slight unloading after the bifurcation. Nevertheless, a continued post-localised deformation could not be observed for the extra-small mesh. Extreme mesh distortions at the intersection of two shear bands cause convergence difficulties resulting in termination around 65% of applied displacement. From the overview of the average behaviours of partitions in the extra-small mesh, it is perceived that partitions with shear bands show enhanced deformation compared with others.
On the contrary, the extra-large mesh with the original NS model can hardly capture this progressive failure. No visible shear band is seen in Figure 15E, and all elements except those in the first two rows show a continuous dilative hardening in Figure 16A. Also, as the distance between two nodes is comparatively high in the extra-large mesh, the boundary movement has more influence on elements in the neighbourhood. Therefore, they tend to show considerable deformation than their respective extra-small mesh partitions. Specifically, elements in the extra-large mesh are unable to apprehend the propagation of shear bands and elastic unloading.
When the diffusion SB model is applied to the extra-large mesh, a considerable improvement in capturing the postlocalised deformation can be seen. Figure 16C illustrates the average stress-strain relationships of elements in the extra-large mesh simulated by the diffusion SB model. It is evident that the stress-strain relationships of elements with shear band simulated by diffusion SB model closely replicate their respective partitions of the extra-small mesh. However, elements in the first two rows of the extra-large mesh apparently undergo greater deformation than their respective partitions of the extra-small mesh. This is because the influence of drained top boundary is higher when the mesh size is larger, and these elements undergo net volumetric strain contradicting the hypothesis of the diffusion SB model. Predictions are satisfactory for elements away from the moving boundary (16, 4, 11 and 7).
In Figure 16B, the deformations of partitions 9, 5, 8, 12, 2, 14, 15 and 3 in the extra-small mesh abruptly terminate after the bifurcation starts. It is because all elements in those partitions are unloading, and their average responses show slight snapbacks. If the extra-small mesh had converged fully, further reloading and unloading could be expected depending on the development of secondary shear bands. Although the diffusion SB model captures this diminished deformation for elements 12, 8, 15 and 3, it inaccurately predicts continuous shear deformations for elements 9, 5, 14 and 2. This is because the impact of permeable and moving boundary causes elements in the first two rows to deform considerably with volumetric changes. As the onset of the shear band is detected based on the deviatoric strain, the diffusion SB model cannot distinguish the difference between concentrated and uniform deformation. Hence, the SB model wrongly assumes that internal shear bands have occurred in those elements. As the proposed method is calibrated for an element with zero volumetric strain, it cannot be recommended for elements close to moving or permeable boundaries.
In conclusion, the diffusion SB model is competent in reducing the stiffness caused by local drainage and correcting the stress increment when the strain increment is given, for bifurcated elements. However, the method used to detect the onset of the shear band is ambiguous for arbitrary boundary conditions. Because the correction is only applied at the constitutive level, the simulation has all the drawbacks of finite element method. The upscaled model with 16 elements is probably not appropriate to assess the validity of the proposed approach. One reason is that there is a considerable impact from moving boundaries. The other is that due to incompressible material behaviour at higher loading rates, larger elements do not perform well.

| Range of applicability
It is evident that the applicability of diffusion SB model is governed by the degree of internal drainage, which depends on boundary constraints, drainage length, loading rate and soil diffusivity. Figure 17 plots the maximum reaction forces achieved within 20% strain by the extra-large and extra-small meshes. The normalised velocity is calculated as VH c v where V is loading velocity and H is the drainage length, which is assumed as the height of the specimen. Presuming that the extra-small mesh predicts the ground truth, Figure 17 also includes the errors of predicting maximum forces by the extra-large mesh with and without the diffusion SB model. The loading strain rate 200 s −1 is ignored here because it closely represents the locally undrained case. It is depicted that when mesh size is increased nearly 56 times, the error of predicting peak is between 61% to 195% by the original NS model. When the diffusion SB model is used with extralarge mesh, this error is between 8% to 52%.
The accuracy of the diffusion SB model is within 20% when the normalised velocity is between 10 to 40 for the considered upscaled analysis with a permeable top boundary. In this velocity region, the mesh sensitivity of the original NS model is reduced by more than 50% by the diffusion SB model. However, it should be emphasised that these values are dependent on the considered specimen size, loading and drainage conditions.

| LIMITATIONS OF MODEL AND PROPOSED EXTENSIONS
The proposed approach is a heavily simplified mathematical model that takes account of internal fluid movements within the constitutive level without changing equilibrium equations. Hence, there are several limitations inherent to this model.
First, because the diffusion SB model is designed for locally undrained (isochoric) material points, it is not valid for Gauss points near permeable boundaries or when the loading rate is slow enough to induce local volume changes. As shown in Section 5.3, its applicability is limited to a certain degree of drainage. Second, the model is calibrated with a numerical finite element test with specific boundary conditions (isochoric biaxial compression). Hence, input parameters, such as shear band angle, the band thickness, onset of the shear band and shear factor, are calibrated for the average material behaviour of biaxial specimen intercepted by a single diagonal shear band. Because these are not intrinsic material parameters, they are only valid for the considered loading path and boundary conditions. Further, scale effects may also affect when parameters calibrated at mesoscale are used for a macro-scale boundary value problem. Third, this model is calibrated for a single diagonal shear band only. Hence, it cannot replicate intricate localisation patterns with tributary and curvilinear bands. Also, change in the loading path can deactivate an existing localisation band and activate another, which cannot be recognised by this model.
Fourth, because the diffusion SB model is validated with another finite element model with a finer mesh, the reliability of the 'ground truth' itself is questionable. By default, it has all the drawbacks of the finite element method, including mesh sensitivity. Numerical singularities can arise, such as structural snapback corresponding to zero energy dissipation. Hence, there is a minimum shear band thickness (smallest mesh size) this model can be calibrated with. So the proposed model is not appropriate for soil with very fine particles.
In an arbitrary condition, the formation of locally drained shear band depends on both micro characteristics of soil (material heterogeneity, grain size, viscosity, permeability) as well as macro properties of the problem (loading path, boundary conditions, problem geometry and dimensions). The conventional bifurcation analysis is incompetent to recognise this kind of instability. Mostly, imperfections of geometry or boundary condition may trigger local drainage, which will ultimately develop into material instability. Thus, the onset of instability is problem specific.
Therefore, if this model is applied for a general boundary value problem, a practical approach to detect the bifurcation condition is necessary. Vardoulakis, 6 Guo 7 and Guo and Stolle 8 theoretically proved that the onset of locally drained shear bands under globally undrained condition is more likely to occur during the softening regime of the mobilised friction angle (following the peak effective stress ratio) because in this region, the jump in local volume change required to trigger a dilative shear band is smaller. As a conservative approach, the criterion for the onset of local drainage can be set as the peak effective stress ratio. The angle of localisation can be determined from the Coulomb orientation (modified for the principal stress rotation). The shear band thickness can be related to the mean grain size of the soil. This method is applicable to any general loading path and boundary condition. However, when compared with the global response of numerical biaxial compression test, this criterion prepredicts the onset of the shear band. Even based on the experimental evidence, locally drained shear bands occur much later to the Coulomb criterion. 29 The proposed constitutive model is extendable for multiple intersecting shear bands. Conceptually, a nested scale model can be developed enriching the zone outside the first localised band with new bands. However, this will be very complicated because it is difficult to predict their onset and propagation. Those details are sensitive to the microstructure of the soil. FEM-DEM approach, as in Guo and Zhao, 36 is more appropriate. Therefore, as far as the current implementation is concerned, the proposed approach is valid only for discretisation, which is sufficiently fine to capture the propagation of the major shear band.

| CONCLUSIONS
A phenomenological constitutive model (diffusion SB model) is presented to apprehend the undrained deformation of saturated dilative sand in the presence of a locally drained shear band. Rate dependence of hydromechanical coupling at the shear band scale is captured through the diffusion equation. Separate constitutive relationships are generated for the shear band and intact material, which are averaged to calculate the macroscopic response. The model can be used at the constitutive level without changing equilibrium equations.
Conventional bifurcation analysis fails to determine the onset of locally drained shear bands within globally undrained condition because it is not an intrinsic material instability. This kind of localisation is affected by soil microstructure and imperfections of the specified problem. Hence, it cannot be recognised by a continuum-based constitutive relationship of the material. Therefore, parameters related to the onset of localisation are derived from a numerical experiment.
First, the proposed diffusion SB model is calibrated with the global response of the extra-small mesh with impermeable boundaries. A small material heterogeneity is added to trigger the local drainage. It is pointed out that undrained dilative hardening predicted by a homogeneous material model is reduced due to the local drainage associated with shear localisation. The individual material behaviours inside and outside the band predicted by the diffusion SB model are qualitatively in good agreement with those of the extra-small mesh. It is illustrated that depending on the loading rate and the shear band thickness, the behaviour inside the band can be fully or partially drained. This degree of drainage, in fact, governs the strength.
Then, the diffusion SB model is applied for an upscaled biaxial compression test, and its accuracy is compared with predictions of an extra-small mesh with the original NS model. Based on upscaled biaxial compression results, it is apparent that applicability of the diffusion SB model depends on the loading rate, permeability, drainage length as well as the remote boundary conditions.
For the considered boundary value problem, the proposed diffusion SB model can reproduce the global reaction forces of a 56 times finer mesh with 20% accuracy limit when the normalised velocity is between 10 and 40. Within this region, the diffusion SB model reduces the mesh sensitivity of conventional NS model by more than 50%. However,these results are problem specific. When the proposed diffusion SB model is applied to a general boundary problem, a more practical approach is necessary to detect the onset of local drainage.
ORCID Hansini Mallikarachchi https://orcid.org/0000-0002-0140-9145 dp i dϵ p where h is the hardening or softening modulus and p i,max is the maximum image pressure, which is a function of state parameter at the image state Ψ i .
The critical void ratio at image state (e c ) i is defined by the relationship where e max and e min are maximum and minimum void ratios, respectively. The hardening of the yield surface is restricted to match the maximum dilatancy of the real sand. According to experimental data, a linear relationship can be obtained between the maximum dilatancy D max and Ψ i .
where χ tc and M tc are the maximum dilatancy coefficient and critical stress ratio obtained from triaxial compression data. This maximum dilatancy is used to size the yield surface so that the real soil behaviour matches with the dilatancy from the normality. Hence, the maximum yield surface size can be calculated from Due to this limiting hardness concept, the NS is capable of properly capturing the soil behaviour in the dry side of the critical state even with an associative flow rule.
The shear modulus G e is calculated as where A is shear coefficient and r is pressure exponent.