Numerical implementation of BBM in FE package for solving unsaturated soil boundary value problems

. Unsaturated soil mechanics is of keen relevance when dealing with soil above the water table or compacted soils. Barcelona Basic Model (BBM) is one of the most widely used constitutive models for describing unsaturated soil behaviour and is available as a soil model in software like PLAXIS and CODE_BRIGHT (developed by UPC). The present study demonstrates a numerical framework to implement BBM within the ABAQUS package as a user-defined soil model using a user subroutine called UMAT (User Material). The backward Euler integration scheme coupled with Newton Raphson iterative algorithm is used to construct the model. A convergence test is conducted to check the accuracy and stability of the integration. UMAT's inability to manage pore pressure led to proxy diffusion to depict suction variation. The heat equation is used as a proxy for diffusion in three-dimensional space and time. An oedometer test is simulated in which the soil is allowed to consolidate and then set to saturate at constant vertical stress, which leads to wetting collapse. The results demonstrate a reasonable behaviour of wetting-induced compression of soil subjected to loading with varying saturation level.


Introduction
The paradigm of unsaturated soil mechanics has been around for quite a while and has proved to be a valuable approach for understanding and describing the engineering behaviour of soils, especially those present above the water table and the compacted ones. The classic soil mechanics let the degree of saturation to move from zero to hundred percent, i.e., the soil mass varies from being totally dry to totally wet, where the void will contain only air and only water/fluid respectively. The main scenarios of practicality are not met by the provided assumption of the traditional soil mechanics. Realistically, the soil mass cannot be entirely saturated since there will always be some percentage of air-filled voids, and even if the same soil is expected to dry out, some hygroscopic water will still be present in the soil mass, preventing complete drying. To comprehend this compressive behaviour of unsaturated soils, several studies have been conducted [1][2][3][4]. One of the most critical phenomena involved in compressive behaviour is wetting-induced collapse. Upon saturation, soil tends to collapse gradually or suddenly, which depends on the factors like the degree of saturation, pre-consolidation stress and the state of stress. Various approaches were suggested to understand the collapse response of the soil [5][6][7] by introducing a term, collapse potential, which was correlated with initial dry density, vertical stress at the wetting path, clay fraction, initial water content and plasticity index (PI). Furthermore, constitutive models were developed by [7][8][9][10][11][12] and used to simulate routine * Corresponding author: arghya@iitk.ac.in laboratory tests, but they were hardly ever used in computer programs because of their model complexities.
Barcelona basic model (BBM) is a widely used unsaturated soil constitutive model that has been used extensively in the area of unsaturated soil mechanics. The BBM was originally introduced by [13] and was further implemented in computer programs to solve boundary value problems [14][15][16][17][18]. Generally, the model was implemented in the software Code_BRIGHT (developed by UPC), finite difference software FLAC, ALLFINE. The present study demonstrates a numerical framework to implement BBM within the ABAQUS package as a user-defined soil model using a user subroutine called UMAT (User Material). The backward Euler integration scheme coupled with Newton Raphson iterative algorithm is used to construct the model.

Barcelona Basic Model (BBM)
BBM is a hydromechanical elastoplastic constitutive model and describes the stress strain relationship of the soil behaviour. The model has been introduced in the space of deviatoric stress ( ), mean net stress ( ) and suction = − . Initially, it was devised in the ( − ) plane. To extend its modeling capabilities into triaxial stress state, it was coupled with modified cam clay model (MCC) in such a way that for the suction of zero value i.e., full saturation, the BBM converts back into MCC. This analysis is based on two independent variables namely net stress − and suction = − where is air pressure and is water pressure and is the stress state of the system. Here, two assumptions are made i.e., (1) Pore air pressure ( ) is zero making the net stress equal to the total stress, which can be generally observed in the field [9] and (2) The stiffness and strength during the wetting process are impacted by matric suction.

Yield Surfaces
The corresponding yield equations for the BBM model is given by Eqs. 1-4.
Here, is preconsolidation stress for unsaturated conditions and governs the compression yielding; is increase in cohesion with suction and governs the tensile yielding of the constitutive behaviour; is the reference pressure, (0) and ( ) is nondimensional stiffness parameter at = 0 and ≠ 0 respectively; is the suction increases, the maximum suction soil has ever witnessed.

Flow Rule and Hardening Law
Non-associated flow rule is used to govern the LC (loading collapse) curve (Eqs. 5), whereas for the SI (suction increase) curve, associated flow rule is implemented (Eqs. 6).
Here, is the deviatoric stress multiplier. Parameters , , , * , contribute to the behaviour of the MCC part of the model, where is the slope of the critical state line and is the shear modulus. , and controls the evolution of the LC curve with suction, where is a material parameter controlling the rate of increase of soil stiffness with suction; is the parameter governing the asymptotic maximum stiffness.
is the matric suction cap and is rate of increase in cohesion with suction. and corresponds to the nondimensional bulk modulus and stiffness parameter for suction changes in virgin state of the soil. This makes total of twelve parameters which are used in this formulation. As the yielding approaches, the hardening parameters and changes in accordance with the total volumetric deformation i.e., the sum of volumetric deformation contributed from both suction and mechanical loading.

Calibration and Numerical Integration
The constitutive models are simply a combination of the various differential equations. In order to get the total stress-strain behaviour using these incremental models, they need to be solved by numerical methods. In this approach, the formulation has been done using the implicit method (predictor-corrector approach) coupled with the newton Raphson iterative scheme. Newton Raphson iterative scheme is helpful in achieving better accuracy and this method is generally called fully implicit formulation.

Formulation
In this formulation, the integration assumes the strains to be the governing state variable in the mechanical loading whereas the suction for the case of the hydraulic loading. The total strain and suction to be induced is initialized, and using an incremental value, each iteration checks for the yield surface. If the state of the system is within the yield zone, the elastic law applies; otherwise, the yield surface updates and leads to plastic deformation.
The supplied incremental ̇ will contribute to incremental stress ̇. The trial stress (̇+̇) will confirm whether the current state is within the yield zone or not. If lying within, the next iteration comes contributing to stress and checks for yield. If lying outside, the hardening parameters and stresses get updated and yield surface shifts accordingly.
is the compliance matrix defining the stress-strain relationship (̇=

Convergence Test
In order to verify the accuracy and stability of the numerical stress integration of the proposed formulation, a convergence test is conducted with various step sizes of suction values, i.e., 02 kPa, 04 kPa, 08 kPa and 16 kPa. The convergence result shows a very satisfactory accuracy as can be seen in the Fig. 4.

Validation of Suction
In this work, the BBM model is implemented in the FE package ABAQUS as a user defined model using the user subroutine UMAT (user material). It can be thought of as a development facility provided by ABAQUS to incorporate a material behaviour not present in the inbuilt library. In general, the UMAT is to do update in the stress and corresponding state variables occurring due to the incremental strain supplied by the ABAQUS. The inability of UMAT to control the pore pressure led to the requirement of a proxy law to supply the input from ABAQUS to suction state variable defined in each step. To fulfil this obligation heat equation is chosen and combined with the UMAT code developed. The heat equation considers/works in two independent spaces, i.e., three-dimensional space and time. The power over the three-dimensional space is manageable since it will be within the obligation of the dimension of the problem defined. However, the control over the time is not that easy. Hence, to device the time regulation, an analogy is created between the 1D consolidation and 1D heat equation. Table 1 represents the side-to-side comparison of 1D consolidation and 1D heat equations and proposed analogy between the variables.
The temperature supplied by ABQUS is treated as suction by specifying the initial suction value as the initial temperature in a predefined field. Furthermore, 0 20 40 60 80 100 120 the suction value to be simulated in the respective step is defined by temperature over the set or surface. The thermal conductivity, specific heat and mass density parameters in the heat law works for the hydraulic conductivity, volume compressibility and specific weight of water of consolidation law, respectively. In the absence of the hydraulic conductivity, it needs to be calibrated by hit and trial until the results fulfil the necessary requirements. To validate the working of this approach of utilizing heat law as a proxy to simulate the suction variation, a test has been conducted for which the material properties is given by Table 2. A fully saturated cube sample ( =0 kPa) is subject to the suction variation from 0 kPa to 100 kPa (path A-C) i.e., drying followed by wetting from 100 kPa to 0 kPa (path C-D) and finally drying till 1000 kPa (path D-E). Fig. 5(a) shows a stress path due to suctioninduced loading on clay paste, adopted from [19]. Point A0 denotes the initial condition at which the clay sample is fully saturated without any stress history, i.e., * and were zero. Resulting from vertical loading of 25 kPa, the hardening took place in both SI and LC, which led to a value of 16.47 kPa and * increased to 22.24 kPa, which is denoted as point A. The clay specimen was dried in the path A-B-C, having an elastic zone (AB) up to the suction value of 16.47 kPa, after which the elastoplastic zone (BC) commenced activating the SI. At 100 kPa, the soil was subjected to wetting/redrying in the elastic domain and followed by continuous drying till the suction reached 1000 kPa. Numerical results in the plastic domain show lower deformation than the experimental results due to suction-dependent plastic compressibility ( ), which is nonlinear in reality.

Oedometer Test
Previous result demonstrated the constitutive response of the model, whereas to establish the ability to phenomenon of wetting collapse, an oedometer is simulated and presented. The parameters used for the simulation are given in Table 3. The rectangular sample size of height 200 mm and width 100 mm is used to simulate the oedometer test. Initially, the lateral sides are constrained in the x direction, whereas the bottom of the sample is constrained in vertical movement, i.e., in the y direction. The top surface is set to saturate in the second step, having an initial suction of 0.2 MPa and saturated preconsolidation stress is * = 0.075 MPa. The intrinsic permeability of the sample is 5*10 -13 m 2 . The unsaturated sample is first loaded till the vertical load reaches 0.2 MPa, followed by wetting till the saturation (Path AB) at constant vertical stress of 0.2 MPa and further loading till 0.4 MPa (Fig. 6). The BBM model is assumed to expand when suction increases, whereas the model contracts when suction is decreased. Thus, as the saturation approaches, if the mean net stress is large enough, the result comprises compressive deformation and leads to a reduction in the void ratio. Zone OA, zone AB, and zone BC represent unsaturated, wetting-induced collapse, and saturated zone, respectively.    Fig. 9, which shows the gradual saturation of the sample over time.

Conclusions
The BBM model is implemented in ABAQUS as a user defined material using the heat law as a proxy to simulate the suction variable. The temperature input for the heat flow is stored in the suction state variable to conduct the simulation. In this approach, the suction formulation is considered separately from the mechanical formulation except for the hardening behaviour. The hardening behaviour is deemed to be dependent on both the hydraulic and the mechanical induced deformation.
The formulation is written in the ABAQUS subroutine UMAT in the language of FORTRAN and successfully validated for the suction induced deformation. The simulation is conducted to establish the volumetric behaviour (which can be described by void ratio) by varying the suction in both drying and wetting paths. The resulting plot has been validated with that present in the literature [21].
Finally, an oedometer test is conducted where a sample of size 100mm diameter and 200mm length is subjected to vertical load. At constant vertical load, the saturation is varied toward saturation. The results demonstrated the wetting induced collapse phenomenon. The outcoming results are satisfying reasonably with the oedometer test results by [20].