Partially saturated flow in a porous medium

The permeability of the fully saturated material is 3.7 × 10⁻⁴ m/sec. The partially saturated permeability is the default model, which assumes that the permeability varies as a cubic function of saturation. The specific weight of the material is 10⁵ N/m³. The initial void ratio is 5.0 throughout the sample. 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 that is set by default to 1.05 times the largest slope of any branch in the absorption/exsorption curves. The initial conditions for pore pressure and saturation are assumed to be those at the beginning of the absorption curve, so the initial saturation is 0.05 and the initial pore pressure is −10 kPa.

The gel particles have a radius of 0.5 mm when completely dry and are capable of swelling to a maximum radius of 1.5 mm when fully exposed to fluid. There are 1.0 × 10⁸ gel particles in each cubic meter of the porous material. The relaxation time constant for swelling of the gel particles is 500 sec. A dummy elastic modulus of 10⁴ N/m² is prescribed to complete the material definition.

The “loading” consists of prescribing a zero pore pressure (corresponding to full saturation) at the center of the sample, node 1 in Figure 1. This is based on the assumption that, in the demand wettability test, the sample has available to it as much fluid as necessary to cause saturation at that point. This boundary condition is held fixed for 600 seconds to model the fluid acquisition process. Then a draining period of 600 seconds is modeled by prescribing a pore pressure of −10 kPa at node 1; this corresponds to a saturation of 10%, which is the least saturation the sample can achieve after it has been wetted (see Figure 2). The draining procedure we use is not physically realistic, since the free fluid saturation cannot drop to 10% instantaneously as we assume. Nevertheless, it serves to illustrate the behavior of the model. Finally, we model the rewetting process over a time period of 800 seconds in the third step of analysis by once again prescribing zero pressure at the center of the sample.

The analysis is performed with a transient soils consolidation procedure (Coupled pore fluid diffusion and stress analysis) 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 volume flux tolerance that controls the accuracy of the solution of the flow continuity equations is set to 1.0 × 10⁻⁸ m³/sec (Convergence criteria for nonlinear problems). This is less than 1% of the flux that occurs during the initial wetting stage of the problem. The problem can also be run with the default tolerance. The results obtained are unchanged. However, in this problem the default tolerances calculated by Abaqus are extremely tight, resulting in additional iterations without benefit to the solution. For this reason we define a less stringent tolerance.

An important issue in these transient, partially saturated flow problems is the choice of initial time step. As in any transient problem the spatial element size and the time step are related, to the extent that time steps smaller than a certain size give no useful information. This coupling of the spatial and temporal approximations is always most obvious at the start of diffusion problems, immediately after prescribed changes in the boundary values. As discussed in Coupled pore fluid diffusion and stress analysis, the criterion 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 }=$ 5.08 mm (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. Near node 1, where we apply the boundary condition, we will approach full saturation conditions early in the transient. If we choose a saturation of 0.9 and the corresponding $ds/du$ in the absorption curve, we can calculate a $\mathrm{\Delta }t$ of about 0.1 sec. We choose to start the analysis with a time increment of 1 sec.

If time increments smaller than the critical value are used, spurious oscillations may appear in the solution (except when reduced-integration, linear, or modified triangular elements are used, in which case Abaqus uses a special integration scheme for the wetting liquid storage term to avoid this problem). If the problem requires analysis with smaller time increments than the critical value, a finer mesh is required. Generally there is no upper limit on the time step, except accuracy, since the integration procedure is unconditionally stable.