Simulation aided testing of hydro-mechanical processes on clay

This paper presents a study focused on the hydro-mechanical behaviour of a plastic clay under partially saturated conditions. Clay remoulded samples were dried using vapour transfer under one-dimensional conditions. Samples underwent an important vertical shrinkage on drying, which progressed along more than one week. To study the time evolution of this phenomenon, simulation aided techniques were used to analyse the progression of suction at local scale and to determine the time required to equilibrate homogeneously the sample. In this case the fully coupled thermo-hydromechanical finite element code Code_Bright [1] was used to analyse this process. The simulation adequately followed the evolution of sample shrinkage measured experimentally. The determination of the suction equalisation time also agreed with the experimental information.


Introduction
Geo-engineering problems such as those related to earthwork constructions (earth dams and embakments) and waste disposal barriers often require a good understanding of the coupled hydromechanical behaviour of geomaterials.Therefore, important developments both in experimental techniques and numerical codes have been undertaken to study coupled phenomena during the last 20 years.Despite the advanced experimental techniques developed to study these processes, their results do not always present the complete picture of understanding, because usually only boundary conditions are applied to macroscopic level.Information on local behaviour and microscopic level usually remains unknown, unless special techniques are used.To partially overcome this problem, additional experimental techniques (tomographic and microstructural) are currently used to better define qualitatively and quantitatively local and microscopic aspects.In some cases, this information at both scales has been included into numerical models through upscaling and homogenisation methods to improve the overall picture of understanding of the physical phenomena.
This paper describes one aspect of this coupled hydro-mechanical phenomena on an unsaturated plastic clay, which is related to the shrinkage undergone on drying [2].A numerical analysis is also performed to complement the hydro-mechanical information of relative humidity application at local level (inside sample).The first part of the paper describes the material used and the test protocols EPJ Web of Conferences followed.The second part of the paper presents the experimental results of the equalisation stage followed by its numerical analysis.Finally, some concluding remarks are addressed.

Material and test protocols
Laboratory tests were conducted on remoulded Boom clay (a deep Tertiary clay formation from Belgium).The liquid limit is 55%, the plastic limit 28%, the specific gravity of solid particles 2.70 and the percentage of clay fraction < 2m is 40% [3].Material passing ASTM #40 was mixed with distilled water to a water content around 34% (initial void ratio around e o =1.07).The relationship between degree of saturation of the material and suction or relative humidity is described through the soil water retention curve.Two procedures were used to determine this curve, namely using highrange tensiometers (CERMES, Paris) up to matric suction of 2 MPa and chilled-mirror dew-point psychrometers (WP4, Decagon devices, Inc, USA) for the higher total suction range.Figure 1 presents the soil water retention curve in terms of degree of saturation (S r ) starting from a remoulded state.The experimental data is fitted by a modified form of van Genuchten's expression [4] proposed by [5].In this expression set by Equation 1, a maximum total suction is considered at a=a r =300MPa for S r =0.Parameter 1/ in Equation 1 is related to the air entry value, n to the slope of the curve and parameter m to the residual degree of saturation.The parameter C(s) is used to correct the curve in the upper suction range.Values of fitted parameters are shown in Figure 1. ( Figure 2 shows a scheme of the adapted Bromhead ring shear apparatus with a glass chamber to allow for sample drying [6].A forced convection system driven by an air pump is used to transfer vapour generated in a vessel containing a saline solution.Vertical displacements of the sample with an initial height of 5 mm, as well as relative humidity and temperature evolutions on drying are automatically registered.A first stage of the test is performed by applying a vertical stress of 100 kPa to the quasi-saturated sample, which develops instantaneous and delayed (due to consolidation) deformations.This first stage is allowed equalising for 24 hours.Afterwards, relative humidity at the top sample boundary is applied through the vapour transfer technique and by using the forced convection circuit.From this moment on, drying process starts leading to a reduction in sample height induced by the loss of water.Sample is usually considered equilibrated on a volume change rate basis of less than 0.2%/day.

Experimental results and numerical interpretation 3.1 Experimental results
Figure 3 shows the first two stages of the test protocol, namely the vertical stress application followed by suction equalisation.Axial strain, as well as relative humidity and temperature at the top sample boundary were measured with time.After the initial vertical stress application, the sample underwent an axial strain of  a =4.92% corresponding to a vertical displacement of u yy =0.244 mm.As shown in Figure 3 during this time relative humidity increased from RH=57% (laboratory conditions) up to RH=95% (corresponding to quasi-saturated sample conditions).After load application, the vessel and the air pump were connected to the glass chamber.Once suction is applied by vapour transfer to the sample boundary, the specimen starts to shrink, undergoing important vertical displacements as time progresses.As observed, the relative humidity decreased sharply to RH=73% after one hour of transferring vapour with the forced convection circuit.Small variations of RH measurements are a consequence of temperature fluctuations.Suction equalisation was assumed at t=192 hours after suction application, when the vertical displacement rate due to shrinkage was less than 0.2%/day.As is shown in Figure 3, maximum axial strain at this time is  a =18.8% for a total suction of s=58MPa estimated from the psychrometric law.However, there is no mean to measure suction inside the sample and know if hydraulic equalisation is correctly achieved.To overcome this problem, a numerical simulation is used to analyse the drying process, which is discussed in the next section.14th International Conference on Experimental Mechanics

Numerical interpretation of the equalisation stage
A series of numerical computations have been performed to simulate the equalisation curve presented in Figure 3 using the fully coupled thermo-hydro-mechanical finite element code CODE_BRIGHT [1].The main challenge of the simulation is to compute numerically the important changes in axial strain and its evolution with time and verify if suction stabilises at the elapsed time used in the experimental procedure.The selected mechanical constitutive model on monotonic drying was the nonlinear surface state model for partially saturated soils [7].This model is selected because of its simplicity to obtain important volumetric changes with a reduced number of parameters.The model requires three parameters (a 1 , a 2, and a 3 ) to define the stress-volume change behaviour.The model also requires the value of the shear modulus G o or the Poisson's ratio  to define the increments of shear strain.Coefficient a 1 is defined by Equation 2 in terms of  which is the slope of the unloading/reloading curve in the e-ln(p') plane (e: void ratio; p': mean net stress).Parameter a 2 is defined by Equation 3 as a function of s  which is the slope of the unloading/reloading curve in the e-ln[(s+p atm )/p atm ] plane (s: suction; p atm : atmospheric pressure).Volumetric strain is calculated by Equation 4. Coupled parameter a 3 was not included in the simulation for simplicity.Values of these coefficients were obtained by back analysis leading to the following values: a 1 =-0.036, a 2 =-0.033 and a constant Poisson's ratio of = 0.35.
    For the hydraulic analysis, the numerical code takes into account changes of water permeability induced by variations of porosity and the degree of saturation.Following this idea, the hydraulic conductivity is defined in the code by Equation 5where, k w is the tensor of intrinsic permeability (affected by porosity changes), k rl is the liquid phase relative permeability (that takes into account degree of saturation) and  l is the liquid dynamic viscosity.In this particular case, a Koseny's model given by Equation 6is used to define the variation of the tensor of intrinsic permeability as function of porosity .In this equation, k wo is the liquid intrinsic permeability at a given reference porosity.A saturated water permeability determination was carried out to define this value.The value of the obtained intrinsic permeability tensor is k wo =2.0E-18m 2 at a reference porosity of  o =0.492.The effect of the degree of saturation is estimated by a power law of the effective degree of saturation S e given by Equation 7with values of n=3, maximum saturation S ls =1.0, and minimum saturation of S rl =1.0.On the other hand, gas permeability is defined in analogous way to the hydraulic conductivity by Equation 8, where k g is the intrinsic permeability tensor for gas, which is related to the liquid permeability tensor through Equation 9.In this case, a value of A=4.0E5 was used in accordance to [3].The parameter k rg is the relative permeability to gas, which is defined by Equation 10.A value of n=2 was adopted from data reported by [3].
The application of relative humidity is controlled by the water mass transfer indicated by Equation 11.In this equation, w g j is the mass-flow of vapour in the gas phase, w g  is the mass of vapour in the gas phase,  g is the gas density and the product of these two quantities is defined as the vapour density w g  given by Equation 12. Evaporation comes from the difference between the vapour density of the soil surface ( w g g   ) and the vapour density in the environment   0 w g g   .The parameter  g is associated with the turbulent exchange function [8] and depends on the movement of the air above the soil surface.Values of this parameter were obtained by [2] using a circulating vapour set-up  g =4.0E-4 m/s.This parameter is very sensitive to the type of circulating system used.
If vapour is forced through the sample with an air pump causing a gradient of gas pressure between the inflow and outflow faces of the sample, advective flow will predominate with high  g values.This means that evaporation takes place rapidly.On the other hand, when the circulation of vapour is applied along the drainage boundaries, evaporation is affected and diffusion takes place at low  g values.In this case, relative humidity is applied along the sample boundaries with a predominant diffusion mechanism. g was obtained by back analysis, giving a value of  g =7.0E-4 m/s, which is in the same order of magnitude to that reported by [2].
The original 3D problem can be simplified into an axi-symmetric simulation due to the symmetry of the annular sample around the vertical axis.For the simulations, the initial conditions of the sample correspond to the state after load application in the experiment.Figure 4 shows the geometry with a sample height of 4.756 mm and 7.5 mm width, together with the boundary conditions and values of the initial state of the sample.In the model, the vertical sides of the sample are considered impervious to water, q w =0, and no horizontal displacements are allowed u xx =0.In the base of the model, a restriction in displacements in the vertical direction is set u yy =0 and at the top of the sample the vertical total stress is applied  v =0.2 MPa.Relative humidity is applied on base and top drainage boundaries trying to follow the experimental measurements.In seek of simplicity a temperature value of T=21.5ºC which corresponds to the mean value of the experimental data is selected for the simulation.
A comparison between the computed and experimental curves is shown in Figure 5.A good agreement is obtained between the two curves.In Figure 6 the evolution of suction at mid-point of the sample is plotted.The figure indicates that suction takes around 168 hours to achieve a stable value.This time is lower than the one determined experimentally (t=192hours).

Concluding remarks
A particular aspect of the coupled hydro-mechanical behaviour of unsaturated plastic clay, related to the shrinkage undergone on drying have been analysed from experimental and numerical standpoints.As shown in the paper, the simulation adequately followed the evolution of sample shrinkage measured experimentally.This simulation allowed determining the equalisation time from a hydraulic viewpoint, when suction at sample mid-height reaches a constant value.This equalisation time was in agreement with the experimental criterion adopted based on volume change rate (<0.2%/day).The results indicate that numerical simulation tools are helpful for the understanding of local processes, allowing for complementary information on the evolution of local (inside specimen) variables and offering additional confidence on the experimental results.

Fig. 2 .
Fig. 2. Scheme of the vapour transfer technique used to dry the sample.

Fig. 3 .
Fig. 3. Evolution of axial strain during stress and suction applications.