1D Infiltration (Plaxis)
Introduction
Validation of element formulations (and implementations) for simulations of partially saturated problems is difficult due to lack of analytical solutions. For this reason, we validate the implementations in numgeo based on a comparison with calculation results obtained with the widely used FE program Plaxis. The boundary value problem (BVP) considered for this purpose is taken from the internal report of Plaxis (Vahid Galavi) "Groundwater flow, fully coupled flow deformation and undrained analyses in PLAXIS 2D and 3D".
The BVP considers a recharge situation of a soil column, which is filled from the bottom in opposite direction of the gravitational force. The soil column has a heigt of 2 m and is displayed in below figure. The displacements are constrained at all nodes (only the flow of water is investigated in this example). The initial pore water pressure is -10 kPa in the entire column leading to an initial degree of saturation of approximately 20 %.
Both the soil-water-retention curve and the dependence of the relative permeability on the effective degree of saturation are modelled using the well known van Genuchten model. The bulk modulus of pore water \(K^w\) and the hydraulic conductivity \(K\) are chosen such as described in the internal Plaxis report: \(K^w = 4.875\cdot 10^5\) and \(K= 1.7604 \cdot 10^{-6}\) m/s. The parameters for the van Genuchten model are \(n^{vG} = 2.286\) and \(\alpha^{vG}=0.224\). The initial void ratio is \(e_0 = 0.5625\)
Left: finite element model of the BVP, Middle: initial distribution of pore water pressure, Right: comparison of soil-water-retention curve and relative permeabiltiy used in Plaxis and numgeo.
Input files
The input files for the benchmark simulation can be downloaded here
Material
For the solid a linear elastic constitutive model is chosen. As no soil deformation is considered in this simulation (an neither was observed in the experiment) this choice is completely arbitrary. The Young's modulus is \(50 \cdot 10^3\) kPa and the Poisson's ratio \(0.3\).
Note that numgeo requires the prescription of the permeability \(K^s\) of the solid and the dynamic viscosity of the pore fluids \(\mu^f\) instead of the hydraulic conductivity \(K^f\), which are related as follows: $$ K^f = \dfrac{K^s \gammaf}{\muf} $$
Therein, \(\gamma^f\) and \(\mu^f\) are the specific weight and the dynamic viscosity of the fluid \(f\), respectively.
The dynamic viscosity of pore water is \(\mu^w=10^{-6}\) kPa\(\cdot\)s and of the pore air \(\mu^a=10^{-8}\) kPa\(\cdot\)s. Assuming a specific weight of 10 kN/m\(^3\) for the pore water, the permeability of the soil is calculated to \(1.7604 \cdot 10^{-13}\) m\(^2\).
Geometry and boundary conditions
The column is discretised with 8-noded rectangular elements (quadratic interpolation). The nodal distance is approximately 0.05 m. For this simulation, changes in pore air pressure are judged as negligible, thus elements based on reduced set of governing equations are used - namely the up-formulation. These elements consider negative pore water pressures as suction \(s=-p^w\) (instead of \(s=p^a-p^w\)).
Initial conditions
The initial pore water pressures is \(p^w_0=-10\) kPa and constant over the height of the column. The initial void ratio is \(e_0 = 0.5625\).
Calculation stages
The simulation is divided into 2 steps in total:
- During the Geostatic step, the self weight of the soil (grains and pore water) is applied without generating any deformation. As stated previously, no deformation of the soil skeleton is expected. We therefore constrain the displacements of all nodes in \(x1\) and \(x2\) direction. In addition, we use boundary conditions to prescribe the pore water pressure for each node. As in the initial conditions, the pore water pressure is prescribed as -10 kPa over the entire column.
- During the transient step we simulate the water supply at the bottom of the soil column. This is done by prescribing the pore water pressure at the bottom nodes. The value is set to 15.0 kPa, which corresponds to a pressure head of 1.5 m. The total step time is 323136 seconds. Due to the strong nonlinearities resulting from the saturation-suction relation and the relative permeability function, we limit the maximum allowed time increment size to 2000 seconds. No ramp is used to increase the pore water pressure from its initial value of -10.0 kPa (*Boundary, instant).
Results
The following figure presents the distribution of the degree of saturation over the column height. The comparison of the simulation results obtained with numgeo and the results presented in the internal Plaxis report show a good agreement.
Distribution of degree of saturation along the column height for different time steps of the simulation. numgeo vs. Plaxis.