Wicking in a partially saturated porous medium

The permeability of the fully saturated material is 3.7 × 10⁻⁴ m/sec. The default model is used for the partially saturated permeability. This assumes that the permeability varies as a cubic function of saturation. The specific weight of the fluid is 10⁴ N/m³. The capillary action in the porous medium is defined by the absorption/exsorption curves shown in Figure 2. These curves give the (negative) pore pressure versus saturation relationship for absorption and exsorption behavior. The transition between absorption and exsorption and vice-versa takes place along a scanning slope which in this example is set by default to 1.05 times the largest slope of any branch in the absorption/exsorption curves. The fluid is assumed to be water, with a bulk modulus of 2 GPa. For the mechanical properties we assume the material is elastic with Young's modulus 50000 Pa and Poisson's ratio 0.0. The dry mass density of the material is 100 kg/m³.

The initial condition for saturation is 5% throughout the column. The initial conditions for pore pressure must have a gradient that is equal to the specific weight of the fluid so that, according to Darcy's law, there is no initial flow. For this purpose we assume the initial pore pressures vary linearly from −12000 Pa at the bottom of the column to −22000 Pa at the top of the column. These initial conditions satisfy the pore pressure/saturation relationship in that they are between the absorption and exsorption curves, as shown in Figure 2. The initial void ratio is 5.0 throughout the column. In the deforming column case the initial conditions for effective stress are calculated from the density of the dry material and fluid, the initial saturation and void ratio, and the initial pore pressures using equilibrium considerations and the effective stress principle. The procedure used is detailed in Geostatic stress state. It is important to specify the correct initial conditions for this type of problem; otherwise, the system may be so far out of equilibrium initially that it may fail to start because converged solutions cannot be found.

The weight is applied by gravity loading. In the case of the deforming column, an initial geostatic step is performed to establish the initial equilibrium state. The initial conditions in the column exactly balance the weight of the fluid and dry material so that no deformation or fluid flow takes place. Then the bottom of the column is exposed to fluid by prescribing zero pore pressure (corresponding to full saturation) at those nodes during a transient soils consolidation step. The fluid will seep up the column until the pore pressure gradient is equal to the weight of the fluid, at which time equilibrium is established.

The transient analysis is performed using automatic time incrementation. The pore pressure tolerance that controls the automatic incrementation is set to a large value since we expect the nonlinearity of the material to restrict the size of the time increments during the transient stages of the analysis and we do not wish to impose any further control on the accuracy of the time integration. The check on displacement and pore pressure changes is relaxed using solution controls. The analysis can also be done with the default tolerances, but Abaqus iterates a lot more without any gain in solution accuracy.

The choice of initial time increment in these transient partially saturated flow problems is important to avoid spurious solution oscillations for some element types (seePartially saturated flow in a porous medium). As discussed in Coupled pore fluid diffusion and stress analysis, the criterion for a minimum usable time increment in partial-saturation conditions is

where $\gamma$ is the specific weight of the wetting liquid, $n^0$ is the initial porosity of the material, k is the fully saturated permeability of the material, $k_s$ is the permeability-saturation relationship, $ds/du$ is the rate of change of saturation with respect to pore pressure as defined in the absorption/exsorption material behavior (Sorption), and $\mathrm{\Delta }\mathrm{\ell }$ is a typical element dimension. For our model we have $\mathrm{\Delta }\mathrm{\ell }=$ 0.1 m (the size of an element side), $\gamma =$ 1.0 × 10⁴ N/m³, $k=$ 3.7 × 10⁻⁴ m/sec, $k_s=s^3$, and $n^0=$ 5/6. Adjacent to where we apply the fully saturated boundary condition, elements will span a region from initial to full saturation early in the transient. A conservative estimate of the minimum time increment is found by choosing the initial saturation of 0.05. From this, we compute $k_s$, $ds/du$, and a value of $\mathrm{\Delta }t$ of about 2700 sec. We find, in practice, that an initial increment of 1000 sec is adequate to avoid oscillations in this problem. For the remaining input files the initial time increment is chosen as 1 second.

In this analysis the prevailing pore pressure in the medium approaches the magnitude of the stiffness of the material skeleton elastic modulus. When reduced-integration elements are used in such cases, the default choice for the hourglass stiffness control, which is based on a scaling of skeleton-material constitutive parameters, may not be adequate to control hourglassing in the presence of the relatively large pore pressure fields. An appropriate hourglass control setting in these cases should scale with the expected magnitude of pore pressure changes over an element and must be defined explicitly by the user.