COMSOL 6.4 - Shallow Foundation on Unsaturated Soil

The hydromechanical behavior of unsaturated soil, such as settlement and heave, is an important topic in geotechnical engineering. Settlement and heave largely depends on the saturation and suction in the soil and on the loading conditions. To study these phenomena, this example looks at a shallow foundation resting on an unsaturated soil. The problem setting and geometry is inspired by the example presented in Ref. 1. In order to demonstrate the settlement and heave of the unsaturated soil, the Modified Cam-Clay model (MCC) and Extended Barcelona Basic model (BBMx) are used. Of these two constitutive models, the latter is more suitable to model unsaturated soils since it includes the effect of suction in its constitutive relationship.

Model Definition

In this example, a 10 m wide and 5 m deep soil stratum is studied. A 1 m wide footing is placed on top of the layer, which applies an incrementally increasing footing pressure. Initially, the ground water level is 3 m below surface. This level then defines the phreatic line, under which the soil is saturated and above which it is unsaturated. In order to study the effects of wetting, the ground water level is increased until the soil is fully saturated after the full footing pressure is applied.

Figure 1: Dimensions, boundary conditions, and pressure load for the unsaturated soil.

Soil Properties

The properties of the soil stratum are given in Table 1.

Table 1: Material properties.
Property	Variable	Value
Poisson’s ratio	ν	0.2
Soil density	ρs	1743 kg/m3
Water density	ρw	1000 kg/m3
Dynamic viscosity of water	μw	0.001 Pa.s
Swelling index	κ	0.006
Swelling index for changes in suction	κs	0.008
Compression index at saturation	λ	0.22
Compression index for changes in suction	λs	0.123
Slope of critical state line	M	1.24
Weight parameter	w	0.4
Soil stiffness parameter	m	50 kPa
Plastic potential smoothing parameter	bb	10
Tension to suction ratio	kb	0.6
Initial yield value for suction	sy	100 kPa
Initial void ratio	e0	1.8
Reference pressure	pref	18 kPa
Initial consolidation pressure	pc0	80 kPa

Constraints and Loads

•

The soil stratum is supported by a rigid and perfectly rough base. Apply a fixed constraint on the lower horizontal boundary.

•

Use roller boundary conditions on the left vertical boundary, symmetry to the right.

•

A pressure head of 3 m is assigned in the Richard’s Equation interface. This pressure head is afterward increased to 5 m.

•

The gravity load is applied using a node. The pore pressure in the saturated region of the soil sample is applied using an node.

•

A boundary load is applied on top of the soil layer to model the weight of the footing, see Figure 1.

Results and Discussion

Figure 2 shows the pore pressure at different ground water levels in the soil. Negative values of the pore pressure indicate regions that are unsaturated.

Figure 2: Pore pressure at different ground water levels.

For saturated soils, both the MCC and the BBMx models use the effective stress principle to account for the pore pressure; while for unsaturated soils, the BBMx model accounts for the negative pore pressure by including suction as a model parameter in the yield function and plastic potential. Unlike the BBMx model, suction is not included in the MCC model, which makes it less suitable for modeling unsaturated soils.

The distribution of volumetric plastic strain obtained with both material models, at a ground water level equal to 3 m, is shown in Figure 3. With the MCC model, a larger region of the soil underneath the footing is subjected to plastic deformation, as compared to the results from the BBMx model. The same observation is made in Ref. 1. However, when the ground water level increases, the reduction in suction results in additional plastic strains with the BBMx model. In contrast, as shown in Figure 4, no significant additional plastic strains are observed when the MCC model is used; negligible changes in plastic strain can, however, be observed due to changes in the effective stress. The distribution of von Mises stress for the fully saturated soil is shown in Figure 5 and Figure 6. It is distinctively different when comparing the two soil models. These results emphasize the importance of including suction as a constitutive parameter in soil models intended for modeling partially saturated soils.

Figure 3: Volumetric plastic strain for a ground water level equal to 3 m.

Figure 4: Volumetric plastic strain for a ground water level equal to 5 m.

Figure 5: Distribution of von Mises stress with the MCC model for a ground water level equal to 5 m.

Figure 6: Distribution of von Mises stress with the BBMx model for a ground water level equal to 5 m.

The curve of footing pressure versus settlement is shown in Figure 7. The first stage of the analysis is done at a water level of 3 m, then, the footing pressure is incrementally increased up to 130 kPa. In the elastic regime both models give the same response, but as plastic strains develop, the two models differ significantly. The BBMx model shows smaller settlement as compared to the MCC model; this indicates that if suction is not taken in to account, the settlement is overestimated.

When the ground water level gradually increases to the surface level, the two material models react differently, see Figure 7. The MCC model shows a reduction in the displacement of the footing, while the BBMx model instead shows a large increase in the footing displacement. With the MCC model, the rise in ground water level causes a reduction in effective stress, which gives an overall positive displacement of the soil surface called heave, see Figure 8. Both the footing and the adjacent soil display heaving with the MCC model. With the BBMx model, the footing shows additional settlement due to the reduction in both suction and effective stress. Note that Figure 8 shows the vertical displacement due to wetting only. The behavior shown in Figure 7 and Figure 8 is also reported in Ref. 1. For the BBMx model, the additional footing settlement due to wetting is referred to as soil collapse due to loss of capillary cohesion.

Figure 7: Footing pressure versus settlement.

Figure 8: Vertical displacement of the stratum surface due to wetting.

Notes About the COMSOL Implementation

The model setup neither includes transient phenomena nor the effect that the deformation has on the pore pressure. This means that it would be possible to model this problem without the Richards’ Equation interface by defining the pore pressure in the saturated region using an analytical function. However, the current setup can be easily extended to more complex multiphysics scenarios.

A linear discretization is used for the Solid Mechanics interface to achieve numerical stability for the nonlinear plasticity problem. A dense mesh is used in the domain to maintain good accuracy.

The Cam-Clay family of soil models, like the MCC or BBMx models, do not define any stiffness at zero stress; hence, numerical simulations that use these soil models always prescribe an initial mean stress equal to the reference pressure at zero strain.

Reference

1. A.A. Abed and P.A. Vermeer, “Numerical Simulation of Unsaturated Soil Behavior,” International Journal of Computer Applications in Technology, vol. 34, no. 1, 2009.

Geomechanics_Module/Soil/settlement_analysis

Modeling Instructions

From the menu, choose .

New

In the window, click

Model Wizard

In the window, click

In the tree, select > > .

Right-click and choose .

In the tree, select > .

Right-click and choose .

Click

In the tree, select > .

Click

Geometry 1

Model parameters are available in the appended text file.

Global Definitions

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 settlement_analysis_parameters.txt.

Footing Pressure

In the toolbar, click

and choose > .

In the window for , type Footing Pressure in the text field.

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

In the table, enter the following settings:


t	f(t)
0	0
1	130
1.1	130

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


Argument	Unit
t	1

In the table, enter the following settings:


Function	Unit
F_P	kPa

Ground Water Level

In the toolbar, click

and choose > .

In the window for , type Ground Water Level in the text field.

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

In the table, enter the following settings:


t	f(t)
0	3
1	3
1.1	5

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


Argument	Unit
t	1

In the table, enter the following settings:


Function	Unit
GWL	m

Initial Suction Profile

In the toolbar, click

and choose > .

In the window for , type Initial Suction Profile in the text field.

In the text field, type InitSuction.

Locate the section. In the text field, type rhow*g_const*(Y-3)*(Y>=3).

In the text field, type Y.

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


Argument	Unit
Y	m

In the text field, type Pa.

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


Plot	Argument	Lower limit	Upper limit	Fixed value	Unit
√	Y	0	5	0	m

Definitions

Variables 1

In the window, expand the > node.

Right-click and choose .

Model variables are available in the appended text file.

In the window for , locate the section.

Click

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

Create half of the geometry by exploiting symmetry.

Geometry 1

Rectangle 1 (r1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 5.

Add a line segment to represent the foundation.

Line Segment 1 (ls1)

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the text field, type 4.5.

In the text field, type 5.

Locate the section. From the list, choose .

In the text field, type 5.

Click

Add points on side boundaries to represent the initial ground water level.

Point 1 (pt1)

In the toolbar, click

In the window for , locate the section.

In the text field, type 3.

Materials

Add a that contains information about the fluid and porous matrix properties together with the structural properties.

Porous Material 1 (pmat1)

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

Continue with setting up the physics. After that, the software automatically detects which material properties are required. First, change the discretization to .

Richards’ Equation (dl)

In the window for , click to expand the section.

From the list, choose .

Unsaturated Porous Medium 1

In the window, under > click .

In the window for , locate the section.

From the list, choose .

Porous Matrix 1

In the window, click .

In the window for , locate the section.

From the list, choose . In the Se text field, type Se.

In the θl text field, type S_res*phi0+Se*(phi0-S_res*phi0).

In the Cm text field, type Cm.

In the κr text field, type k_rel.

In the θr text field, type S_res*phi0.

Pressure Head 1

In the toolbar, click

and choose .

Select Boundaries 1 and 3 only.

In the window for , locate the section.

In the Hp0 text field, type GWL(para)-Y.

Symmetry 1

In the toolbar, click

and choose .

Select Boundary 6 only.

Continue with setting up the interface. Change the discretization to .

Solid Mechanics (solid)

In the window, under click .

In the window for , click to expand the section.

From the list, choose .

Modified Cam-Clay Model (MCC)

In the toolbar, click

and choose .

In the window for , type Modified Cam-Clay Model (MCC) in the text field.

Select Domain 1 only.

Locate the section. From the Γ(θ) list, choose .

Find the subsection. From the M list, choose .

In the pref text field, type pref.

In the pc0 text field, type pc0.

External Stress 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the pA text field, type PorePressure.

In the pref text field, type 0.

From the αB list, choose . In the associated text field, type 1.

Go to the material node and assign the required material properties.

Materials

Porous Material 1 (pmat1)

In the window, under > click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Density	rho	rhos	kg/m³	Basic
Poisson’s ratio	nu	Nu	1	Critical state model
Swelling index	kappaSwelling	kappa	1	Modified Cam-clay
Initial void ratio	evoid0	e0	1	Critical state model
Permeability	kappa_iso ; kappaii = kappa_iso, kappaij = 0	k	m²	Basic
Slope of critical state line	M	MB	1	Critical state model
Compression index	lambdaComp	lambda	1	Modified Cam-clay

Solid Mechanics (solid)

Extended Barcelona Basic Model (BBMx)

In the toolbar, click

and choose .

In the window for , type Extended Barcelona Basic Model (BBMx) in the text field.

Locate the section. From the list, choose .

Select Domain 1 only.

From the Γ(θ) list, choose .

Find the subsection. From the M list, choose .

In the s0 text field, type InitSuction(Y).

In the s text field, type Suction.

In the pref text field, type pref.

In the pc0 text field, type pc0.

External Stress 1

In the toolbar, click

and choose .

In the window for , locate the section.

From the list, choose .

In the pA text field, type PorePressure.

In the pref text field, type 0.

From the αB list, choose . In the associated text field, type 1.

Go to the material node and assign the required material properties.

Materials

Porous Material 1 (pmat1)

In the window, under > click .

In the window for , locate the section.

In the table, enter the following settings:


Property	Variable	Value	Unit	Property group
Plastic potential smoothing parameter	bb	bB	1	Barcelona Basic
Swelling index for changes in suction	kappaSwellingsu	kappa_s	1	Barcelona Basic
Weight parameter	wb	wB	1	Barcelona Basic
Initial yield value for suction	sy0	syB0	Pa	Barcelona Basic
Soil stiffness parameter	mb	mB	Pa	Barcelona Basic
Initial void ratio	evoid0b	e0	1	Barcelona Basic
Tension to suction ratio	kb	kB	1	Barcelona Basic
Compression index for changes in suction	lambdaCompsu	lambda_s	1	Barcelona Basic