ANALISI DI UN TAGLIO VERTICALE CON IL METODO · PDF fileIncontro Annuale dei Ricercatori di...

6
Incontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012 Cuomo, Cascini, Iannone, Dufour, Darve, Prime 1 ANALISI DI UN TAGLIO VERTICALE CON IL METODO FEMLIP Sabatino Cuomo (1) , Leonardo Cascini (1) , Anna Iannone (1) , Frederic Dufour (2) , Felix Darve (2) , Noemie Prime (2) (1) Università di Salerno, (2) Grenoble-INP / Université Joseph Fourier / Laboratoire 3SR [email protected] , [email protected] , [email protected] , [email protected] , [email protected] , [email protected] Sommario La presente nota illustra le potenzialità di un metodo numerico innovativo FEMLIP (Finite Element Method with Lagrangian Integration Points) che consente la simulazione di grandi deformazioni in analisi di tipo tenso-deformativo. In particolare, dopo una breve introduzione degli aspetti teorici, il metodo FEMLIP è applicato al caso di una coltre superficiale omogenea sede di un taglio verticale. I risultati conseguiti sono, da una parte, congruenti con quelli forniti da metodi più tradizionali per quanto riguarda le condizioni generalmente riconosciute come tipiche dei fenomeni di instabilità di versante e, dall’altra, delineano una molteplicità di scenari, conseguenti alle condizioni cosiddette di equilibrio limite, di sicuro interesse sia in campo scientifico e sia da un punto di vista più squisitamente applicativo. Introduzione I fenomeni di instabilità di versante possono essere caratterizzati da molteplici fasi (pre- rottura, rottura, post-rottura e propagazione) che si differenziano, tra l’altro, per l’entità degli spostamenti del terreno (che possono variare tra pochi mm ad alcuni km nel caso delle frane tipo flusso). Ne consegue che la modellazione geomeccanica di tali fenomeni è, usualmente, riferita alle singole fasi, analizzate separatamente e con approcci differenti. Per l’analisi delle fasi di pre-rottura e rottura si utilizzano prevalentemente approcci di tipo lagrangiano e i metodi numerici più utilizzati sono: i) FEM (Finite Element Method) nell’ipotesi di piccole deformazioni, ii) DEM (Discrete Element Method), e più recentemente iii) SPH (Smooth Particle Hydrodynamics). Tuttavia, importanti limitazioni possono scaturire dall’applicazione di ciascuno di essi in quanto: i) il metodo FEM perde efficacia quando le deformazioni e/o gli spostamenti superano certi limiti, ii) il metodo DEM richiede lungi tempi di calcolo quando il numero di particelle cresce, iii) il metodo SPH pone alcune limitazioni per quanto attiene le condizioni al contorno. Per l’analisi della fase di propagazione si ricorre, invece, ad approcci di tipo euleriano o misti qauli: i) FEM con formulazione euleriana, ii) ALE (Arbitrary Lagrangian Eulerian). Tali approcci non consentono, tuttavia, di portare in conto, in modo rigoroso, le interfacce tra materiali diversi e la storia tenso-deformativa dei singoli materiali. Un’importante differenza tra i metodi in precedenza richiamati riguarda l’associazione tra i punti di calcolo e il materiale (fig. 1) che è caratteristica dei metodi FEM, DEM e SPH; negli altri metodi, tale corrispondenza viene, viceversa, a mancare. Conseguentemente, nel primo caso non si riesce a simulare il flusso di un materiale durante la fase di propagazione; nel secondo non possono adeguatamente analizzarsi materiali il cui comportamento dipende dalla sua storia tenso-deformativa. Da quanto detto scaturisce la necessità di implementare e validare innovativi metodi numerici (fig. 1) che consentano di simulare, in maniera unitaria, il passaggio dalle piccole alle grandi

Transcript of ANALISI DI UN TAGLIO VERTICALE CON IL METODO · PDF fileIncontro Annuale dei Ricercatori di...

Page 1: ANALISI DI UN TAGLIO VERTICALE CON IL METODO · PDF fileIncontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012 Cuomo, Cascini, Iannone, Dufour, Darve,

Incontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012

Cuomo, Cascini, Iannone, Dufour, Darve, Prime 1

ANALISI DI UN TAGLIO VERTICALE CON IL METODO FEMLIP

Sabatino Cuomo (1), Leonardo Cascini (1), Anna Iannone (1), Frederic Dufour (2), Felix Darve (2), Noemie Prime (2)

(1) Università di Salerno, (2) Grenoble-INP / Université Joseph Fourier / Laboratoire 3SR [email protected], [email protected], [email protected], [email protected],

[email protected], [email protected]

Sommario

La presente nota illustra le potenzialità di un metodo numerico innovativo FEMLIP (Finite Element Method with Lagrangian Integration Points) che consente la simulazione di grandi deformazioni in analisi di tipo tenso-deformativo. In particolare, dopo una breve introduzione degli aspetti teorici, il metodo FEMLIP è applicato al caso di una coltre superficiale omogenea sede di un taglio verticale. I risultati conseguiti sono, da una parte, congruenti con quelli forniti da metodi più tradizionali per quanto riguarda le condizioni generalmente riconosciute come tipiche dei fenomeni di instabilità di versante e, dall’altra, delineano una molteplicità di scenari, conseguenti alle condizioni cosiddette di equilibrio limite, di sicuro interesse sia in campo scientifico e sia da un punto di vista più squisitamente applicativo.

Introduzione I fenomeni di instabilità di versante possono essere caratterizzati da molteplici fasi (pre-rottura, rottura, post-rottura e propagazione) che si differenziano, tra l’altro, per l’entità degli spostamenti del terreno (che possono variare tra pochi mm ad alcuni km nel caso delle frane tipo flusso). Ne consegue che la modellazione geomeccanica di tali fenomeni è, usualmente, riferita alle singole fasi, analizzate separatamente e con approcci differenti. Per l’analisi delle fasi di pre-rottura e rottura si utilizzano prevalentemente approcci di tipo lagrangiano e i metodi numerici più utilizzati sono: i) FEM (Finite Element Method) nell’ipotesi di piccole deformazioni, ii) DEM (Discrete Element Method), e più recentemente iii) SPH (Smooth Particle Hydrodynamics). Tuttavia, importanti limitazioni possono scaturire dall’applicazione di ciascuno di essi in quanto: i) il metodo FEM perde efficacia quando le deformazioni e/o gli spostamenti superano certi limiti, ii) il metodo DEM richiede lungi tempi di calcolo quando il numero di particelle cresce, iii) il metodo SPH pone alcune limitazioni per quanto attiene le condizioni al contorno. Per l’analisi della fase di propagazione si ricorre, invece, ad approcci di tipo euleriano o misti qauli: i) FEM con formulazione euleriana, ii) ALE (Arbitrary Lagrangian Eulerian). Tali approcci non consentono, tuttavia, di portare in conto, in modo rigoroso, le interfacce tra materiali diversi e la storia tenso-deformativa dei singoli materiali. Un’importante differenza tra i metodi in precedenza richiamati riguarda l’associazione tra i punti di calcolo e il materiale (fig. 1) che è caratteristica dei metodi FEM, DEM e SPH; negli altri metodi, tale corrispondenza viene, viceversa, a mancare. Conseguentemente, nel primo caso non si riesce a simulare il flusso di un materiale durante la fase di propagazione; nel secondo non possono adeguatamente analizzarsi materiali il cui comportamento dipende dalla sua storia tenso-deformativa. Da quanto detto scaturisce la necessità di implementare e validare innovativi metodi numerici (fig. 1) che consentano di simulare, in maniera unitaria, il passaggio dalle piccole alle grandi

Page 2: ANALISI DI UN TAGLIO VERTICALE CON IL METODO · PDF fileIncontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012 Cuomo, Cascini, Iannone, Dufour, Darve,

Incontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012

Cuomo, Cascini, Iannone, Dufour, Darve, Prime 2

deformazioni. A tale fine, nella presente nota si illustrano brevemente gli aspetti teorici del cosiddetto metodo FEMLIP (Finite Element Method with Lagrangian Integration Points) che è applicato al caso di una coltre superficiale nella quale viene realizzato un taglio verticale non sostenuto (Cuomo et al., 2012).

FEMLIP

Mezzo discontinuo

I punti materiali e quellidi calcolo sono legati

I punti materiali e quellidi calcolo sonoindipendenti

Numero di particellemolto elevato

Discretizzazione di un dominio continuo in sub-

domini usualmente chiamatielementi

FEM

SPH

DEM

Mezzo continuo

Figura 1. Differenti alternative per la modellazione geomeccanica di frane tipo flusso.

Analisi FEMLIP di un taglio verticale

Una breve illustrazione del metodo FEMLIP Il metodo degli elementi finiti con punti di integrazione lagrangiana (Moresi et al, 2003) deriva dal cosiddetto Particle In Cell method (Sulsky et al., 1995) e si basa sulla dissociazione cinematica tra i nodi di calcolo di una maglia agli elementi finiti euleriana (fissa) ed i punti materiali che: i) contengono le informazioni relative al materiale, ii) possono muoversi all’interno del dominio di calcolo, iii) sono utilizzati come punti di integrazione per ciascun elemento della maglia. Ad esempio, per una generica configurazione del materiale (step n di figura 2), i punti materiali sono utilizzati come punti di integrazione per ciascun elemento e dalla risoluzione delle equazioni reggenti nei nodi della maglia è possibile ottenere il campo di velocità del(i) materiale(i). La peculiarità di FEMLIP è rappresentata dal fatto che il campo di velocità è interpolato sui punti materiali che vengono congruentemente spostati, ad ogni step, all’interno del dominio (fig. 2c), consentendo l’analisi delle successive configurazioni deformate (step n+1 di figura 2d). FEMLIP coniuga, dunque, le potenzialità di una maglia euleriana fissa per la descrizione di grandi spostamenti e le potenzialità di un approccio lagrangiano per la modellazione del comportamento meccanico di un geomateriale qual’è il terreno. È importante, inoltre, rilevare

a) b) c) d)

Figura 2. Uno schema relativo al metodo FEMLIP.

Page 3: ANALISI DI UN TAGLIO VERTICALE CON IL METODO · PDF fileIncontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012 Cuomo, Cascini, Iannone, Dufour, Darve,

Incontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012

Cuomo, Cascini, Iannone, Dufour, Darve, Prime 3

che, per un fissato step di calcolo e per una fissata configurazione del materiale, il metodo FEMLIP discretezza le equazioni reggenti (nei nodi) in modo del tutto analogo ad un metodo FEM tradizionale e, quindi, può prendere a riferimento e integrare, in linea di principio, qualunque modello costitutivo attualmente disponibile. Il metodo FEMLIP è implementato nel codice di calcolo ELLIPSIS (Moresi et al., 2003) che risolve le classiche equazioni di conservazione della massa e della quantità di moto in condizioni quasi-statiche, assumendo quali incognite la velocità e lo stato tensionale dei materiali. In particolare, ciascun materiale è schematizzato come mezzo monofase, per il quale la velocità di deformazione totale è data dalla somma di tre componenti: i) elastica, ii) viscosa e iii) plastica. Con riferimento al legame costitutivo, nel codice è implementato il cosiddetto modello "Plasol" (Barnichon, 1998) che si caratterizza per: i) un criterio di rottura di Van Eekelen (VE) (van Eekelen, 1980), non molto differente dal criterio di Drucker-Prager (DP), ii) incrudimento di tipo isotropo, iii) legge del flusso plastico non associata.

Analisi numeriche Di seguito si propone il caso di una coltre superficiale (figura 3) - con inclinazione di 30°, spessore di 3 m e con proprietà meccaniche analoghe a quelle dei terreni piroclastici della regione Campania (Bilotta et al., 2005) - nella quale si realizza un taglio verticale non sostenuto di altezza rispettivamente pari 2,5 m o 1,5 m. Il terreno è schematizzato come materiale elastico perfettamente plastico con criterio di rottura di Van Eekelen per rendere possibile il confronto dei risultati di FEMLIP con quelli forniti da metodi semplificati; si assume, inoltre, che la legge del flusso plastico sia non associata, con angolo di dilatanza pari a 10°. Durante tutto il processo di deformazione del pendio, vale a dire dalla fase di pre-rottura fino a quando la massa mobilitata si arresta, sono assunte condizioni drenate e pressioni neutre pari a zero. Infine, il dominio di calcolo è schematizzato attraverso una griglia di elementi quadrangolari bi-lineari di dimensione pari a 0,1 m. Per valutare le potenzialità di FEMLIP nel descrivere il processo di deformazione del pendio è stata condotta un’analisi parametrica con valori dei parametri di resistenza a rottura riportati in tabella 1. Per ciascun caso sono state inizialmente applicate le forze di gravità ed è stato poi introdotto un taglio verticale, rimuovendo la corrispondente porzione del dominio di calcolo. In figura 4 si riportano i risultati ottenuti per il caso 7 di tabella 1 che si caratterizza per la formazione di una zona all’interno della quale le velocità di deformazione plastica si incrementano nel tempo e determinano uno scenario (0-4, in figura 4a) nell’ambito del quale il materiale si propaga verso valle (4-6, in figura 4a) dando luogo ad un nuova configurazione di equilibrio (6). In figura 4b si riporta, inoltre, il diagramma di spostamenti di un punto (P) posto alla base del taglio verticale che, congruentemente con quanto esposto in precedenza, è caratterizzato da una fase di accelerazione (0-4, fase di rottura) e da una successiva fase di decelerazione (4-6, fase di propagazione) fino al completo arresto. Analoghe analisi FEMLIP sono state condotte per un’ampia casistica i cui risultati possono essere efficacemente descritti facendo riferimento a un parametro adimensionale c’/ sat×hcut,

30°

2.5 m g

a) b)

20 m

1 m

3 mz

x

P (10, 1)

0

A

B

C

Condizioni al contorno.:0A, BC: xdispl = 0, 0B: xdispl = zdispl = 0

g

Figura 3. Schema geometrico del pendio considerato (modificata da Cuomo et al., 2012).

Page 4: ANALISI DI UN TAGLIO VERTICALE CON IL METODO · PDF fileIncontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012 Cuomo, Cascini, Iannone, Dufour, Darve,

Incontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012

Cuomo, Cascini, Iannone, Dufour, Darve, Prime 4

Caso 1 2 3 4 5 6 7 8 9 10 11 12 13hcut (m) 1.5 1.5 1.5 1.5 1.5 1.5 2.5 2.5 2.5 2.5 2.5 2.5 2.5c’ (kPa) 1 1 2 2 3 3 2 2 3 3 4 5 5ϕ‘ (°) 32 36 32 36 32 36 32 36 32 36 36 32 36g =13 kN/m3, E = 5000 kPa, n = 0.3, y = 10°

Tabella 1. Configurazioni geometriche e parametri meccanici assunti per la modellazione geomeccanica.

0 1.30E+003

Legend

Equivalent Plastic strain rate

0

4

6

5

2 3

1

(0-4: rottura, 4-6: propagazione)

0.00

0.10

0.20

0.30

0.40

0.50

0 30 60 90 120xd

ispl (

m)

time step

12

3

4 5 6

Figura 4. Velocità di deformazione plastica (a), curva spostamenti-tempo (b) per il caso 7 di tabella 1 (modificata da Cuomo et al., 2012).

che, in questo caso, che può considerarsi una misura del rapporto tra la resistenza del terreno e lo stato tensionale mobilitato all’interno del pendio. La Figura 5 evidenzia che per valori decrescenti del fattore c’/ sat×hcut si individuano quattro principali scenari: i) pendio stabile (caso 5), ii) rottura di tipo localizzata che mobilita un volume di terreno di modesta entità con spostamenti ridotti (caso 3) iii) rottura localizzata che instabilizza volumi di maggiori dimensioni sede di spostamenti più rilevanti (caso 9), iv) rottura di tipo generalizzata in un volume di grandi dimensioni che, a differenza dei casi precedenti, non perviene ad una nuova configurazione geometrica di equilibrio. È appena il caso di notare che il campo di deformazioni plastiche e di spostamenti evolve progressivamente nel tempo a conferma della necessità di analizzare il processo di instabilità del pendio come un unico processo di tipo progressivo (Cuomo et al., 2012).

0.000

0.002

0.004

0.006

0.008

0.010

0 50 100 150

xdisp

l·c' / g

tot·h

cut2

time step

5397

35 79

0 1.30E+003

Legend

Equivalent Plastic strain rate

Case 5 3 9 7

c’/gtot · hcut 0.154 0.103 0.092 0.062

tipo di rottura no localizzata localizzata diffusa

volume mobilitato - piccolo medio grande

nuova configurazione geometrica del pendio

- si si no

cuttot hc /'

a)

b)c)

Figura 5.Differenti scenari ottenuti al variare dei parametri geometrici e meccanici

(modificata da Cuomo et al., 2012).

Page 5: ANALISI DI UN TAGLIO VERTICALE CON IL METODO · PDF fileIncontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012 Cuomo, Cascini, Iannone, Dufour, Darve,

Incontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012

Cuomo, Cascini, Iannone, Dufour, Darve, Prime 5

Confronto con metodi tradizionali

Confronto con il metodo all’equilibrio limite (LEM) I risultati ottenuti con il metodo FEMLIP sono stati confrontati con quelli ottenibili con metodi più tradizionali, quali i metodi dell’equilibrio limite e il metodo FEM tradizionale; tale confronto ha avuto la duplice finalità di validare i risultati FEMLIP, ove possibile, e di consentire una stima delle differenze tra i risultati di FEMLIP e quelli forniti da metodi più semplificati. La Figura 6 mette in luce che le analisi all’equilibrio limite (svolte con l’ausilio dei metodi di Janbu, Mongenstern e Price e facendo riferimento ad un criterio di resistenza tipo Drucker-Prager), forniscono nel caso esaminato un fattore di sicurezza superiore all’unità per molteplici combinazioni dei parametri geometrici e meccanici; analogamente, le analisi FEMLIP individuano campo di spostamenti moderati senza che siano simulati processi di rottura del pendio. Al fine di confrontare in modo quantitativo tali risultati, si è fatto riferimento a un parametro geometrico-meccanico (c’/ sat×hcut×tanf’) che contempla tutti i principali fattori che governano le condizioni di stabilità di un pendio di fissata inclinazione. In figura 6, in particolare, si riportano i casi di tabella 1, suddivisi in due gruppi corrispondenti rispettivamente all’insorgere o meno di fenomeni di rottura. Con riferimento a tale grafico si osserva che: i) i punti sono ben interpretati da linee di regressione lineari e con bassa pendenza, ii) i metodi all’equilibrio limite (LEM) forniscono valori più elevati, con differenze di circa il 10-30%, del parametrio c’/ sat×hcut×tanf’. Se ne deduce, rispettivamente che: i) la coesione gioca, come è noto, un ruolo fondamentale per il mantenimento delle condizioni di stabilità, ii) i risultati dei metodi LEM sono leggermente cautelativi rispetto a quelli di FEMLIP in virtù dei diversi stati tensionali e delle diverse ipotesi assunte a riferimento nei suddetti metodi. Si sottolinea, inoltre, che le grandezze adimensionali in precedenza citate sono analoghe a quelle individuate da Taylor (1937) e Michalowski (2002) per la costruzione delle cosiddette "carte di stabilità", che si riferiscono al caso di un taglio verticale/inclinato con piano campagna orizzontale. I risultati conseguiti sono da considerarsi, pertanto, di un certo interesse nell’ambito della letteratura scientifica attualmente disponibile sull’argomento.

R² = 0.9491

R² = 0.8143R² = 0.8920

R² = 0.6384

0.00

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.55 0.6 0.65 0.7 0.75

c' / g

tot·h

cut·t

anϕ

' (-)

tan ϕ' (-)

no failure - LEMno failure - FEMLIPfailure - LEMfailure - FEMLIP

Figura 6.Confronto tra i risultati di FEMLIP e i metodi dell’equilibrio limite (LEM), (Cuomo et al., 2012).

Confronto con metodo FEM La Figura 7, infine, confronta i risultati di FEMLIP con quelli ottenuti tramite un codice FEM commerciale Sigma/W (GeoSlope, 2005), facendo riferimento alle stesse configurazioni geometriche, ad una analoga schematizzazione spaziale del dominio (maglia non-strutturata di

Page 6: ANALISI DI UN TAGLIO VERTICALE CON IL METODO · PDF fileIncontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012 Cuomo, Cascini, Iannone, Dufour, Darve,

Incontro Annuale dei Ricercatori di Geotecnica 2012 - IARG 2012 Padova, 2-4 Luglio 2012

Cuomo, Cascini, Iannone, Dufour, Darve, Prime 6

elementi lineari triangolari di dimensione non superiore a 0,2 m) e ad un analogo modello costitutivo (criterio di rottura tipo Drucker-Prager, con legge del flusso plastico non associata). Dal confronto (fig. 7) si osserva che i metodi FEM e FEMLIP forniscono risultati molto prossimi tra loro per elevati valori della coesione; congruentemente, i risultati differiscono in misura sempre maggiore al diminuire della coesione stessa. Inoltre, quando si raggiungono le condizioni di rottura, cioè quando i metodi LEM forniscono fattori di sicurezza inferiori all’unità, il metodo FEMLIP è ancora in grado di fornire gli spostamenti nel campo delle grandi deformazioni, mentre il metodo FEM tradizionale non contempla tale possibilità.

0.0

0.5

1.0

1.5

2.0

2.5

3.0

0.00 0.05 0.10 0.15 0.20 0.25

z (m

)

x-Displacement (m)

case 93

Results at vertical cut- Classical FEM (a,b)- FEMLIP (1,2, 3)

case 12 case 11a,1 2,b

LEM

Case c’(kPa)

ϕ‘(°)

Factorof safety(-)

12 5 32 1.05

11 4 36 1.01

9 3 32 0.79

Figura 7. Spostamenti orizzontali del punto P di figura 3b calcolati tramite FEMLIP ovvero FEM e confronto con

i risultati dei metodi all’equilibrio limite (LEM) (modificata da Cuomo et al., 2012).

Considerazioni conclusive Nella presente comunicazione è stata illustrata una prima semplice applicazione di un innovativo metodo numerico denominato FEMLIP al caso estremamente schematico di una coltre superficiale sede di un taglio verticale non sostenuto. I risultati conseguiti appaiono particolarmente interessanti in quanto delineano differenti scenari di rottura e di propagazione del volume instabilizzato che possono essere messi in relazione a grandezze adimensionali che dipendono dalla geometria del problema e dalle caratteristiche meccaniche dei materiali coinvolti.

Bibliografia

Barnichon, J.D. (1998). Finite Element Modelling in Structural and Petroleum Geology. PhD thesis, Université de Liège, 264 p.

Bilotta, E., Cascini, L., Foresta, V., Sorbino, G. (2005). Geotechnical characterization of pyroclastic soils involved in huge flowslides. Geotechnical and Geological Engineering, 23, 365-402.

Cuomo, S., Prime, N., Iannone, A., Dufour, F., Cascini, L., Darve, F. (2012). Large deformation FEMLIP drained analysis of a vertical cut. Acta Geotechnica, DOI: 10.1007/s11440-012-0179-2.

Michalowski, R.L. (2002). Stability Charts for Uniform Slopes, ASCE's Journal of Geotechnical and Geoenvironmental Engineering, 4, 351-355.

Moresi, L. N., Dufour, F. and Muhlhaus, H.-B. (2003). A Lagrangian integration point finite element method for large deformation modelling of viscoelastic geomaterials. Journal of Computational Physics 184, 476-497.

Sulsky, D., Zhou S.-J. and Schreyer H.L. (1995). Application of a particle-in-cell method to solid mechanics. Comput. Phys. Comm. 87, 236-252.

Taylor, D.W (1937). Stability of earth slopes. J. Boston Soc. Civil Eng., 243. Reprinted in: Contributions to Soil Mechanics 1925 to 1940, Boston Society of Civil Engineers, 337–386.

Van Eekelen H.A.M. (1980). Isotropic yield surfaces in three dimensions for use in soil mechanics. International Journal for Numerical and Analythical Methods in Geomechanics, 4(1), 89-101.