.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "examples/poromechanics/terzaghi_1d.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_examples_poromechanics_terzaghi_1d.py: Terzaghi 1D consolidation (small strain) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ A column of saturated porous material with a linear elastic skeleton is loaded by a constant pressure on its top face. The top face is drained (``PorePressure = 0``), the bottom and lateral faces are sealed, and the lateral displacements are blocked (oedometric, i.e. 1D, consolidation). When the load is applied the fluid cannot escape instantly: it carries the load and the pore pressure jumps to its undrained value. The fluid then drains through the top face and the load is progressively transferred to the skeleton. The pore pressure at the sealed bottom face is compared with the classical Terzaghi series solution. The model uses Biot poroelasticity: the total Cauchy stress is :math:`\sigma = \sigma_{\mathrm{eff}}(\varepsilon) - \alpha p \mathbf{I}` and the fluid mass balance is .. math:: \frac{1}{M}\dot p + \alpha\, \dot\varepsilon_v - \nabla \cdot \left( \frac{k}{\mu_f} \nabla p \right) = 0. See :doc:`/poromechanics` for details. .. GENERATED FROM PYTHON SOURCE LINES 27-33 .. code-block:: Python import numpy as np import matplotlib.pyplot as plt import fedoo as fd .. GENERATED FROM PYTHON SOURCE LINES 34-40 Material parameters and analytical solution ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ With a single drainage path of length :math:`L`, the closed-form solution uses the oedometric modulus :math:`M_{oed} = K + 4G/3`, the consolidation coefficient :math:`c_v = (k/\mu_f) / (\alpha^2/M_{oed} + 1/M)` and the time factor :math:`T = c_v t / L^2`. .. GENERATED FROM PYTHON SOURCE LINES 40-64 .. code-block:: Python L = 1.0 # column height [m] E, nu = 1.0e6, 0.3 # drained skeleton [Pa] M = 1.0e8 # Biot modulus [Pa] (nearly incompressible fluid) alpha = 1.0 # Biot coefficient k = 1.0e-7 # permeability mu_f = 1.0 # fluid viscosity sigma0 = 1.0e3 # applied surface load [Pa] G = E / (2.0 * (1.0 + nu)) K = E / (3.0 * (1.0 - 2.0 * nu)) M_oed = K + 4.0 * G / 3.0 c_v = (k / mu_f) / (alpha**2 / M_oed + 1.0 / M) p0 = alpha * M * sigma0 / (M_oed + alpha**2 * M) # undrained pore pressure settlement_drained = -sigma0 * L / M_oed def terzaghi_bottom_pressure(T, nterms=300): """Normalized pore pressure p(z=0, t) / p0 at the sealed face.""" m = np.arange(nterms) Mm = (np.pi / 2.0) * (2 * m + 1) return np.sum((2.0 / Mm) * ((-1.0) ** m) * np.exp(-(Mm**2) * T)) .. GENERATED FROM PYTHON SOURCE LINES 65-74 Mesh, constitutive laws and weak form ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The skeleton law is a standard mechanical law. The fluid properties are grouped in :py:class:`~fedoo.constitutivelaw.PoroFluidProperties`. The :py:class:`~fedoo.weakform.PoroMechanics` weak form couples the momentum balance (unknowns ``Disp`` and the skeleton ``Pressure`` multiplier) with the fluid mass balance (unknown ``PorePressure``). :py:class:`~fedoo.weakform.PoroMechanicsSimple` is the two-field alternative, recommended when the model has a free-traction face. .. GENERATED FROM PYTHON SOURCE LINES 74-102 .. code-block:: Python fd.ModelingSpace("3D") nz = 20 mesh = fd.mesh.box_mesh( nx=2, ny=2, nz=nz + 1, x_min=0, x_max=0.1, y_min=0, y_max=0.1, z_min=0, z_max=L, elm_type="hex8", ) skeleton = fd.constitutivelaw.ElasticIsotrop(E, nu) fluid = fd.constitutivelaw.PoroFluidProperties( permeability=k, fluid_viscosity=mu_f, biot_coefficient=alpha, biot_modulus=M, initial_porosity=0.5, ) wf = fd.weakform.PoroMechanics(skeleton, fluid, bulk_modulus=K, nlgeom=False) poro_assembly = fd.Assembly.create(wf, mesh) .. GENERATED FROM PYTHON SOURCE LINES 103-107 Loading and boundary conditions ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The surface load is applied as a step: setting ``initial_pressure`` equal to the final pressure removes the default linear ramp over the time step. .. GENERATED FROM PYTHON SOURCE LINES 107-127 .. code-block:: Python top = mesh.find_nodes("Z", L) bottom = mesh.find_nodes("Z", 0.0) side_x = np.unique( np.concatenate([mesh.find_nodes("X", 0.0), mesh.find_nodes("X", 0.1)]) ) side_y = np.unique( np.concatenate([mesh.find_nodes("Y", 0.0), mesh.find_nodes("Y", 0.1)]) ) load = fd.constraint.Pressure.from_nodes(mesh, top, sigma0, initial_pressure=sigma0) pb = fd.problem.NonLinear(poro_assembly + load) pb.set_nr_criterion("Displacement", tol=1e-4, max_subiter=30) pb.bc.add("Dirichlet", side_x, "DispX", 0.0) pb.bc.add("Dirichlet", side_y, "DispY", 0.0) pb.bc.add("Dirichlet", bottom, "DispZ", 0.0) pb.bc.add("Dirichlet", top, "PorePressure", 0.0) # drained top face .. rst-class:: sphx-glr-script-out .. code-block:: none Dirichlet boundary condition: var = 'PorePressure' n_nodes = 4 value = 0.0 .. GENERATED FROM PYTHON SOURCE LINES 128-132 Time integration ~~~~~~~~~~~~~~~~ The time increments are driven by hand to record the bottom pore pressure, the top settlement and a few pressure profiles along the column. .. GENERATED FROM PYTHON SOURCE LINES 132-161 .. code-block:: Python dt = 0.2 nb_steps = 60 snapshot_steps = [0, 4, 14, 29, 59] pb.initialize() pb.tmax = dt * nb_steps pb.dtime = dt pb.set_start() column = np.where((mesh.nodes[:, 0] == 0.0) & (mesh.nodes[:, 1] == 0.0))[0] column = column[np.argsort(mesh.nodes[column, 2])] p_bottom = np.zeros(nb_steps) uz_top = np.zeros(nb_steps) profiles = {} for step in range(nb_steps): pb.time = step * dt convergence, nb_nr, res = pb.solve_time_increment() if not convergence: raise RuntimeError(f"Step {step}: Newton-Raphson failed (res={res:g})") pb.set_start() p = pb.get_dof_solution("PorePressure") p_bottom[step] = np.mean(p[bottom]) uz_top[step] = np.mean(pb.get_dof_solution("DispZ")[top]) if step in snapshot_steps: profiles[step] = p[column].copy() .. GENERATED FROM PYTHON SOURCE LINES 162-164 Comparison with the Terzaghi solution ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ .. GENERATED FROM PYTHON SOURCE LINES 164-199 .. code-block:: Python t = (np.arange(nb_steps) + 1) * dt T = c_v * t / L**2 T_fine = np.linspace(1e-3, T[-1], 300) series = [terzaghi_bottom_pressure(Ti) for Ti in T_fine] print(f"Undrained pore pressure: fedoo {p_bottom[0]:.4g} Pa, analytical {p0:.4g} Pa") print( f"Final settlement: fedoo {uz_top[-1]:.4g} m, " f"drained analytical {settlement_drained:.4g} m" ) fig, axes = plt.subplots(1, 2, figsize=(11, 4.5)) ax = axes[0] ax.plot(T_fine, series, "k-", label="Terzaghi series") ax.plot(T, p_bottom / p0, "o", ms=4, label="fedoo") ax.set_xlabel(r"time factor $T = c_v t / L^2$") ax.set_ylabel(r"$p(z=0) / p_0$") ax.set_title("Pore pressure at the sealed bottom") ax.legend() ax.grid(True) ax = axes[1] z = mesh.nodes[column, 2] for step, p_z in profiles.items(): ax.plot(p_z / p0, z, label=f"T = {T[step]:.2f}") ax.set_xlabel(r"$p / p_0$") ax.set_ylabel("z [m]") ax.set_title("Pore pressure profiles") ax.legend() ax.grid(True) fig.tight_layout() plt.show() .. image-sg:: /examples/poromechanics/images/sphx_glr_terzaghi_1d_001.png :alt: Pore pressure at the sealed bottom, Pore pressure profiles :srcset: /examples/poromechanics/images/sphx_glr_terzaghi_1d_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Undrained pore pressure: fedoo 982.6 Pa, analytical 986.7 Pa Final settlement: fedoo -0.0007297 m, drained analytical -0.0007429 m .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.866 seconds) .. _sphx_glr_download_examples_poromechanics_terzaghi_1d.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: terzaghi_1d.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: terzaghi_1d.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: terzaghi_1d.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_