Numerical Implementation of the Barcelona Basic Model Based on Return-Mapping Integration

Abstract

This paper implemented the Barcelona basic model (BBM) into the OpenGeoSys (OGS) platform for numerical modeling of the coupled hydro-mechanical (HM) behavior of unsaturated soil. Within the implicit integration approach in the OGS, the integration rule of the BBM was developed first and the integration form of the BBM under a return mapping algorithm was built. The closest point projection method was used for calculating the return mapping directions with the associative flow rule. The numerical simulation results show that the BBM is feasible in fitting the experimental results. The numerical integration algorithm can reflect the elastic–plastic mechanical behavior of materials, and improve the calculation accuracy. The material exhibits obvious elastic–plastic characteristics during numerical simulation and experiment, and the water absorption process can lessen the mixture’s compression stiffness while enhancing its recovery stiffness.

1. Introduction

The hydro-mechanical (HM) performance of unsaturated swelling clay soil is a serious concern for many engineering applications including infrastructure, transportation, and engineering buffers [1,2,3]. Take the high-level radioactive waste (HLW) repository as an instance, the typical swelling material, bentonite–sand mixture, is preferred as the backfill or buffer material to seal the repository with an attempt to close the excavating fracture and isolate nuclide migration by its swelling behavior and low permeability [4,5,6]. Evidence from the laboratory suggests a complex HM coupling process in unsaturated swelling soil [7,8,9], and great effort has been devoted to the study of the constitutive model for unsaturated swelling soil. Among these constitutive models, the Barcelona basic model (BBM) and its expanded formulation, named the Barcelona expansive model (BExM), obtained great success in describing the features of the HM coupling processes in unsaturated swelling soil [10,11,12,13].

Numerical simulations are commonly used to analyze the coupling process of multiple physical fields, as carrying out a long-term experimental investigation is often not possible [14]. Many studies have focused on numerical simulations of the mechanical properties of bentonite-based materials during water absorption and expansion, and have made many achievements. By taking the ideal elastic model, the Drucker–Prager plastic model, and the wetting expansive model, the stress–strain curve of bentonite was described in a simulation using ABAQUS [15]. According to the finite element method, the nonlinear swelling pressure development of bentonite under saturation conditions is predicted [16]. Ballarini et al. [17] simulated and analyzed the thermo-hydro-mechanical (THM) coupling process of buffer materials by OpenGeoSys, and the sensitive parameters were identified, which has guiding significance for the quantitative selection of materials in construction. Teams from four countries coded with ABAQUS, FRACON, THAMES, and ROCMAS to conduct numerical simulation studies on the coupling process of bentonite-based buffer materials in the international cooperation project DECOVALEX II [18,19,20,21].

Since the BBM can describe many typical features of the mechanical behavior of unsaturated soils [22], including the swelling or collapse strain caused by wetting, the BBM is widely used in the simulation of HLW buffer materials in countries such as Europe and Japan [23]. Xu et al. [24] used BExM and implemented in TOUGHREACT-FLAC3D, which is a multi-phase reactive transport and geomechanics simulator, to simulate the coupled thermal-hydraulic-mechanical-chemical (THMC) processes, with a focus on the influence of chemical changes on mechanical process over time. Lee et al. [25] employed the thermo-elasto-plastic version of the BBM to characterize the mechanical behavior for the numerical simulations of full-scale engineered barriers experiment (FEBEX) and obtained approximate results in simulations including temperature, relative humidity, total stress, and displacement.

Although research related to the mechanical constitutive model of the coupling process was widely performed, studies that evaluated the influence of integral algorithms on numerical simulation are limited. The integration method of constitutive equation is directly related to the accuracy and stability of calculation [26,27,28]. Different from the general unsaturated soil constitutive model, the integration algorithm of expansive unsaturated soil is more complicated due to the consideration of plastic characteristics and the effect of plastic deformation on water retention characteristics [29,30]. At present, implicit integration algorithms [31,32] and explicit integration algorithms [33] are commonly used in the field of numerical simulation. TOUGHREACT-FLAC3D uses an explicit integration algorithm, which does not require iteration and has high calculation efficiency. However, it does not strictly meet the consistency conditions in the integration process, which ultimately leads to the stress point possibly deviating from the yield surface, resulting in error accumulation and limiting the algorithm’s accuracy [32,34,35]. In contrast, the implicit integration algorithm corrects plastic stress using a nonlinear equation iterative method, which is widely used due to its high calculation accuracy and efficiency [33,36]. The application of the return mapping algorithm allows the deviation of the initial stress to the yield surface to be predicted elastically, and the implicit integration algorithm can correct the deviated stress plastically and pull the deviated stress back to the yield surface to complete the calculation [37,38,39,40].

In this paper, an implicit algorithm is used to study the numerical calculation of the elastic–plastic mechanical model of expansive unsaturated soil. The elastic–plastic mechanical response of bentonite/quartz sand mixture in the HM coupling process is analyzed by numerical simulation, and the applicability of the theoretical model is verified by experiments. The study of the elastic–plastic mechanical properties of bentonite/quartz sand mixtures in the HM coupling process of the HLW repository can provide theoretical and model support for the study of HM coupling and other multi-physical field coupling processes in the geological disposal of HLW and the implementation of engineering. At the same time, it has a certain guiding significance for the future application of this kind of materials in geological engineering and geotechnical slope engineering.

2. Implicit Integration of the BBM

As a traditional model to characterize the elastic–plastic mechanical behavior of unsaturated soil in geotechnical engineering, the BBM was developed based on the CAM clay model, proposed by Alonso et al. [41]. Based on plasticity theory, the critical state model, and the modified CAM clay model, the BBM is suitable for describing the mechanical behavior of unsaturated soils, expansive clay, low plastic sand, silty clay, and other materials [42,43]. Taking pore suction as a variable, the BBM was used to study the influence of soil suction on soil elastic–plastic mechanical behavior, and can be used to characterize the shear strength and consolidation brought on by hydroscopic expansion, collapsibility, and pore suction change in unsaturated rock and soil [44]. The BBM, like other elasto-plastic constitutive models, divides soil deformation into elastic and plastic deformation, employs three stress variables $(p,q,s)$, and employs yield surfaces $(LC,SI,CLS)$ to divide stress space into plastic and elastic space. Figure 1 presents the three-dimensional representation of the yield surface in the thermo-elasto-plastic BBM three-dimensional yield surface in $p-q-s$ space and $p^{\prime }-q-T$ space, where $p^{\prime }$ is the net mean stress (i.e., total stress minus gas-phase pressure), q is the deviatoric stress (or shear stress), s is suction, and T is temperature [45].

Due to the low permeability of bentonite materials, the seepage–stress coupling time is long during the encapsulation process [46]. The numerical calculation time step is an important factor affecting the accuracy of the calculation results and the time cost. Therefore, the implicit integration algorithm based on the return mapping is used to establish the numerical integration expressions of the BBM, which is embedded in the OpenGeoSys numerical simulation software to reduce the calculation steps and improve the calculation accuracy. In this section, the equations of the BBM were introduced firstly, afterward, the integration form of the BBM was developed according to the approach of Borja and Lee [47]. The consistent tangent modulus is a very important part of the constitutive model integration, and the determination of Tangent Modulus of Stress–Strain Curve is given in Appendix A.

2.1. Equations of the BBM

The stress state in the BBM is described by the net mean stress in the form of

$${\sigma }^{\prime }=\sigma -u_a\mathit{I}$$

(1)

in which the volumetric stress p and deviatoric stress q are defined as

$$p=\frac{1}{3}tr\left({\sigma }^{\prime }\right),q=\sqrt{\frac{3}{2}}\parallel \xi \parallel ,$$

(2)

here, $\xi ={\sigma }^{\prime }-1/3tr\left({\sigma }^{\prime }\right)\mathit{I}$.

The elastic strain of volumetric and shear components are expressed as

$$d{ϵ}_v^e=\frac{\kappa }{v}\frac{dp}{p}+\frac{{\kappa }_s}{v}\frac{ds}{\left(s+p_{a}t\right)}$$

(3)

$$d{ϵ}_s^e=\frac{dq}{3G}$$

(4)

It is notable that the BBM has two yield surfaces in the form of

$$f_1\left(p,q,s,p_0^{*}\right)\equiv q^2-M^2\left(p+p_s\right)\left(p_0-p\right)=0$$

(5)

and

$$f_2\left(s,s_0\right)\equiv s-s_0=0$$

(6)

where

$$p_s=k_{s}s$$

(7)

$$\frac{p_0}{p^c}={\left(\frac{p_0^{*}}{p^c}\right)}^{\left[\lambda \left(0\right)-\kappa \right]/\left[\lambda \left(s\right)-\kappa \right]}$$

(8)

$$\lambda \left(s\right)=\lambda \left(0\right)\left[\left(1-r\right)exp\left(-\beta s\right)+r\right]$$

(9)

The yield surface of the BBM is determined by the hardening parameters $p_0^{*}$ and $s_0$, which depend on the total plastic volumetric strain increment $d{ϵ}_v^p$ in the form of

$$\frac{d{p}_0^{*}}{p_0^{*}}=\frac{v}{\lambda \left(0\right)-\kappa }d{ϵ}_v^p$$

(10)

$$\frac{d{s}_0}{s_0+p_{at}}=\frac{v}{{\lambda }_s-{\kappa }_s}d{ϵ}_v^p$$

(11)

2.2. Equations of Associative Flow Rule

The evolution of plastic strain can be expressed with the associative flow rule in the form of

$${\dot{\epsilon }}_{vm}^p=\dot{\varphi }\frac{\partial f_1}{\partial {\sigma }^{\prime }}$$

(12)

where $\dot{\varphi }$ is a factor. Submitting Equation (2) into Equation (12), we get

$$\Delta {ϵ}_v^p=\Delta \varphi tr\left(\frac{\partial f_1}{\partial {\sigma }^{\prime }}\right)=\Delta \varphi \frac{\partial f_1}{\partial p}$$

(13)

The term of $\partial f_1/\partial {\sigma }^{\prime }$ can be obtained by the chain rule as Borja and Kavazanjian [48]

$$\frac{\partial f_1}{\partial {\sigma }^{\prime }}=\frac{\partial f_1}{\partial p}\frac{\partial p}{\partial {\sigma }^{\prime }}+\frac{\partial f_1}{\partial q}\frac{\partial q}{\partial {\sigma }^{\prime }}=\frac{1}{3}\left(\frac{\partial f_1}{\partial p}\right)\mathit{I}+\sqrt{\frac{3}{2}}\left(\frac{\partial f_1}{\partial q}\right)\widehat{\mathit{n}}$$

(14)

where $\widehat{\mathit{n}}$=$\xi /\parallel \xi \parallel$. The derivation of $f_1$ respecting to p, q, $p_0$, and $p_s$ can be deduced from Equation (5) as

$$\frac{\partial f_1}{\partial p}=M^2\left(2p+p_s-p_0\right)$$

(15)

$$\frac{\partial f_1}{\partial q}=2q$$

(16)

$$\frac{\partial f_1}{\partial p_s}=M^2\left(p-p_0\right)$$

(17)

$$\frac{\partial f_1}{\partial p_0}=M^2\left(p+p_s\right)$$

(18)

2.3. Closest Point Projection of the Constitutive Relation

The returning mapping tensor of stress can be defined as

$${\sigma }_{n+1}^k={\sigma }_{n+1}^{tr}-{\mathbb{C}}^e:\delta {\epsilon }^p$$

(19)

where ${\sigma }_{n+1}^{tr}$ is the trial stress. For the volumetric and deviatoric parts of ${\sigma }_{n+1}^k$, which has a formulation of

$$p\equiv p_{n+1}^k=\frac{1}{3}tr\left({\sigma }_{n+1}^k\right)=p_{n+1}^{tr}-K\Delta {ϵ}_v^p$$

(20)

$$q\equiv q_{n+1}^k=\sqrt{\frac{3}{2}}\parallel {\xi }_{n+1}^k\parallel =\sqrt{\frac{3}{2}}\parallel {\xi }_{n+1}^{tr}-2\mu \Delta {\gamma }^p\parallel$$

(21)

where K and $\mu$ are the elastic moduli defined as

$$K=\frac{1+e}{\kappa }p$$

(22)

$$\mu =\frac{3K\left(1-2\upsilon \right)}{2\left(1+\upsilon \right)}$$

(23)

According to Borja and Lee [47], p and q can be determined by

$$p=p_{n+1}^{tr}-K\Delta \varphi \frac{\partial f_1}{\partial p}=p_{n+1}^{tr}-K\Delta \varphi M^2\left(2p+p_s-p_0\right)$$

(24)

$$q=q_{n+1}^{tr}-3\mu \Delta \varphi \frac{\partial f_1}{\partial q}=q_{n+1}^{tr}-6\mu \Delta \varphi q$$

(25)

The updation of hardening parameters $p_0^{*}$ and $s_0$ is complemented by the hardening equations of Equations (10) and (11). The integration of Equations (10) and (11) over a finite time increment yields

$$\begin{array}{cc}\hfill {\left(p_0^{*}\right)}_{n+1}& ={\left(p_0^{*}\right)}_{n}exp({\vartheta }^{*}\Delta {ϵ}_v^p)\hfill \\ & ={\left(p_0^{*}\right)}_{n}exp({\vartheta }^{*}\Delta \varphi \frac{\partial f_1}{\partial p})\hfill \\ & ={\left(p_0^{*}\right)}_{n}exp\{{\vartheta }^{*}\Delta \varphi M^2\left(2p+p_s-p_0\right)\}\hfill \end{array}$$

(26)

$$\begin{array}{cc}\hfill {\left(s_0\right)}_{n+1}+p_{at}& =\left({\left(s_0\right)}_n+p_{at}\right)exp({\vartheta }^s\Delta {ϵ}_v^p)\hfill \\ \hfill & =\left({\left(s_0\right)}_n+p_{at}\right)exp({\vartheta }^s\Delta \varphi \frac{\partial f_1}{\partial p})\hfill \\ \hfill & =\left({\left(s_0\right)}_n+p_{at}\right)exp\{{\vartheta }^s\Delta \varphi M^2\left(2p+p_s-p_0\right)\}\hfill \end{array}$$

(27)

in which

$${\vartheta }^{*}=\frac{v}{\lambda \left(0\right)-\kappa }$$

(28)

$${\vartheta }^s=\frac{v}{{\lambda }_s-{\kappa }_s}$$

(29)

$$\int (\nabla q)wd\Omega =n\int (\nabla qw)d\Omega -n\int \nabla q\nabla wd\Omega$$

(30)

3. Numerical Implementation

To validate the proposed numerical algorithm, this section builds a numerical model based on compression experiments of bentonite/quartz sand mixture. The elastic–plastic mechanical model of expansive unsaturated soil embedded in the OGS was used in the simulation. The THMC processes in porous media are numerically simulated using the finite element approach, which is based on a flexible numerical framework [49,50]. Using the BBM, the elastic–plastic mechanical behavior of a bentonite/quartz sand mixture during compression was simulated, and the findings of the numerical simulation were compared with the experimental results.

3.1. Experiment Benchmark

According to the compression experiment of MX-80 bentonite/quartz sand mixture completed by Wang et al. [51], the ratio of MX80 bentonite/quartz sand mixture used in the experiment was 70:30, and the sample suction was adjusted by dialysis method. Three standard samples with diameters of 38 mm and heights of 10 mm were prepared in the experiment and named SO-02, SO-03, and SO-04. In order to consider the influence of the gap between the sample and the container on the compression process, sample SO-01 with a diameter of 35.13 mm and a height of 10 mm was also prepared.

The specimens were placed in the dialysis device, as shown in Figure 2, and initially loaded 0.1 MPa pressure to ensure full contact between the specimens and the stress sensor, followed by a steam cycle to allow the specimens to reach a preset suction force. The initial dry density of samples SO-02~04 was 1670 kg/m${}^3$, and steam circulation was used to achieve suction forces of 4.2, 12.6, and 38 MPa. Before compression, the samples were allowed to expand fully until stability under the premise of maintaining the axial load of 0.1 MPa. The initial density of sample SO-01 was 1970 kg/m${}^3$, whose diameter was smaller than that in the dialysis container. After absorbing water and filling voids, the dry density of sample SO-01 was close to the initial dry density of the other three samples of 1670 kg/m${}^3$. Sample size and initial suction parameters in the compression test are shown in

Table 1. Physical parameters and dimensions of sample for suction controlled compression.

Tests	${\mathit{\rho }}_{\mathit{d}\mathbf{I}}$ (Mg/${\mathbf{m}}^3)$	D (mm)	s  (MPa)
SO-01	1.97	35.13	0.0
SO-02	1.67	38.00	4.2
SO-03	1.67	38.00	12.6
SO-04	1.67	38.00	38.0

.

After the specimen reached the preset suction, the axial compressive stress was gradually loaded up to 50 MPa by stress control at a loading rate of 0.5 MPa/h, and then we unloaded the stress at the same rate. The change of axial displacement with time was measured automatically by the experimental system during the experiment. Through calculation, the change data of volume strain with time in the process of experimental compression recovery is obtained, and then the change data of void ratio with time can be calculated. The relation between void ratio and volumetric strain is

$$K=\frac{1+e}{\kappa }p$$

(31)

3.2. Software Utilization

OGS is a scientific open-source project developed in the 1980s, which is dedicated to solving the coupling problems of temperature–seepage–stress and other physical fields in HLW packaging, geothermal development, energy storage, and other fields. It has been applied in many fields involving the coupling problems of temperature–seepage–stress chemistry and other physical fields. The integration algorithm of the unsaturated soil elastic–plastic model was used in the numerical simulation, and the software implementation of the BBM with implicit integration has been shown in Figure 3.

According to the experimental process described above, the numerical calculation model adopts an axisymmetric model. As shown in Figure 4 for specimen SO-01, the vertical displacement of the bottom side of the specimen was fixed, and set the bottom and side as the water inflow boundary condition until the model reached the preset suction. Set the water inflow boundary condition at the bottom of the model and fix the lateral displacement of the lateral boundary when the sample expands to the pore closes. The suction of the sample would reach to 0.01 MPa when saturated. When the water absorption of sample SO-02~04 reaches the preset suction, the vertical displacement of the bottom boundary and the lateral displacement of the side boundary of the model are fixed, and the bottom boundary of the model is set as the water inflow boundary, as shown in Figure 5. In the process of compression and recovery, the boundary conditions of samples SO-01~04 are the same. The axial pressure is gradually loaded to 50 MPa on the top boundary of the sample for 100 h and then unloaded the pressure within the same time. The time step for both compression and recovery is 1 hour. In the process of water absorption, the final axial strains of samples SO-01~04 are approximately 18%, 6.8%, 5.4%, and 1.2%, respectively. Therefore, in the compression experiment, the initial porosity of samples SO-01~04 is 0.97, 0.75, 0.73, and 0.69, respectively.

In the compression experiment under suction control, parameters related to the numerical calculation of water absorption and expansion process are shown in

Table 2. Numerical parameters for hydromechanical modeling of bentonite/sand mixture.

Mechanical Parameters
M (-)	1.07
$P_{os}$ (MPa)	0.66
$P_{ref}$ (MPa)	0.01
$\lambda$ (-)	0.1272
${\lambda }_s$ (-)	0.03
$kappa$	0.081
$kapp{a}_s$	0.008
r (-)	0.75
$\beta$ (MPa${}^{-1}$)	0.25
$e_o$ (-)	0.67
k (-)	0.1
$s_o$ (MPa)	20.2
Seepage Parameters
$n_e$ (-)	0.389
a (-)	1.2
b (-)	6.2
$K_{int}$ (m${}^2$)	$7\times {10}^{-21}$
$\alpha$ (-)	22
A (-)	2
B (-)	4

, and parameters related to the numerical calculation of the compression and recovery process are shown in

Table 3. Numerical parameters in suction controlled compression.

Paramter	Unit	SP01	SP02	SP03	SP04
M	-	1.07	1.07	1.07	1.07
$P_{0s}$ (-)	MPa	0.67	1.052	1.052	1.052
$P_c$ (-)	MPa	0.01	0.01	0.01	0.01
${\lambda }_0$	-	0.1352	0.104	0.106	0.12
$Kappa$	-	0.0687	0.034	0.0308	0.024
$Kapp{a}_s$	-	0.03	-	-	-
r	-	0.75	0.75	0.75	0.75
$\beta$	-	$0.25\times {10}^{-6}$	$0.25\times {10}^{-6}$	$0.25\times {10}^{-6}$	$0.25\times {10}^{-6}$
$e_o$	-	0.97	0.75	0.73	0.69
$k_s$	-	0.1	0.1	0.1	0.1
${\sigma }_x$	MPa	−1.52	−1.68	−1.443	−1.22
${\sigma }_y$	MPa	−1.52	−1.68	−1.443	−1.22

. In the process of water absorption and expansion, the initial stress value of the model was calculated to reach the fixed suction force. The compression coefficient and recovery coefficient in the process of compression and recovery are obtained by fitting the experimental data, which correspond to the slope of elastic V-ln (P) and the slope of saturated V-ln (P) in Table 3. In order to calculate the void ratio during the compression experiment, the displacement and stress at the node of the symmetric axis on the top boundary of the model during the numerical calculation are output. The model volume change was calculated according to the output displacement, and then the void ratio axial stress curves of the samples SO-01~04 were calculated.

3.3. Result Analysis

Experimental results and numerical simulation results are shown in Figure 6. It can be seen from Figure 6 that the four groups of samples show obvious plastic processes during compression and recovery. The yield stress of samples SO-01~04 increased gradually, indicating that the water absorption process reduced the yield stress of the material. According to the fitted compression parameters during the compression process, the compression parameters of bentonite decreased gradually during the water absorption process. It can be seen from the elastic parameters fitted during the recovery process that the absorption of bentonite reduces the linear elastic modulus of the material.

Figure 7 and Figure 8 show the displacement diagram of samples SO-01 and SO-02 respectively. As can be seen from Figure 7 and Figure 8, the axial displacement of the model in the numerical calculation is small within the compression pressure of 0.5 MPa. When the compression load reaches 10 MPa, the model has significant displacement change and reaches the maximum displacement at the compression load of 50 MPa, and the model is mainly manifested as the deformation of the upper part of the model under compression. Then the displacement of the upper part of the sample decreases gradually with axial pressure decreasing. By comparing Figure 7 and Figure 8, it can be found that the sample with smaller suction has a larger axial displacement, indicating that the water absorption process reduces the stiffness of the bentonite sample.

Figure 9 and Figure 10 reflect the strain cloud diagrams of samples SO-02 and SO-03 under six different compression pressures, respectively. It can be seen from Figure 10 that the strain generated by the model at 0.5 MPa is very small, and when the compression stress reaches 2.5 MPa, obvious strain occurs in the area near the symmetry axis at the upper end of the model. When the compression load reaches 10 MPa, the strain in the upper part of the model near the symmetry axis is large, while the strain in the bottom part of the model is close to 0. In the process of gradually unloading the compression load, the position near the symmetry axis at the bottom of the sample has obvious deformation. When the compressive stress is released to 0.5 MPa, sample SO-03 still maintains a large strain in the center area at the bottom, indicating that the model has a large stiffness under the condition of a large suction force.

4. Results

Based on the implicit integration algorithm of return mapping, the numerical integration expression of the BBM is established, and implement the BBM into the OGS platform for numerical modeling of the coupled HM behavior of unsaturated soil. The numerical basis for studying the elastic–plastic mechanical behavior of bentonite/quartz sand mixture in the coupled process of seepage and stress is provided. The corresponding experimental research is conducted and compared the results with the numerical simulation, the results show that:

(1)

Numerical calculation of void ratio along with the change of axial compression stress results can fit the experimental results well, showing that the numerical calculation described the elastic–plastic mechanics response of bentonite/quartz sand mixture parameters with rationality, and the unsaturated soil elastic–plastic model of the numerical integration algorithm can well reflect the bentonite/elastic–plastic mechanical behavior of the quartz sand mixture.

(2)

In the initial stage of the compression process, the strain distribution of the numerical model is uniform, showing an obvious elastic stage, and then the local strain increases, showing an obvious plastic process.

(3)

Both numerical calculation and experimental results show that the water absorption process can reduce the compression stiffness of bentonite/quartz sand mixture and increase the recovery stiffness.