E3S Web of Conferences 415, 01006 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341501006

Constitutive modelling of phase transitions in granular soils
Claudio di Prisco1* and Pietro Marveggio1
1Politecnico di Milano, Dipartimento di Ingegneria Civile e Ambientale, 20133 Piazza Leonardo da Vinci, 32, Milano, Italy

Abstract. The use of numerical codes for simulating the propagation phase of fast landslides requires thedefinition of constitutive relationships suitable for the description of the mechanical behaviour of granularsoils, under both solid and fluid conditions. In this paper, the authors briefly illustrate a constitutive
modelling approach interpreting the material as a bi-phasic continuum, according to which the solid phase
may behave like a solid, under quasi-static conditions, and a fluid, under collisional/dynamic conditions.
The model presented seems to satisfactorily reproduce a series of results obtained by using a discrete element
numerical code on mono and polydisperse granular specimens. The constitutive model has been
implemented in a Material Point Method code to solve boundary value problems, in which strain rate large
values are expected to influence the material response.

1 Introduction
The mechanical behaviour of soils is very fascinating fornumerous reasons: for the apparent simplicity of theirmicrostructure, for the continuous evolution of theirfabric when subject, under quasi-static conditions, tostress perturbations (self-adaptation) and, in particular,
for their capability of behaving like either solids orviscous fluids, according to the current void ratio andthe effective confining pressure imposed [1].
Soils behave like solids when permanent force chains at
the microstructural level develop among particles: inthis case, the soil skeleton gives arise and themechanical energy is stored only elastically by thematerial. The system elastic stiffness is governed, on
one side, by the material constituting grains and, on theother, by the microstructural geometry [2]. In contrast,
granular soils behave like fluids when, because of grain
agitation and/or large solid dilution, force chains cannot
develop, grains collide among each other and energy ismainly stored as fluctuating kinetic energy [3].
An additional peculiarity of granular soils derives fromthe presence, within the representative elementaryvolume (REV) of interconnected voids that, undersaturated conditions, are totally filled with water. For
this reason, in the geotechnical community, focusing on
the analysis of problems in which granular soils behave
like solids, multi-phasic continuum approaches,
according to which in the REV two or three phases(solid, fluid and, under partially saturated conditions,gas) are superimposed. In contrast, in fluid mechanics
community, suspensions are commonly interpreted as aone phase material, whose rheology is governed by thesolid concentration, which is assumed not to evolve
during the processes analysed [4].
This second approach is characterized by a fundamental
limitation: the evolutions/phase transitions from solid to

fluid regimes and vice versa cannot be described. In
other words, this approach cannot be employed todescribe fluidization and solidification processes.

2 Constitutive modelling
From these premises starts the work developed at
Politecnico of Milan in the last decade. This is aimed atdefining a class of constitutive models capable ofsimulating/interpreting the mechanical behaviour ofgranular soils under both fluidized and solid conditions,
to be implemented in numerical codes dealing with large
displacements, with the ambition of numerically
reproducing
inception,
propagation
anddeposition/arrest of fast/catastrophic landslides.
The proposed approach is inspired by the analysis
and the discussion of the results obtained by performinga large number of discrete element numerical tests on
both triaxial and simple shear cubic specimens [5-8].
Both ideal monodisperse and polydisperse granularassemblies were tested, changing loading paths, initialconditions and strain rates. These test results haveallowed to describe the mechanical behavior of thematerial in both quasi-static and dynamic/collisional
regimes, under both evolving and steady conditions andto better understand the mechanisms governing thetransition from one regime to the other.
2.1 Dry conditions
The interpretation of the results and the availabilityof many microstructural information enabled the authors
to conceive a model capable of describing the material
under steady conditions in both quasi-static andcollisional regimes under dry simple shear conditions* Corresponding author: claudio.diprisco@polimi.it© The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0(https://creativecommons.org/licenses/by/4.0/).

E3S Web of Conferences 415, 01006 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341501006[9,10]. The model, accounting for the double nature of
the microstructural dissipative mechanisms (collisionsand permanent force chains evolution), is based on thedefinition of three state variables: the void ratio, thegranular temperature and the fabric, describing this last
one the directional characteristics of the system.
According to this model, obeying to a parallel scheme,
the critical state is only the limit, for granulartemperature tending to zero, of the steady state: the total
stress applied to the REV is therefore the addition of twocontributions, the quasi static and the collisional. The
former one is assumed to obey to the well-establishedrelations of critical state theories for granular materials[11, 12], whereas the latter one to the kinetic theories of
granular gasses [13].
The model is capable of correctly capturing thedependence of stress level 𝜇𝜇 = 𝜏𝜏/𝜎𝜎 (where 𝜏𝜏 and 𝜎𝜎
stand for shear and normal stresses, respectively) on
inertial number (𝐼𝐼 = 𝛾𝛾̇ 𝑑𝑑 /�𝜎𝜎/𝜌𝜌, where 𝛾𝛾̇ stands for shear
strain rate, 𝑑𝑑 for the mean value of grain diameter and 𝜌𝜌
for grain density) and also of void ratio 𝑒𝑒 on 𝐼𝐼, in casegrains are not assumed to be rigid. This means that thisapproach is capable of overpassing a limitation ofstandard 𝜇𝜇 - 𝐼𝐼 rheologies [14] that, in case of deformablegrains, are not capable of correctly reproducing thedependence of void ratio on I when both constantvolume and constant normal stress tests are performed[15].
The model allows to classify three regimes: thequasi-static, where the collisional contribution isnegligible and the energy is mainly stored elastically bypermanent force chains, the collisional, where the quasistatic contribution is nil and the energy is mainly stored
as kinetic fluctuating energy, and an intermediate, wherethe twos are comparable, even if the material mechanicalbehavior seems not to be rate-dependent.

As was previously suggested, state variables govern
the occurrence of one regime with respect to the other.
In particular, for sufficiently large void ratios (𝑒𝑒 > 𝑒𝑒𝑐𝑐 ),force chains cannot develop and the collisional regime
prevails even for small values of 𝑇𝑇. For 𝑒𝑒 < 𝑒𝑒𝑐𝑐 , both twocontributions are in general non negligible and theregime is called intermediate (Figure 1), but forincreasing values of 𝑇𝑇, the dynamic/collisional
contribution prevails again (inertial regime).
The same model was then extended to threedimensional states of stress and this required a careful
discussion about the shape of the critical state locus in
the deviatoric plane and a redefinition of kinetic-theory
functions. Again, the numerical DEM results allowed toemphasize the role of stress anisotropy (term very
common in the hydraulic and physics communities) inaffecting the 𝜇𝜇-I and the e-I relations [15].
One goal of the research is to study and simulate
phase-transition, that is the passage from one regime toanother: a phenomenon taking place when themicrostructure evolves. The numerical tests performed
to describe evolving conditions are of two differenttypes: heating and cooling tests. The former type ischaracterized by a continuous increase in the energy rate
provided to the specimen until a steady value, the latter
one, in contrast, by a continuous decrease. In the case ofheating tests, the initial conditions may be static orsteady, whereas, in the case of cooling, initial conditions
are only steady.

Fig. 2. Stored energy ratio (elastic over kinetic fluctuating)and mechanical coordination number Zc. Dotted line
represents the steady limit. From [12]. Data from [6].
Fig. 1. Phase diagram: quasi-static, collisional andintermediate regimes.

In Figure 2, numerical DEM results of a certain
number of both heating and cooling tests are illustrated
in the 𝑍𝑍𝑐𝑐 –𝐸𝐸𝑒𝑒𝑒𝑒 /𝐸𝐸𝑐𝑐 plane, where 𝑍𝑍𝑐𝑐 stands for the2E3S Web of Conferences 415, 01006 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341501006

mechanical coordination number [16], whereas 𝐸𝐸𝑐𝑐 and𝐸𝐸𝑒𝑒𝑒𝑒 , for the fluctuating kinetic and elastic energy,
respectively. In the same plane, even the steady state line
is plotted. According to this graphical representation,
the solid mechanical response is placed on the right,
whereas the fluid phase on the left. As is evident,
according to the initial conditions imposed and to thechosen type of test, phase transition may take place (forinstance test “I2-H5”) or not (for instance test “I1-H1”).
The intersection point of the two phases corresponds tothe inertial regime [15].
The constitutive model has been thus conceived byconsidering steady state condition as a mathematicalattractor. As far as the quasi-static contribution isconcerned, this implies the employment of the criticalstate theory and its interpretation within the frameworkof anisotropic strain hardening elastic-plasticity bymeans of a suitable definition of the dilatancy rule
[17,18].
When constant pressure 𝑝𝑝 heating tests are simulated
and large strain rates imposed, the material dilationinduces a progressive shrinkage of the yield locus until
its nullification. When this nullification takes place,phase transition occurs and the material fluidizes. Froma thermodynamic point of view, this implies thenullification of the elastic energy stored by permanent
force chains and, consequently, the disappearance ofsoil skeleton.
In contrast, when cooling p-constant tests are
simulated starting from initial steady conditions
characterized by very large strain rates, the continuous
reduction in strain rate may cause a re-solidification of
the material, that is a progressive increase in the yield
locus size.
What was also numerically analyzed is the evolutionof fabric under quasi static, intermediate and collisionalconditions. DEM numerical results clearly show that thesecond invariant of fabric increases in the collisional
regime: collisions are not randomly distributed but arealigned along preferential directions, justifying theincrease in 𝜇𝜇 for increasing 𝐼𝐼 values. Such anobservation seems to suggest that more advanced kinetic
theories should account for fabric tensor (and inparticular its deviator) to better reproduce themechanical behavior of granular gases.
When constant volume heating tests are simulated, if
the initial material void ratio is sufficiently small, themodel is capable of reproducing the materialliquefaction and the transformation of the material into
a “viscous fluid”, characterized by a dependence of themechanical response on the strain rate imposed.

presence modifies the granular collisions, causing theonset of lubrication forces [17].
From a macroscopic viewpoint, water does not affectthe constitutive modelling of the solid/quasi staticcontribution (under quasi static conditions the effective
stress principle continues to be valid and the response of
the two phases to be uncoupled) whereas underdynamic/collisional conditions the hydro-mechanicalcoupling dominates the system response, in particularfor sufficiently small strain rates (Newtonian regime). In
Figure 3 simple shear model predictions under steadystate conditions for different imposed 𝑒𝑒 values areshown. It is evident for very small 𝛾𝛾̇ ∗ = 𝛾𝛾̇ 𝑑𝑑 /�𝐸𝐸𝑝𝑝 /𝜌𝜌(being 𝐸𝐸𝑝𝑝 the particle Young modulus) values thetransition from fluid to solid regimes, testified by anullification of the dependance of 𝜏𝜏 ∗ = 𝜏𝜏/𝐸𝐸𝑝𝑝 on 𝛾𝛾̇ ∗ . In
case of sufficiently large 𝑒𝑒 values, the model is also
capable of simulating the passage from a Newtonian toa Bagnoldian regime, when 𝛾𝛾̇ ∗ is increased. The dashedline in Figure 3 corresponds to the inertial regime, whichis reached even by dense mixtures in case 𝛾𝛾̇ ∗ issufficiently large.

Fig. 3. Model prediction of saturated steady simple shear
conditions: constant volume lines. Data from [17].

3 Numerical implementation
More recently, the model was implemented in a
Material Point Model code to solve boundary valueproblems. An interesting application was thatconcerning the impact of granular masses on rigid
obstacles, during which very large strain rates occur andthe material experiences fast phase transitions(solidification and fluidization). The comparison withDEM numerical results is very instructive, allowing thecomprehension of the dynamic physical phenomena
governing the system response [19,20].

2.2 Saturated conditions
When the granular material is saturated, the hydromechanical coupling gives rise to additional interaction
mechanisms that make more complex the constitutive
modelling. On one side, the presence of grains, behaving
like obstacles for the water flux, cause an increase in theequivalent water viscosity, that becomes, therefore, afunction of solid concentration. On the other side, thewater viscosity reduces the granular agitation, damping
the fluctuating energy of grains. Moreover, the water

The authors acknowledge Prof. Diego Berzi, Dr. Dalila
Vescovi, Dr. Irene Redaelli and Dr. Matteo Zerbi for the work
done on this subject and for all the results obtained, that haveallowed to develop the approach briefly outlined in this paper.

3

E3S Web of Conferences 415, 01006 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341501006

References
1.

D. Muir Wood. Canadian Geotechnical Journal
44(11), 1329-1350 (2007).
2. P. Lade, R. Nelson. Int. J. Numer. Anal. Meth.
Geomech. 11(5), 521-542 (1987).
3. Campbell CS. Rapid granular flows. Annu Rev
Fluid Mech 22, 57- 90 (1990).
4. É. Guazzelli, O. Pouliquen. Journal of Fluid
Mechanics, 852, 1 (2018).
5. J. Sun, S. Sundaresan. J Fluid Mech. 682, 590- 616(2011).
6. S. Chialvo, J. Sun, S. Sundaresan. Phys Rev E Stat
Nonlin Soft Matter Phys. 85 (2 Pt. 1), 021305.
(2012)
7. I. Redaelli, C. di Prisco. Int J Numer Anal
Methods Geomech. 43(16). 2525-2550 (2019)
8. I. Redaelli, C. di Prisco. Granul Matter. 23(2), 134 (2021).
9. D. Vescovi, C. di Prisco, D. Berzi. Int J Numer
Anal Methods Geomech. 37(17), 2937-2951(2013).
10. D. Berzi, JT Jenkins. Soft Matter; 11(24), 47994808 (2015).
11. A. Schofield, P. Wroth. Critical State Soil
Mechanics. McGraw-Hill; 1968.
12. M. Jefferies. Géotechnique.43(1), 91-103 (1993).
13. Garzó V, Dufty JW. Phys Rev E. 59(5), 58955911 (1999)
14. MiDi GDR. Eur Phys J E. 14(4), 341- 365 (2004).
15. P.Marveggio, I. Redaelli, C. di Prisco. Challengesand Innovations in Geomechanics. IACMAG
2021. Lecture Notes in Civil Engineering 125,491–497 (2021).
16. T. Shire, C. O’Sullivan, K. Hanley. J. Geotech.
Geoenviron. Eng., 140 (12), p. 04014072 (2014)
17. M. Manzari, Y. Dafalias. Geotechnique. 47(2),255-272 (1997)
18. P. Marveggio, I. Redaelli, C. di Prisco. Int J
Numer Anal Methods Geomech. 46, 2415– 2445
(2022).
19. D. Vescovi, P. Marveggio, C. di Prisco.
Géotechnique. 70(7), 608-620. (2020).
20. M. Zerbi, P. Marveggio, C. di Prisco. Underrevision
21. P. Marveggio, M. Zerbi, C. di Prisco. In
preparation

4