In transient analysis of unsaturated flow or consolidation analysis with saturated flow, small time increments may lead to spurious, non-physical oscillations in the pore water pressure evolution and lead to convergence difficulties, especially for hydro-mechanically coupled material behaviour. Therefore, a critical time increment depending on the current conditions should not be undercut at any time of the analysis.
Within normal transient analysis, for avoidance of a time increment smaller than the critical one, the initial time increment and the minimum time increment of the respective simulation step should be chosen as \(\Delta t_{0} > \Delta t_{\text{crit}}\) and \(\Delta t_{\min} > \Delta t_{\text{crit}}\) respectively, and the critical time increment must be determined analytically in advance by the user. Criteria for critical time increments exist in both, saturated flow in consolidation analysis and unsaturated flow analysis upon infiltration.
During unsaturated flow, the (minimum) critical time increment is not constant, but inversely proportional to the non-constant diffusion coefficient \(D\) (in m\(^2\)/s), which itself is a function of the porosity \(n\), the specific weight of water \(\gamma^w\), the current hydraulic conductivity as a product of saturated hydraulic conductivity \(k^w\) and relative permeability \(k^{r,wf}\), and the rate of change of saturation with respect to pore water pressure \(\text{d} S^w/ \ \text{d} p^w\):
\[
\Delta t_{\text{crit}} \sim \frac{1}{D} \quad \text{with} \quad D = \frac{k^{r,wf} \, k^w}{\gamma^w \, n \dfrac{\text{d} S^w}{\text{d} p^w}}.
\]
The conservative choice may be to consider the expected maximum value of \(\Delta t_{\text{crit}}\) during the analysis. However, this effectively cancels out convergence-based time stepping in the analysis, whenever the current critical time increment is smaller than the defined maximum. In some cases, this completely negates convergence and the analysis terminates. To account for the variation of \(D\) upon infiltration, which may span over several orders of magnitude, the geometric mean value of the diffusion coefficient can be used to compute a representative value of the critical time increment \(\Delta t_{\text{crit}}\) over the course of infiltration1. It computes as \(D_{\text{mean}}=\sqrt{D_{\min} D_{\max}}\), wherein \(D_{\min}\) and \(D_{\max}\) correspond to the expected minimum and maximum value of the diffusion coefficient, which also must be determined analytically in adavnce following above Equation.
The critical time increment is evaluated within numgeo and is available as element output. It can be:
used to check the choice of \(\Delta t_{\min}\) in transient simulations,
Unsaturated flow: The critical time increment was derived from the mass balance equation for unsaturated flow diffusion problems 23 and seepage analysis 45 and is propsed as:
\(\Delta t_{\text{unsat}}\) is evaluated for unsaturated conditions, i.e. \(\dot S^e > 0\) and \(S^e \leq 0.999\).
Fully saturated conditions: The critical time increment \(\Delta t_{\text{sat}}\) is estimated based on an analytical expression of the one-dimensional consolidation problem of a homogeneous elastic saturated porous medium 6:
\(\Delta t_{\text{sat}}\) is only evaluated for saturated conditions (\(\dot S^e = 0\) and \(S^e \geq 0.999\)) in the presence of solid-phase displacement (\(|| v^s > 1e-4 ||\)).
Therein are:
\(\gamma^w\) the specific weight of water (in kN/m\(^3\))
\(n\) is the porosity
\(\text{d} S^w/ \ \text{d} p^w\) the rate of change of saturation with respect to pore water pressure (in 1/kPa, returned by the hydraulic model )
\(k^w\) the saturated, hydraulic conductivity of the soil (in m/s) in the governing direction, thus flow for \(\Delta t_{\text{unsat}}\) (parallel to the pore water pressure gradient \(p^w_i\)) and displacement for \(\Delta t_{\text{sat}}\) (parallel to the solid velocity \(v^s\)):
\(\Delta l_e\) the representative element dimension (in m), evaluated with a projection of the coordinates \(x_j^n\) of all nodes \(n\) of the finite element with the pore water pressure as a degree of freedom, in the governing direction. For elements with linear interpolation of the pore water pressure, \(\Delta l_e\) corresponds to the element length, whereas in needs to be cut in half for quadratic inteprolation of the pore water pressure.
\(M_{el}\) the constrained modulus of the material, approximated from the current material tangent stiffness \(\frac{\partial \sigma_{ij}}{\partial \varepsilon_{kl}}\) assuming incremental linear elasticity.
Consideration in time-stepping
The critical time increment of the whole FE model \(\Delta t_{\text{crit}}\) in each calculation increment is determined as the critical time increment of the critical element \(\Delta t_{\text{crit}}^{el_\text{crit}}\). The element is identified in the following manner:
In each element, the maximum value of the critical time increments \(\Delta t_{\text{sat}}\) and \(\Delta t_{\text{unsat}}\) in all integration points is determined:
Furthermore, the critical increment of pore water pressure is computed in each element based on the maximum value of \(\Delta t_{\text{unsat}}\) and the corresponding rate of pore water pressure \(\text{d} p^w /\text{d} t\):
The critical element \(el_{\text{crit}}\) is identified in the element assembly as that element with the largest critical increment of pore water pressure \(\Delta p^{w,el}_{\text{crit}}\) for unsaturated conditions or largest critical time increment \(\Delta t_{\text{crit}}^{el}\) for saturated conditions:
The critical time increment of the critical element \(\Delta t_{\text{crit}}^{el_\text{crit}}\) is returned as the critical time increment of the whole FE model \(\Delta t_{\text{crit}}\):
An exemplary schematic evolution of the time increment \(\Delta t\) upon an infiltration is shown in Figure 1. The permissible range for \(\Delta t\) is the overlapping area of the critical range (red) between \(\Delta t_{\mathrm{crit}}\) and \(\Delta t_{\max}\) and the convergence based range (grey). First, \(\Delta t_{\mathrm{crit}}\) increases and thus does \(\Delta t\). When \(\Delta t_{\mathrm{crit}}\) decreases, \(\Delta t\) increases further until \(\Delta t_{\max}\) is reached.
Figure 1: Exemplary evolution of time increment \(\Delta t\) in analysis type "transient-critical".
Output variables
The following element output variables are available:
Keyword
Variable
DT-UNSAT
Average value of \(\Delta t_{\text{unsat}}\)
DT-SAT
Average value of\(\Delta t_{\text{sat}}\)
DELTA-T
Critical time increment of the element \(\Delta t_{\text{crit}}^{el}\)
DELTA-PW
Critical increment of pore water pressure in the element \(\Delta p^{w,el}_{\text{crit}}\)
Furthermore, the critical time increment of the whole FE model \(\Delta t_{\text{crit}}\) and the respective element number \(el_{\text{crit}}\) are provided in the statistics file of the calculation.
H. R. Thomas and Z. Zhou. Minimum time-step size for diffusion problem in fem analysis. Int. J. Numer. Meth. Engng., 40(20):3865–3880, October 1997. ↩
Muthusamy Karthikeyan, T.-S. Tan, and Kok-Kwang Phoon. Numerical oscillation in seepage analysis of unsaturated soils. Canadian Geotechnical Journal, 38:639–651, January 2001. doi:10.1139/t01-018. ↩
