COMSOL 6.4 - The Mandel

Cryer’s problem is a three-dimensional consolidation benchmark where a porous sphere is subjected to a uniform boundary pressure. The pore pressure at the center of the sphere rises due to the Mandel–Cryer effect, which is captured by a bidirectional coupling between Darcy’s Law and Solid Mechanics interfaces through use of the Poroelasticity multiphysics interface.

In its initial state, the sphere is not deformed and only subject to ambient atmospheric pressure. The sudden load increase on the surface is instantly affecting the pore pressure, while the solid needs time to deform. The pore fluid is squeezed out of the outer layers, first, so they act like a tightening belt. This leads to a pressure increase in the center, where it even exceeds the value of the applied load. Due to drainage the pressure then gradually decreases toward a new equilibrium state.

Model Definition

In the case of small deformations, there is an analytical solution for Cryer’s problem that describes the pressure p (normalized with the pressure load p0) at the center of the sphere as a function of the bulk modulus K, shear modulus G, Biot’s coefficient αB, permeability κ, fluid viscosity μf, and storage porosity Sp:

(1)

Here zi are the roots of the nonlinear equation

(2)

Furthermore, a is the radius of the sphere, and the parameters m and β and the coefficient of consolidation cv are defined as

(3)

The model is set up using the multiphysics interface, where the interfaces and are coupled using the biphasic poroelasticity model. Thus, the solid and fluid constituents of the porous medium are assumed to be incompressible and the storage coefficient Sp is zero. The coupled equation system for solid mechanics and Darcy’s law is

(4)

The deformation gradient F is in general defined as the gradient of the current coordinates with respect to the original coordinates, the Jacobian J is the ratio between the current volume and the initial volume, S is the second Piola–Kirchhoff stress tensor, and fV is a volume force. Moreover, pA is the absolute pressure, which equals the sum of pore pressure and reference pressure pA = pf + pref, and Qm is the mass source term in the Darcy’s Law equation and vd is the Darcy velocity.

A hyperelastic material is used, which is based on a Neo-Hookean material model.

The original geometry is a sphere but due to symmetry it is sufficient to model a 2D axisymmetric component. A symmetry boundary condition is applied at the bottom edge and a constant boundary load is applied on the surface (see Figure 1).

Figure 1: Simplified geometry and prescribed load of the Mandel–Cryer experiment.

The simulation experiment is performed for two different permeability models: for the Kozeny–Carman permeability model and for constant permeability.

Results and Discussion

Figure 2 illustrates the setup: the full spherical geometry and the boundary load, as well as the resulting pore pressure normalized with the pressure load p0 at the surface of the sphere. The pore pressure increases toward the center of the sphere and reaches a maximum at about 21 s of simulated time which exceeds the outer pressure load by up to 50%.

Figure 2: Normalized pore pressure (surface) and prescribed load (arrows) for a pressure load of 1000 Pa after 21 s in the case where the Kozeny–Carman permeability model was used.

Figure 3 and Figure 4 show the pressure, normalized with the external pressure load, as a function of time for the Kozeny–Carman permeability model (Figure 3) and for constant permeability (Figure 4). In both cases, the Mandel–Cryer effect is correctly represented by the model and the model matches the analytical solution for small deformations.

Figure 3: Normalized pore pressure over time at the sphere center for different prescribed loads in the case of Kozeny–Carman permeability.

Figure 4: Normalized pore pressure over time at the sphere center for different prescribed loads in the case of constant permeability.

The volume ratio — that is, the ratio of the current pore volume to the initial pore volume — at the center of the sphere is displayed in Figure 5 and Figure 6; for the first case, which matches the analytical solution, the volume hardly changes, whereas for higher values of boundary pressure large deformation occurs during consolidation.

Figure 5: Volume ratio over time at the sphere center for different loads and Kozeny–Carman permeability values.

The two different permeability models seem to differ only in the cases of higher pressure loads where the volume ratio and the normalized pressure decrease more rapidly with time when a constant permeability model is used.

Figure 6: Volume ratio over time at the sphere center for different loads and constant permeability.

References

1. J. Choo, “Large deformation poromechanics with local mass conservation: An enriched Galerkin finite element framework,” Int. J. Numer. Methods Eng., vol. 116, no. 1, pp. 66–90, 2018; DOI doi.org/10.1002/nme.5915.

2. A. Verruijt, “Theory of Consolidation,” An Introduction to Soil Dynamics, Theory and Applications of Transport in Porous Media, vol. 24, pp. 65–90, Springer, Dordrecht, 2009; doi.org/10.1007/978-90-481-3441-0_4.

Porous_Media_Flow_Module/Verification_Examples/mandel_cryer_effect

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Click .

Click

In the tree, select > .

Click

Global Definitions

Start by defining the model parameters. You can load them from an external file:

Parameters 1

In the window, under click .

In the window for , locate the section.

Click

Browse to the model’s Application Libraries folder and double-click the file mandel_cryer_effect_parameters.txt.

In the table, enter the following settings:


Name	Expression	Value	Description
mv	1/(K0+4/3*G0)	3.3333E-7 1/Pa	Vertical compressibility of the soil in the case of lateral confinement
kappa0	cvmufmv	1E-12 m²	Initial permeability
dp	sqrt(180kappa0((1-poro0)^2/poro0^3))	1.8974E-5 m	Particle diameter

Geometry 1

Circle 1 (c1)

In the toolbar, click

In the window for , locate the section.

In the text field, type R0.

In the text field, type 90.

Click

Click the

button in the toolbar.

Materials

In the next step, define the porous material as empty material node; as soon as the physics has been defined the material node’s context menu will show you which properties are needed for the simulation. You can then just fill in the values defined in the list.

Porous Material 1 (pmat1)

In the window, under right-click and choose > .

Solid Mechanics (solid)

Now set up the physics. Start with the node.

Hyperelastic Material 1

In the toolbar, click

and choose .

Select Domain 1 only.

Symmetry Plane 1

In the toolbar, click

and choose .

Select Boundary 2 only.

Boundary Load 1

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

From the list, choose .

In the p text field, type p0.

Darcy’s Law (dl)

In the window, under click .

In the window for , click to expand the section.

From the list, choose .

Change the shape function to linear to reduce the computational time and to increase numerical stability, especially for higher pressure load values.

Click to collapse the section.

Symmetry 1

In the toolbar, click

and choose .

Select Boundary 2 only.

Outlet 1

In the toolbar, click

and choose .

Select Boundary 3 only.

In the window for , locate the section.

From the list, choose .

Porous Medium: Kozeny-Carman

In the toolbar, click

and choose .

Select Domain 1 only.

In the window for , type Porous Medium: Kozeny-Carman in the text field.

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the list, choose .

In the dp text field, type dp.

Multiphysics

Poroelasticity 1 (poro1)

In the window, expand the > > node, then click .

In the window for , locate the section.

From the list, choose .

From the pref list, choose .

Materials

Having defined the physics, now fill in the empty expressions in the node.

Porous Material 1 (pmat1)

In the window, expand the > > node, then click .

In the window for , locate the section.

Click

Fluid 1 (pmat1.fluid1)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rhof	kg/m³	Basic
Dynamic viscosity	mu	muf	Pa·s	Basic

Solid 1 (pmat1.solid1)

In the window, right-click and choose .

In the window for , locate the section.

In the θs text field, type 1-poro0.

Porous Material 1 (pmat1)

In the window, click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Lamé parameter μ	muLame	G0	N/m²	Lamé parameters
Lamé parameter λ	lambLame	Lambda0	N/m²	Lamé parameters
Density	rho	rhos	kg/m³	Basic

Mesh 1

Free Quad 1

In the toolbar, click

Size

In the window, click .

In the window for , locate the section.

From the list, choose .

Click

Study 1

Step 1: Time Dependent

In the window, under click .

In the window for , locate the section.

In the text field, type 10^{range(-4,4/50,0)}* R0^2/cv. This input corresponds to the expected time dependency as shown in Equation 1, enforcing small time intervals in the beginning which are then exponentially growing. The end time of the study time interval is calculated as the square of the sphere radius divided by the coefficient of consolidation.

In the window, click .

In the window for , locate the section.

Clear the checkbox.

In this case, three different pressure loads are investigated. Therefore, add a parametric sweep.

Parametric Sweep

In the toolbar, click

In the window for , locate the section.

Click

In the table, enter the following settings:


Parameter name	Parameter value list	Parameter unit
p0k (Pressure over bulk modulus)	0.0001 0.25 0.5

In the window, click .

In the window for , type Study with Kozeny-Carman Permeability in the text field.

In the toolbar, click

Global Definitions

To compare the simulation results with the analytical solution as in Figure 3, you can import the analytical data from an external file provided in the applications database and create an interpolation function.

Analytical Solution

In the toolbar, click

and choose > .

In the window for , locate the section.

From the list, choose .

Click

Browse to the model’s Application Libraries folder and double-click the file mandel_cryer_effect_p_analytical.txt.

In the text field, type Analytical Solution.

Locate the section. Click

In the text field, type p_analytical.

Locate the section. In the table, enter the following settings:


Function	Unit
p_analytical	1

In the table, enter the following settings:


Argument	Unit
t	s

Results

In the window, expand the node.

Grid 1D: Analytical Solution

In the window, expand the > node.

Right-click > and choose > .

In the window for , type Grid 1D: Analytical Solution in the text field.

Locate the section. From the list, choose .

From the list, choose .

Locate the section. In the text field, type t.

In the text field, type -33.26333.

In the text field, type 366.29667.

Click to expand the section. In the text field, type 10000.

From the list, choose .

Normalized Pressure (Kozeny-Carman Permeability)

In the toolbar, click

In the window for , locate the section.

From the list, choose .

Click to expand the section. From the list, choose .

In the text area, type Normalized Pressure at center using Kozeny-Carman permeability.

Locate the section. Select the checkbox.

In the text field, type Normalized Pressure (Kozeny-Carman Permeability).

Locate the section.

Select the checkbox. In the associated text field, type t \cdot \[\mathrm{c_{v} / R_{0}^{2}}\], which is LaTex-syntax to write mathematical formulas and expressions.

Select the checkbox. In the associated text field, type \[\mathrm{p / p_{0}}\].

Locate the section. From the list, choose .

Point Graph 1

Right-click and choose .

Select Point 1 only.

In the window for , locate the section.

In the text field, type p/p0.

Locate the section. From the list, choose .

In the text field, type t/(R0^2/cv).

Click to expand the section. From the list, choose .

Find the subsection. From the list, choose .

From the list, choose .