POLITECNICO DI TORINO
I Facolt`a di Ingegneria Corso di Laurea in Ingegneria Aerospaziale
Tesi di Laurea Magistrale
Stability and applications of the peridynamic method
Relatori: Prof. Marco Di Sciuva Prof. Paolo Maggiore Candidato: Matteo Polleschi
Luglio 2010
Ringraziamenti Il mio ringraziamento va ai miei relatori, al Professor David Steigmann dell’Universit`a di Berkeley che mi ha accolto nel mio periodo statunitense, ai miei compagni in questa lunga avventura e tutti coloro che mi sono vicini, fisicamente e con lo spirito, nonostante la mia anima bipolide e cosmopolita.
1
Sommario La meccanica strutturale `e portata ad affrontare con modelli numerici sempre pi` u affidabili la realt`a che ci circonda; dai primi modelli agli elementi finiti, limitati dalla scarsa potenza computazionale a disposizione, si `e passati a modelli estremamente accurati in grado, con indicazioni certificate, di sostituire la validazione sperimentale delle strutture aeronautiche e spaziali. Si va quindi dai calcoli statici a quelli dinamici e di buckling; una branca da non sottovalutare sia per scopi di sicurezza che di ottimizzazione di prestazioni, materiali e processi di produzione `e quella della meccanica della frattura e dell’impatto. Questi casi sono pi` u complessi da trattare in quanto, in seguito all’evento dinamico posto in essere, si creano delle discontinuit`a nel materiale; essendo le equazioni della meccanica classica basate su derivate spaziali, in presenza delle suddette discontinuit`a queste grandezze non sono pi` u calcolabili dando origine a delle singolarit`a. Per ovviare a questo esistono vari metodi, dall’affinare la mesh dove il materiale diventa discontinuo all’introdurre delle ipotesi sul comportamento della frattura. Esiste un modello di meccanica del continuo che ovvia a queste limitazioni, in quanto si basa su equazioni integro-differenziali che non vanno a coinvolgere nessuna derivata spaziale, utilizzando invece gli integrali, che sono in grado di trattare anche questo tipo di problema. Questa teoria si basa su degli studi sul comportamento dei cristalli di Rogula e Kunin (intorno al 1982), poi ripresi e rielaborati dal fisico statunitense Stewart Silling nella sua pubblicazione [1], e prende il nome dalle due parole greche perÐ e dÔnamic, cio`e “forza dell’intorno”; `e un argomento di ricerca molto recente, ma ha tutte le potenzialit`a per essere utilizzato in larga scala. Il lavoro principale svolto per questa tesi `e consistito nella stesura di un codice c++ in grado di risolvere problemi di frattura ed impatto, implementando un algoritmo di integrazione misto esplicito-implicito per ottenere la stabilit`a numerica (che affligge il modello da una parte in quanto rischia di far “esplodere ”il calcolo, dall’altra perch´e prendendo troppo margine si rallenta il calcolo stesso); successivamente, lo si `e applicato ad alcuni casi semplici a sufficienza per essere lanciati senza avere a disposizione un hard1
ware particolarmente performante, ma significativi per dimostrare che questo metodo `e in grado di descrivere la propagazione delle cricche senza alcuna ipotesi aggiuntiva, come diretta conseguenza di geometrie, materiali, carichi e condizioni al contorno. Segue un riassunto dei vari capitoli. Capitolo 1: fondamenti teorici L’equazione del moto del modello peridinamico `e: Z ρ(x)¨ u(x, t) = f (u(x0 , t) − u(x, t), x0 − x)dVx0 + b(x, t) (1) R
dove le variabili hanno il significato di1 : • ρ(x): densit`a in x • x e x0 : punti (particelle) nella configurazione di riferimento • u: vettore spostamento • f : pairwise force function • b: campo di forze esterno Questa equazione `e una forma della seconda legge di Newton. Come si pu`o notare, l’equazione `e integro-differenziale, e la differenziazione `e effettuata nel tempo, dunque non vi `e la presenza delle derivate spaziali che causano problemi sulle discontinuit`a. Si introduce il vettore posizione relativa ξ ξ = x0 − x e il vettore spostamento relativo η η = u(x0 , t) − u(x, t) = u0 − u La quantit`a η + ξ descrive quindi la posizione relativa attuale tra le due particelle. Nella formulazione pi` u generale, l’integrazione `e effettuata su tutto il corpo. Tuttavia `e conveniente definire come in figura [1] una regione all’interno del dominio delimitata da un raggio δ tale che: if kξk ⇒ δ =⇒ f (η, ξ) = 0 ∀η 1
I vettori sono indicati in grassetto
2
(2)
R f δ
x'
x
Figura 1: Dominio ed orizzonte
Questa quantit`a `e chiamata orizzonte, a ricordare che `e la regione all’interno della quale la particella “vede”altre particelle esercitare una forza su di essa. Per quanto riguarda il valore di δ, `e prassi utilizzare un valore che `e tre volte la spaziatura tra i nodi in cui viene discretizzato il dominio. Se si utilizza un valore inferiore, le fratture tendono a seguire un percorso innaturale, lungo le righe o colonne della mesh; al contrario, se si utilizza un orizzonte troppo ampio, non solo si appesantisce molto il calcolo, ma si causa un’eccessiva dispersione ondosa ed il solido inizia a comportarsi come un fluido. Il modello peridinamico `e non locale, nel senso che le particelle interagiscono a distanza finita, come nella dinamica molecolare. Esistono molte teorie di dinamica molecolare, ma questa `e l’unica ad eliminare le derivate spaziali. In questa tesi si `e utilizzato l’approccio basato sui collegamenti particellaparticella, che ha alcune limitazioni (come il non poter utilizzare modelli costitutivi dalla teoria classica ed essere limitato a materiali con coefficiente di Poisson ν = 0.25): per risolvere queste limitazioni `e stato creato il modello dei Peridynamic States, che permette il collegamento della densit`a di energia di deformazione con il cambiamento di volume locale, e non solo l’interazione tra due sole particelle. Le propriet`a di legame non dipendono solo dall’allungamento del legame stesso, ma anche dagli allungamenti di tutti gli altri legami connessi alla particella che si sta considerando; questo permette di definire modelli costitutivi pi` u generali. Il termine di forza nell’equazione 1 `e un funzionale che dipende dalla funzione vettoriale f , chiamata la pairwise force function (PFF). Essa fornisce
3
la forza per unit`a di volume al quadrato dovuta ad una particella x0 su una particella x. Nella meccanica del continuo classica questo funzionale `e la divergenza del tensore di stress. Le propriet`a del materiale ed il suo comportamento sono descritte dalla PFF. Con una serie di considerazioni basate sull’equilibrio e della conservazione della quantit`a di moto e del momento angolare, si ottiene che la forza tra due particelle dipende solo dalla loro posizione relativa attuale, e che f `e una funzione dispari. L’allungamento di un legame `e dato da: s=
y − |ξ| |η + ξ| − |ξ| = |ξ| |ξ|
(3)
dove y = |η + ξ| `e il modulo del vettore posizione relativa attuale. Si possono caratterizzare un gran numero di materiali in funzione dell’allungamento; in questo lavoro si `e scelto un materiale proporzionale (lineare), che abbia un comportamento di questo tipo: f ∝ cs
(4)
dove c `e una costante elastica che dipende dal bulk modulus K (caratteristica macroscopica del materiale) secondo la formula: c=
18K πδ 4
(5)
Il danneggiamento del materiale `e introdotto a livello di legame: se l’allungamento supera una certa soglia, il legame si rompe e non restituisce pi` u alcuna forza. Anche il massimo allungamento del legame si pu`o legare a parametri macroscopici: r 5G0 s0 = (6) 9Eδ Entrambe queste definizioni derivano dal lavoro [5] di Silling ed Askari. La teoria della peridinamica non considera, nella sua formulazione di base, forze tra particelle che all’inizio del calcolo erano pi` u lontani dell’orizzonte; questo impedirebbe la trattazione dell’impatto (se non sotto certe ipotesi limitative, assegnando una forte velocit`a iniziale alla parte impattata). Per ovviare a questo problema sono state introdotte delle forze a corto raggio, di natura esclusivamente repulsiva (e ancora una volta dirette secondo la congiungente delle particelle), la cui intensit`a `e data da 0 |y − y| fs = min 0, cs −1 (7) 2rs 4
dove y = x + u e y0 = x0 + u0 sono le posizioni correnti di x ed x0 , cs e rs costanti positive. Le forze repulsive compaiono quando le particelle sono pi` u vicine di 2rs , che `e chiamato il raggio del nodo. Di solito si assume che cs = 15c. Capitolo 2: metodo numerico Il dominio `e suddiviso in nodi equispaziati, e nell’equazione del moto l’integrale diventa una somma finita: X (8) ρu¨ni = f (unp − uni , xp − xi )Vp + bni p
essendo uni = u(xi , tn )
(9)
Vp `e il volume del nodo p, che per una griglia di calcolo uniforme `e semplicemente pari a x3 . La somma `e effettuata su tutti i nodi p tali che |xp −xi | < δ. L’equazione del moto discretizzata `e facile da implementare, anche nella sua forma non lineare. In ogni caso, la versione linearizzata diventa: X C(unp − uni )(xp − xi )Vp + bni (10) ρu¨ni = p
e per il termine d’accelerazione si utilizzano le differenze finite centrate un+1 − 2uni + uin−1 i n ¨ ρui = ∆t2
(11)
L’analisi di stabilit`a di von Neumann porta al requisito per il passo di integrazione di s 2ρ (12) ∆t < P p Vp |C(xp − xi )| Questo requisito, per`o, da una parte pu`o non essere valido per prolemi che si discostino dalla linearit`a, e dall’altra porta ad un calcolo che non `e completamente efficiente: il passo di integrazione non risulta essere ottimizzato, e si va a sprecare del tempo durante il calcolo. Per questo motivo, si `e introdotto un metodo che permetta di valutare “in itinere”la bont`a del passo di integrazione; questo si basa su una valutazione dell’errore che `e fondata sul moto delle particelle, o meglio sulla loro posizione: la nuova posizione calcolata non pu`o discostarsi esageratamente da quella precedente (ragionando sulla somma dei valori assoluti degli spostamenti). In questo modo le posizioni dei nodi tendono a variare in maniera graduale, evitando le “esplosioni”tipiche degli errori numerici. 5
Alla base della correttezza dell’analisi numerica che si va a sviluppare c’`e comunque la precisione di macchina, storicamente indicata come : essa `e definita come il pi` u piccolo numero tale che la differenza tra 1 e 1 + `e non nulla. Quando si ha a che fare con numeri troppo piccoli, vicini alla precisione di macchina, si presenta un alto rischio di errore numerico, che porta ad avere errori nel computo delle forze anche di diversi ordini di grandezza, fornendo una soluzione con forze “numeriche”molto intense e che non sono equilibrate. Questo `e dovuto al calcolo dell’allungamento, per il quale spesso la differenza tra valori cos`ı piccoli dall’essere vicini all’ `e la criticit`a principale. Per ovviare a questi problemi, bisogna controllare tutti i valori che potrebbero avvicinarsi troppo alla precisione di macchina e, nel caso questo si verifichi, portarli ad annullamento. Ritornando al metodo numerico, l’approccio utilizzato in questo lavoro si ispira al lavoro del Professor Zohdi dell’Universit`a della California, Berkeley; la sua prima formulazione `e stata sviluppata per la soluzione di problemi multicampo termochimici. Il metodo che ne risulta `e misto esplicito-implicito: l’integrazione procede esplicita attraverso i passi temporali, ma ad ogni passo d’integrazione si attiva un ciclo implicito che ha il compito di verificare se il passo d’integrazione stesso sia quello ottimale per rimanere nei limiti della stabilit`a, ma non perdendo in efficienza computazionale. A tal proposito, si utilizzano le seguenti quantit`a: PN n+1,K − ui n+1,K−1 || def i=1 ||ui (13) ωK = PN n+1,K − u n || ||u i i i=1 che `e l’errore normalizzato, ed inoltre: def
ZK =
ωK T OL
1 T OL pKd ω0 1 ωK pK ω0
(14)
def ΦK =
L’algoritmo utilizzato consiste in questa sequenza di operazioni: • per tutti gli N nodi,
6
(15)
– calcolare la nuova posizione come ui n+1,K ≈ ∆tu˙ni + uni
∆t2 f (ui n+1,K−1 ) + mi
– calcolare le nuove forze (del ciclo interno, da memorizzare in variabili temporanee separate) • calcolare gli errori definiti precedentemente • se la tolleranza `e rispettata (ZK ≤ 1) e (K < Kd ), allora – incrementare il tempo t = t+∆t e ripartire con la prima istruzione – costruire un nuovo time step ∆t = ΦK ∆t • se la tolleranza non `e rispettata (ZK > 1) e K = Kd allora – costruire un nuovo time step ∆t = ΦK ∆t – ricominciare il calcolo dal tempo t Il metodo `e quindi caratterizzato dai parametri T OL, Kd e p. Valori tipici sono di T OL = 0.001, Kd = 20 e p = 3, ma, come si potr`a notare nei calcoli presentati, nel caso della peridinamica questi valori possono essere diversi. In sostanza, il metodo controlla quanto i nodi (particelle) si muovono: se la differenza tra due istanti successivi soddisfa la tolleranza, l’intervallo temporale viene aumentato per migliorare l’efficienza; se invece si `e fuori tolleranza, il ∆t viene ridotto per soddisfare la stabilit`a. Capitolo 3: descrizione dell’architettura software L’obiettivo di questa tesi `e lo sviluppo di un codice peridinamico stabile, ma per ottenere uno strumento flessibile e riutilizzabile `e necessario considerare tutte e tre le fasi dell’analisi: pre-processamento, soluzione e post-processamento. Da un punto di vista operativo, si `e realizzata ed integrata una piattaforma software in grado di svolgere tutte e tre queste mansioni, completando il solutore in c++ con altri software/linguaggi di scripting preesistenti e basati sulla filosofia open source. L’obiettivo finale, da un punto di vista software, `e la creazione di una suite integrata per il calcolo peridinamico che sia in grado di affrontare la pi` u grande variet`a possibile di geometrie e materiali, fornendo risultati di comoda interpretazione. Questo `e stato sviluppato mediante le varie componenti che andiamo ad elencare.
7
Pre-processamento La geometria viene generata con il componente CAD della suite open source Salom`e; la griglia di calcolo viene creata basandosi su questa geometria, grazie al meshatore ancora incluso in Salom`e. Questo passo `e necessario per ottenere una nuvola di nodi che vada a discretizzare il continuo, ma nella struttura di file che viene prodotta come output `e presente anche la definizione degli elementi. A tal proposito si utilizza Impact, un software libero per calcolo esplicito di simulazioni dinamiche; non essendo interessati alle sue capacit`a di processamento, si sfrutta la sua comoda interfaccia grafica per eliminare i suddetti elementi, lasciando solamente i nodi, necessari al calcolo peridinamico. Sempre in questo ambiente si applicano i carichi, si assegnano i vincoli e le velocit`a iniziali. Un ultimo passo viene svolto per rendere il file di output di Impact leggibile dal solutore: `e questa l’ultima fase, detta di parsing, in cui viene prodotto il file con la struttura utilizzata in input dal codice c++. Analizzando le varie fasi, si intuisce come sia possibile trattare un numero virtualmente infinito di geometrie, grazie alla combinazione cad-meshatore; particolare attenzione va posta nel cercare di ottenere una mesh equispaziata, in quanto l’irregolarit`a della griglia porta ad instabilit`a pi` u o meno accentuate. E’ questo uno dei filoni di ricerca immediati in campo peridinamico. Soluzione Il cuore della tesi `e lo sviluppo del solutore, che implementa il metodo esplicito-implicito, il parsing dell’input e la produzione dei risultati. Il codice `e scritto in c++, linguaggio creato da Bjarne Stroustrup ai Bell Labs nel 1983; la scelta `e basata sia sulla sua diffusione nel calcolo scientifico grazie alle alte performance, ma anche per la sua versatilit`a. Esso implementa due librerie: Armadillo per l’algebra lineare (sviluppata in un centro di ricerca australiano, e dotata di licenza open source GNU LGPL), e VTK per la parte grafica (sviluppata da Kitware e fornita con licenza open source). Entrambe sono multipiattaforma, e non utilizzando il codice librerie specifiche di sistema operativo `e dunque portabile. Lo sviluppo `e stato comunque effettuato su piattaforma GNU/Linux (Ubuntu e Mandriva). La struttura del codice `e prevalentemente basata sulla programmazione ad oggetti, infatti nodi e bonds sono sviluppati come classi con propri metodi; l’algoritmo `e schematizzato nelle figure 2 e 3.
8
Start Read input parameters
Read input file
Create nodes
Create links
Peridynamic iterations no
Iterations performed? yes Stop Figura 2: Diagramma di flusso del programma principale
9
Velocity Verlet step 1 and 2
Stable time step calculation
Bond force update
Nodal force computation
Velocity Verlet step 3
Damping correction for velocities
Broken links dump
Output solution at current time step Figura 3: Diagramma di flusso delle iterazioni
10
Post-processamento Per consentire sia un’analisi in tempo reale del calcolo che sta avvenendo, sia la possibilit`a di produrre report, il solutore produce due tipi di file diversi. Il primo tipo sfrutta le librerie VTK implementate per fornire un file nel formato VTK polydata, da aprire all’interno di software per post-processamento con VisIt o Paraview. Questo permette di verificare, durante l’esecuzione di un calcolo, che l’esito vada come previsto, confermando le condizioni iniziali ed al contorno. L’altro output prodotto `e un file geometria per GMSH, un meshatore open source con capacit`a di pratica visualizzazione e produzione immagini (anche vettoriali). La frequenza con cui viene prodotto l’output `e un parametro di lancio del calcolo, in quanto inficia le prestazioni in maniera non trascurabile (sia per via del tempo di scrittura, sia perch`e in quel momento non si possono utilizzare le capacit`a di calcolo parallelo). Capitolo 4: stato dell’arte In questo capitolo viene descritto brevemente qual `e lo stato dell’arte della ricerca sulla fluidodinamica, in particolare quali software sono stati sviluppati e su che cosa si basano. Successivamente vengono presentati alcuni risultati ottenuti in altre pubblicazioni, che confermano ed affiancano il lavoro svolto in questa tesi: sono calcoli di fratture e impatti svolti su piattaforme hardware di potenza rilevante e quindi in grado di trattare un alto numero di nodi, che mostrano come la peridinamica sia in grado di modellare tali problemi in maniera affidabile. I due programmi attualmente esistenti sono due: EMU e LAMMPS. Per quanto riguarda il primo, si fonda sulla bond-based versione della peridinamica (quindi non utilizza i peridynamic states) ma non `e liberamente accessibile: solo un limitato numero di utenti e centri di ricerca degli Stati Uniti ne `e in possesso, e dunque si pu`o solamente verificare attraverso le pubblicazioni per le quali viene utilizzato. LAMMPS invece `e un software molto ampio di dinamica molecolare, liberamente accessibile (fornito dai Sandia Labs), all’interno del quale vi `e un modulo dedicato alla peridinamica. La sua esecuzione ha tempi di calcolo estremamente ottimizzati, ma ha la limitazione di poter descrivere soltanto geometrie semplici. Con entrambi i software, comunque, sono state effettuate simulazioni che mostrano un perfetto accordo con la realt`a che scaturisce esclusivamente da vantaggi del metodo, senza bisogno di ulteriori affinamenti o postprocessamenti. Capitolo 5 e 6: simulazioni e conclusioni Nell’ultima parte di questo scritto sono presentati i risultati delle simulazioni. Va detto che si `e 11
reso necessario affrontare casi relativamente semplici, limitando il numero dei nodi per via del fatto che l’architettura hardware a disposizione (laptop biprocessore) non era sufficiente per modelli di pi` u larga scala (basti pensare che le pubblicazioni sugli stessi argomenti riportano calcoli fatti su cluster di supercomputer), tuttavia i risultati ottenuti sono in ogni caso apprezzabili. L’obiettivo di rendere il metodo stabile e di verificare le potenzialit`a della peridinamica di cogliere fenomeni dinamici come la frattura senza ulteriori ipotesi aggiuntive pu`o dirsi certamente raggiunto: in ogni calcolo che si va a lanciare, qualsiasi sia il time step da cui si parte, questo si assesta immediatamente al valore che soddisfa la stabilit`a, andando a raffinare ed ottimizzare immediatamente dopo. L’effetto del metodo si pu`o conferare visivamente osservando le simulazioni della membrana oscillante: per via della dissipazione introdotta, questa ritorna nella posizione iniziale di equilibrio senza denotare problemi numerici. Inoltre, si pu`o notare come l’evoluzione della frattura fragile (conseguenza del modello di materiale implementato) sia ben descritta dall’evolvere del calcolo: sia nel caso della piastra forata che del provino, infatti, il comportamento descritto `e quello che corrisponde alla teoria ed ai casi sperimentali. Anche l’impatto `e un campo promettente, sebbene sia computazionalmente ancor pi` u dispendioso: infatti, in questo caso, compaiono anche le forze a corto raggio che dipendono dalla posizione corrente di tutte le particelle, aspetto che implica il ricalcolo delle vicinanze di ogni nodo della griglia con una certa frequenza.
12
Indice 1 Peridynamic theory 1.1 Equation of motion . . . . . . 1.1.1 Pairwise force function 1.1.2 Damage . . . . . . . . 1.1.3 Short range forces . . . 1.1.4 Peridynamic states . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
16 16 18 21 21 22
2 Numerical method 24 2.1 Direct numerical method . . . . . . . . . . . . . . . . . . . . . 24 2.2 Time-adaptive method . . . . . . . . . . . . . . . . . . . . . . 25 3 Computational code 3.1 Solver . . . . . . . . 3.1.1 Parallelization 3.2 Pre-processing . . . . 3.2.1 Salome . . . . 3.2.2 Impact . . . . 3.2.3 Python . . . . 3.3 Post-processing . . . 3.3.1 VisIt . . . . . 3.3.2 GMSH . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
28 28 32 33 34 35 36 36 36 37
4 State of the art 4.1 Software . . . . . . . . . . . . . . . . . . 4.1.1 LAMMPS . . . . . . . . . . . . . 4.1.2 EMU . . . . . . . . . . . . . . . . 4.2 Simulations . . . . . . . . . . . . . . . . 4.2.1 Perforation of ductile targets . . . 4.2.2 Perforation of reinforced concrete 4.2.3 Perforation of a fabric . . . . . . 4.2.4 Ball explosion . . . . . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
. . . . . . . .
38 38 38 44 45 45 45 46 46
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
13
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
. . . . . . . . .
4.2.5
Time rate influence . . . . . . . . . . . . . . . . . . . . 47
5 Simulations 5.1 Membrane . . . 5.2 Plate with hole 5.3 Traction test . . 5.4 Impact . . . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
48 48 60 65 68
6 Conclusions
72
Bibliografia
73
14
Introduction Structural mechanics has often to deal with discontinuities, such as cracks and interfaces’ delamination. The classical models are mainly based upon equations, like equilibrium ones, involving partial derivatives with respect to the spatial coordinates; consequently, they’re not suitable to describe discontinuities, wherein derivatives are not defined. Currently used solution methods require knowledge of the crack path, to redefine the mesh alongside. This is obviously not the most general approach. A revolutionary approach was proposed by Stewart Silling in 1998 [1], with a brand new theory called Peridynamics. The name comes from two greek words, perÐ and dÔnamic, meaning “force of the surroundings”, from the kind of particles’ interaction. The funding aspect of his theory is that the model is based upon integral equations, which can be directly applied on every kind of calculation dominium. The core work of this thesis is the implementation of a time-adaptive peridynamic solver using c++ programing language, applying it to simple but relevant cases to test its stability and precision. In the next chapters an overview of the method will be presented, along with calculations’ results and evaluation of industrial applicability.
15
Capitolo 1 Peridynamic theory 1.1
Equation of motion
Considering a body occupying a domain R, the equation of motion for every particle (at point x in the reference configuration and time t) inside the domain is given by: Z ρ(x)¨ u(x, t) = f (u(x0 , t) − u(x, t), x0 − x)dVx0 + b(x, t) (1.1) R
where the symbols define1 • ρ(x): density at x • x and x0 : points in the reference configuration • u: displacement vector • f : pairwise force function • b: body force density field This is the fundamental integro-differential equation of peridynamic theory, based on Newton’s second law. It involves an integral part and a differentiation w.r.t. time, thus avoiding spacial derivatives. It is convenient to introduce the relative position vector ξ ξ = x0 − x and the relative displacement vector η 1
All vector quantities throughout this work are in bold
16
R f δ
x'
x
Figura 1.1: Domain and horizon
η = u(x0 , t) − u(x, t) = u0 − u The quantity η+ξ is therefore the current relative position vector between two particles. As shown in figure 1.1, it is convenient to define a region inside the domain, delimited by the distance δ, such that: if kξk ≥ δ =⇒ f (η, ξ) = 0 ∀η
(1.2)
The quantity δ is called the horizon, as it is the region in which a particle “sees ”forces from the surrounding ones. A typical value for δ is three times the grid spacing ∆x; different values are to avoid. Whether an horizon comparable to the grid spacing is used, cracks follow unnatural paths (along rows or columns of the mesh), and a smaller timestep is needed; conversely, if a larger one is used, besides longer computation time, there will be an eccessive wave dispersion and the solid will behave like a fluid. The peridynamic method for continuum mechanics is nonlocal, as particles interact at a finite distance with each other, like in molecular dynamics. This is not the only nonlocal theory, but it is the only one eliminating spatial derivatives within its formulation. A particular attention must be put when chosing the peridynamic approach to use: this work is based upon the “bond approach ”, which has the following limitations (as stated in ref. [8]): • can not directly use a constitutive model from standard theory 17
ξ+η
u
x
u'
R
ξ
x'
Figura 1.2: Variables
• assumes bonds act indipendently of each other, implying only a poisson coefficient ν = 0.25 can be used • can have permanent deformation of bonds, but can’t enforce plastic incompressibility • can model discrete particles, but only using two-body potentials For getting over these contraindications, a new approach has been developed: peridynamics states. This field of research is however not part of this thesis, while could be an interesting evolution for the application of the numerical stability method.
1.1.1
Pairwise force function
The specific force term in equation 1.1 is a functional depending on a vectorvalued function f , called the pairwise force function (PFF). The PFF gives the force per unit volume squared on a particle at x due to a particle at x0 . In conventional continuum-mechanics theory, the functional is the divergence of the stress tensor. All material properties and its behaviour are defined by the form of the pairwise force function. The dependence of f on x−x0 implies homogenity of the material; the lack of a temporal term but the current one avoids the dependence on deformation history, also. 18
To satisfy Newton’s third law the PFF must respect f (−η, −ξ) = −f (η, ξ) ∀η, ξ
(1.3)
meaning that f (η, ξ) is an odd function of (η, ξ). Moreover, the conservation of angular momentum in the absence of external forces implies that f (η, ξ) must also satisfy (if it is not conserved, two particles at rest would move even in absence of external forces)2 : (η + ξ) × f (η, ξ) = 0 ∀η, ξ
(1.4)
The force between any two particles has therefore to be parallel to their current relative position η + ξ. It is then possible to rewrite the pairwise force function in this form: f (η, ξ) = F (η, ξ)(η + ξ) ∀η, ξ
(1.5)
where F (η, ξ) is a scalar function, and must be an even function of(η, ξ), as f (η, ξ) is odd. A material is said to be microelastic whether (necessary and sufficient condition) a scalar-valued function w(η, ξ) exists, such that f (η, ξ) =
∂w (η, ξ) ∂η
(1.6)
The micropotential is the potential energy density of a bond (and its unit is energy per volume squared. For a microelastic material the general form of the pairwise force function is f (η, ξ) = H(p, ξ)(η + ξ) p = |η + ξ|
(1.7)
1w (p, ξ) (1.8) pp The basic consequence is that the force between any two points depends only upon their current relative position; a pratical physical interpretation is that every node is connected to the other ones inside its horizon through a spring. This spring in the most general case is nonlinear, but the model can be linearized as well. The stretch s of the bond can be defined as H(p, ξ) =
s= 2
|η + ξ| − |ξ| y − |ξ| = |ξ| |ξ|
× is the cross product
19
(1.9)
force rupture
stretch rupture
Figura 1.3: linear elastic
where y = |η + ξ| is the current relative position vector norm. A number of material behaviors can be defined from the pairwise force definition, such as: • proportional material • viscoelastic material • viscoplastic material
20
A proportional material is characterized by a PFF like f ∝ cs
(1.10)
where c is the spring constant. The spring constant, as described in [5], is obtained from the classical mechanics bulk modulus K, and particles’ horizon δ: 18K (1.11) πδ 4 In this thesis a proportional material model has been implemented. A further development of this work could lead to a structured analysis of the best pairwise force function able to describe any known material. c=
1.1.2
Damage
The main ability of peridynamics is to catch crack paths without needing any information but shape, material and pairwise force function. Damage is accounted at bond level: after trepassing a certain level of stretch, it breakes i.e. does not offer interaction force anymore. Following again [5], the critical stretch for bond failure is r 5G0 (1.12) s0 = 9Eδ The material behavior implemented in this thesis produces brittle fracture results; other models have been successfully implemented in other works cited as reference.
1.1.3
Short range forces
The basic formulation of peridynamic theory for continuum mechanics considers only interactions through bonds: particles that are initially separated (because their distance exceeds δ) can not interact in any way during the following timesteps. This limitation allows to trat crack propagation due to loads, but neglects impact between different bodies. To encompass this interactions, short range forces are here defined: these forces do not depend on the distance between the particles (the computational nodes) in the reference configuration, but only in the current state. The direction is again along the distance between the particles, and the module is given by: 0 |y − y| −1 (1.13) fs = min 0, cs 2rs 21
(a) Microelastic
(b) Fragile microelastic
(c) Microplastic
22
where y = x + u and y0 = x0 + u0 are the current positions of x and x0 , cs and rs positive constants. Repulsion forces are nonzero when the particles are closer than 2rs , which is called the node radius. Usually it is assumed that cs = 15c.
1.1.4
Peridynamic states
These models allow only one elastic constant, c, but isotropic linear elastic materials are actually characterized by two constants. An elastic solid that involves only two-particle interactions (a “Cauchy crystal ”) always has a Poisson ratio of 0.25. The theory of peridynamic states is a refinement of the peridynamic theory that allows dependence of strain energy density on local volume change, and not only two-particle interactions. It can also be implemented to model fluids. The bond properties depend not only on the stretch in the bond itself, but also on the stretch of all bonds connected to its node. This permits more general constitutive models.
23
Capitolo 2 Numerical method 2.1
Direct numerical method
The nodes form a grid. The method is mesh free: there are no elements or other geometrical connections between the nodes. The discretized form of the equation of motion replaces the integral by a finite sum: X f (unp − uni , xp − xi )Vp + bni (2.1) ρu¨ni = p
being uni = u(xi , tn )
(2.2)
Vp is the volume of node p, which for a uniform rectangular grid is simply x . The sum is taken over all nodes p such that |xp −xi | < δ. The discretized equation of motion is easy to implement in its fully nonlinear form. However, for purposes of error analysis and numerical stability analysis, it is necessary to work with a linearized version. The discretized form of the linearized peridynamic model is X ρu¨ni = C(unp − uni )(xp − xi )Vp + bni (2.3) 3
p
An explicit central difference formula is used for acceleration un+1 − 2uni + uin−1 i (2.4) ∆t2 Through the use of standard error analysis, the spatial and temporal discretization errors can be shown to be O(∆x2 ) + O(∆t2 ). In addition, von Neumann stability analysis leads to a time step requirement of ρu¨ni =
24
s ∆t <
2ρ p Vp |C(xp − xi )|
P
(2.5)
The value of the horizon δ may depend on the physical nature of the application being modeled. For example, at the nanoscale, δ would be determined by the distance over which physical interaction between atoms or molecules occurs. However, in macroscale calculations, it can be chosen according to convenience, since, as shown previously, the parameters such as c that determine the bulk elastic properties of the material can be fitted to experimental data for any value of δ. In practice, for macroscale modeling, the value δ = 3δx usually works well. Values much smaller than this typically result in undesirable grid effects (cracks grow along the rows or columns of the grid). Values much larger than this may result in excessive wave dispersion and require very large computer run times. Boundary conditions are somewhat different in the peridynamic method than in the classical theory. The variational formulation of the peridynamic equations does not lead to natural boundary conditions, which in the classical theory are traction boundary conditions. Instead, forces at the surface of a body must be applied as body forces b within some layer under the surface. In practice, this layer usually has thickness. Similarly, displacement boundary conditions must be prescribed within a layer of finite thickness under the surface. Although the above analysis assumes a regular lattice of nodes, the method is not restricted to such structured grids. Through the use of irregular grids, complex geometries may be modeled, and the full advantages of a meshless method are realized.
2.2
Time-adaptive method
The time stability is a big issue in peridynamics: when dealing with large deformations, complex geometries and fast dynamical phenomena the time step requirement is someway weak. Moreover, for the sake of conservativity, it can happen to choose a time interval which is unconvenient from a computational point of view: the time required to define the dominium behavior in an interval can happen to be out of scale. To overcome these drawbacks, a new time integration has been implemented, which permits to evaluate the integration’s correctness. This is performed by evaluating an error which is related to the motion of the particles, i.e. the new position can not be too far from the previous one. This 25
way a situation where the nodes tend to “explode ”is avoided, and position updates are kinda “smooth ”. A fundamental aspect which is not going to be heavily analysed, but that is crucial for a correct cycle of simulations, is machine precision, usually called . It is defined as the smallest number such that the difference between 1 and 1 + is nonzero. If the code has to deal with too little numbers, numerical errors risk to appear: the value of number summations is spoiled, therefore new “numerical ”and not physical forces appear; this is the consequence of the little difference between nodes’ positions in particular, with direct effect on the resulting bond stretch. This issue has to be solved by carefully checking every value to be far from the machine precision, and bringing it to zero whether it is too close. Let’s get back to the adaptive method. This approach comes from a solution to continuum thermo-chemical multifield problems developed by professor Zohdi of the University of California, Berkeley. It has been applied to peridynamics in a mixed explicit-implicit way: the integration through the timesteps is explicit, but at every time step an implicit cycle is performed to check whether that time interval ∆t is the best to achieve both stability and computational efficiency. The implicit iterations at every time step are called “Global fixed-point iterations ”. In the following, we recall that i indicates the nodes (total amount is N ), and K is the internal step counter. Before listing the algorithm, it is convenient to define some quantities to be used: PN n+1,K − ui n+1,K−1 || def i=1 ||ui (2.6) ωK = PN n+1,K − u n || ||u i i i=1 is the normalized error, and also: def
ZK =
ωK T OL
1 T OL pKd ω0 1 ωK pK ω0
(2.7)
def ΦK =
26
(2.8)
The global fixed-point iteration is then performed as: • for all the N nodes, – compute the new position as ui n+1,K ≈ uni
∆t2 f (ui n+1,K−1 ) + ∆tu˙ni + mi
– compute the new (internal cycle) interaction forces (storing them in temporary variables) • compute the error measures defined above • if tolerance met (ZK ≤ 1) and (K < Kd ), then – increment time t = t + ∆t and start from the beginning – construct new time step ∆t = ΦK ∆t • if tolerance not met (ZK > 1) and K = Kd then – construct new time step ∆t = ΦK ∆t – restart from time t The method is therefore characterized by three parameters, T OL, Kd and p. Typical values are T OL = 0.001, Kd = 20 and p = 3 but, as will be shown in the simulation results chapter, different values can be chosen for peridynamics applications. To summarize, the method controls how much do the nodes (particles) move: if the difference between a step and the other satisfies the tolerance, the time step is increased to improve efficiency; on the converse, if the particles move outside tolerance the step is reduced to keep the calculation stable.
27
Capitolo 3 Computational code The main operational goal of this work is the development of a calculation code implementing peridynamics theory (on a bond based approach and including the time-stable algorithm). Like every analysis suite, the solution of a problem is obtained after three steps, usually performed by different softwares: pre-processing, solving and post-processing. The core of the thesis is the creation from scratch of a stable solver, able to process as much as different geometries possible while presenting handful results. Moreover, for creating an useful chain of tools, pre and post-processing open source tools have been tested and integrated; this chapter presents a detailed description of them as well.
3.1
Solver
The solver is the core of this thesis work. It has been implemented from scratch using the c++ programming language; the compiler is GNU g++, from GCC (GNU Compiler Collection). GCC is the standard compiler in GNU/Linux environments; Richard Stallman started to write it in 1985, while the first release was in 1987: the aim was to be the compiler for the GNU Project, to have a compiler that was free software. At the current development stage, a c++ GNU compiler is available for both Linux and Windows platform; this means that, if OS-specific libraries are not included, the code is portable. The code for this thesis has been written and compiled on a Linux platform (Ubuntu 9.04 and Mandriva 2010), but it can be ported on windows as well.
28
To perform algebraic operations, file input and handful graphical output, two libraries have been included: Armadillo and VTK. Armadillo Armadillo (currently at version 0.9.5) is an open source (LGPL licensed) linear algebra library written in c++ and compatible with Matlab and Octave also. It is primarily developed at NICTA (Australia’s ICT Research Centre of Excellence), with contributions from around the world. It has functions for Matrices, Vectors, Scalars, Decompositions, Inverses, input and output; it has been chosen because it is simple to implement, utilize and has no license problems as it can be statically linked to the executable (i.e. can be part of the executable itself, and does not need to be present in the executing system). VTK VTK is an open-source, freely available software system for 3D computer graphics, image processing and visualization. It is produced by Kitware, which produces the postprocessor Paraview as well. TK supports a wide variety of visualization algorithms including: scalar, vector, tensor, texture, and volumetric methods; and advanced modeling techniques such as: implicit modeling, polygon reduction, mesh smoothing, cutting, contouring, and Delaunay triangulation. VTK has an extensive information visualization framework, has a suite of 3D interaction widgets, supports parallel processing, and integrates with various databases on GUI toolkits such as Qt and Tk. VTK is cross-platform and runs on Linux, Windows, Mac and Unix platforms. In this work it is used for producing simply accessible output files, to be read inside VisIt for real-time processing and animations production. Computational algorithm The code is composed by a main file and a library, which collects functions and classes. This software makes use of object oriented programming: nodes and bonds are classes, the majority of functions are method of the classes themselves. A brief overview of the algorithm is following:
29
Start Read input parameters
Read input file
Create nodes
Create links
Peridynamic iterations no
Iterations performed? yes Stop Figura 3.1: Main program flowchart
30
Velocity Verlet step 1 and 2
Stable time step calculation
Bond force update
Nodal force computation
Velocity Verlet step 3
Damping correction for velocities
Broken links dump
Output solution at current time step Figura 3.2: Peridynamic iterations flowchart
31
3.1.1
Parallelization
A dynamical problem can involve a considerable amount of time-steps to find the desired solution; numerical stability is a constraint that has to be considered as well. Moreover, regarding peridynamics, a large amount of nodes and especially of links between them raises the calculation power needed. This is someway similar to molecular dynamics, whether the above mentioned “collocation method ”is used. The most simple plate fracture calculation with a medium quality “mesh ”(i.e. spacing between the nodes) can take up to some hours, depending on the material properties. For industrially valuable simulation the number of nodes rises up of some order of magnitude, therefore large clusters of processors must be used. This implies a modification into the code, to exploit all the processors at the same time; two solutions (libraries) are available to solve this issue: OpenMP and MPI. MPI MPI is the acronym for Message Passing Interface: it is basically a communication protocol that supports point-to-point and collective communication. It is the standard for high performance computing. The basic concept is that every processor is a separate unit, and communicates with the other ones through messages; the balancing of this information exchange is fundamental for computational efficiency, as transmission time is subtracted to the calculation one. For a meshless method, three approaches can be chosen (coming from molecular dynamics, which is computationally similar to the peridynamic method): • node decomposition • bond decomposition • spatial decomposition In node decomposition algorithms, every processor is assigned a determined number of nodes, and calculates position, velocities and forces of those nodes only; in bond decomposition, the processors “own ”a bond, instead of a node, and exchange information about the latters. Eventually, in spatial decomposition algorithms the domain is divided into a number of 3d boxes, and the communication is just about nodes in the boxes’ surroundings. 32
OpenMP OpenMP is the acronym for Open Multi-processing, and it is an application programming interface for shared memory applications. It does not involve communication between processors, as it is assumed that they share the memory on which they store the variables. Whenever a non-conceptually serial operation is performed, the processors access different zones of the memory itself and perform their calculations just on that limited domain; besides they act independently, the data are freely accessible to every processor, in every moment. An MPI-based software architecture (carefully coded) is the de-facto standard for high performance super computers; but during this work there was not enough power available to justify such an implementation. The hardware for all the computations was a simple dual-processor laptop for which the MPI implementation would not have caused a considerable decrease of running time; therefore, because of both its code-level simplicity and efficiency for low number of processors platforms, OpenMP has been chosen. The speedup of the program is limited by the number of serial operations; in particular, the output writing phase is the performance limiter. If the output file is rarely produced (which is a software parameter), it is possible to observe both the processors working at 100%, i.e. the speedup tends to 2, which is the maximum available for a dual processors hardware.
3.2
Pre-processing
As described before, the solver needs a structured input file including node locations, initial velocity, external forces and particle’s mass; in a first step, this was accomplished by a part of the solver itself for simple rectangular plate domains. As the purpose of the work widened to a general purpose simulation suite, this gross generator has been put aside by the integration of a number of tools able to treat any geometry, keeping open source software as driving philosophy. The pre-processing phase is described in figure 3.3: the geometry is created with Salome’s geometry module and then meshed (again inside Salome) with the mesh module; this is a critical step, as the regularity of the mesh itself is a fundamental parameter of peridynamics numerical stability (one of the new research topics in this field is indeed spacial discretization). The mesh is then exported in I-DEAS format (unv) to be processed by Impact: inside this powerful finite element software the mesh is “cleaned ”by
33
Operation
Software
Geometry Salom`e Meshing UNV file Element deletion Impact Constraints, external loads and initial velocities
Python script
Final parsing
Figura 3.3: Pre-processing flowchart
eliminating the elements, just keeping the nodes (the code is meshless); next step is to apply constraints, external loads and initial velocity to the nodes and exporting the created calculation domain. However, the Impact output file is not in the required format, therefore a last parsing phase (by python scripting) is necessary before launching a calculation on it. The following sections present a detailed description of the above mentioned tools.
3.2.1
Salome
Salome is an open source pre and post processing platform for numerical simulations, GNU LGPL licensed therefore freely downloadable from the website (www.salome-platform.org). The project was born from a joint venture of 9 different companies interested in a multi-purpose numerical platform: Open Cascade, CEA, EDF, EADS, Bureau Veritas, Principia, Cedrat, LIP6 and LEG. The software is being developed since years, as the current version (june 2010) is 5.1.4 and 34
the first publicly available release (vers. 2.2.8) is dated july 2004. It runs on many different operative systems, from the native GNU/Linux environment (compiled versions for Mandriva and Debian distributions) to a brand new implementation on MS Windows. It is based upon Open Cascade libraries, and its code infrastructure is mainly C++ and python based; the main feature is just python scripting interaction: this allows a powerful parametric design and batch programming for more complex software chains. The architecture is made up of many modules: • Kernel : distributed components management, study management, general services • GUI : graphical user interface • Geometry: create, edit, import/export CAD models • Mesh : mesh a CAD model using a standard meshing algorithm or any external mesher (NETGEN, GHS3D, GHS3D parallel, BLSURF, Hexotic) • Med : MED data files management • Post-processor : dedicated viewer to analyze the results of solver computations (scalar, vectorial) • Supervisor : an easy way to describe a computational schema involving multi-solver coupling. Note: this module is not supported since SALOME version 5.0. • YACS : alternative supervision module (replaces Supervisor module in SALOME series 5x)
3.2.2
Impact
Impact is a finite element program for dynamic phenomena such as collisions, forming operations, chassis movements, etc. It was designed with the purpose to be a free and very simple alternative to commercial packages, being modular and supporting object programming. The current version is 0.7, therefore it is not a completely reliable tool still, but its preprocessor is very simple and effective and perfectly suites the need of manipulating the “meshless mesh ”and applying loads, constraints and initial velocities using a graphical interface. 35
A perhaps strange choice is the language: despite it is not well considered for high-performance computing, the software is implemented in Java. This is quite unusual, as usually C++ and Fortran are used, but this option has been chosen for the sake of programming simplicity and portability of the code itself; moreover, there are some publications like [11] describing how java can also be a fast and reliable environment for finite element calculations. It is an explicit time stepping program, so it is best suited for dynamics rather than statics.
3.2.3
Python
Python is a high-level programming language created in 1991 by Guido van Russum; the current developer of the platform is the Python Software Foundation. As the other tools used for this thesis, it is open source software (Python license). Python can be used as an object-oriented language, or as a structured one; the most important feature is that it has almost infinite bindings with external libraries, and covers a wide range of different applications. In this work, it has simply been used as ascii parser, for the manipulation of the output from Impact.
3.3
Post-processing
As previously mentioned, the solver outputs a certain amount of files, each one regarding a precise iteration time step; the availability of a bunch of output file formats aims to perform different tasks on the data. In fact, one needs to check the correctness of a running calculation as well as to publish clear and almost self-explanatory results. The former goal is achieved with VTK polidata files for visualization into VisIt or Paraview, the latter with geo files, which is the geometry input for GMSH. Both programs are to be described in the following sections.
3.3.1
VisIt
VisIt is an open source parallel visualization tool for scientific data. Its development began in the summer of 2010, by the Department of Energy (DOE) Advanced Simulation and Computing Initiative (ASCI) to visualize and analyze the results of terascale simulations. It is a professional software not only because of the quality of outputted plot, but for the ability of working with parallel data and across networks: 36
this is a very powerful feature, as it is possible for research teams spread all over the world to perform a remote visualization without the need of slow conversions and communications of massive data. In this thesis work it has been user for real time result visualization, to check the correctness of both code and input parameters; depending on the development step, it is so fast and well-conceived that it can save up to 50% of the debugging time.
3.3.2
GMSH
GMSH is a CAE (computer assisted engineering) platform based upon a flexible library for code reusability; it is written in C++ and again open source licensed (GNU GPL - LGPL). Although it is a complete (but no-frills) engineering suite including a professional mesher, in this work it has just been used for visualization purposes: in fact, this is the only simple software able to export vector images, which guarantee the higher quality for typography.
37
Capitolo 4 State of the art In this chapter an overview of the works that have been carried on in the previous years; since its creation, peridynamics has been growing in its theoretical foundations, peridynamic-based software programs and practical applications.
4.1 4.1.1
Software LAMMPS
LAMMPS is an acronym for Large-scale Atomic/Molecular Massively Parallel Simulator. The current version of LAMMPS is 15 Jan 2010. It has potentials for soft materials (biomolecules, polymers) and solid-state materials (metals, semiconductors) and coarse-grained or mesoscopic systems. It can be used to model atoms or, more generically, as a parallel particle simulator at the atomic, meso, or continuum scale. The software runs on single processors or in parallel using message-passing techniques and a spatial-decomposition of the simulation domain. The code is designed to be easy to modify or extend with new functionality. The code is distributed as an open source code under the terms of the GPL. The current version can be downloaded here. Links are also included to older F90/F77 versions. The last major release is also available on SourceForge. LAMMPS is distributed by Sandia National Laboratories, a US Department of Energy laboratory. Funding for LAMMPS development has come primarily from DOE (OASCR, OBER, ASCI, LDRD, Genomes-to-Life) and is acknowledged here. The LAMMPS WWW site is hosted by Sandia Laboratories.
38
General features: • runs on a single processor or in parallel • distributed-memory message-passing parallelism (MPI) • spatial-decomposition of simulation domain for parallelism • open-source distribution • highly portable C++ • optional libraries used: MPI and single-processor FFT • easy to extend with new features and functionality • runs from an input script • syntax for defining and using variables and formulas • syntax for looping over runs and breaking out of loops • run one or multiple simulations simultaneously (in parallel) from one script Particle and model types: • atoms • coarse-grained particles (e.g. bead-spring polymers) • united-atom polymers or organic molecules • all-atom polymers, organic molecules, proteins, DNA • metals • granular materials • coarse-grained mesoscale models • extended spherical and ellipsoidal particles • point dipolar particles • rigid collections of particles • hybrid combinations of these 39
Force fields • pairwise potentials: Lennard-Jones, Buckingham, Morse, Yukawa, soft, class 2 (COMPASS), tabulated • charged pairwise potentials: Coulombic, point-dipole • manybody potentials: EAM, Finnis/Sinclair EAM, modified EAM (MEAM), Stillinger-Weber, Tersoff, AI-REBO, ReaxFF • coarse-grained potentials: DPD, GayBerne, REsquared, colloidal, DLVO • mesoscopic potentials: granular, Peridynamics • bond potentials: harmonic, FENE, Morse, nonlinear, class 2, quartic (breakable) • angle potentials: harmonic, CHARMM, cosine, cosine/squared, class 2 (COMPASS) • dihedral potentials: harmonic, CHARMM, multi-harmonic, helix, class 2 (COMPASS), OPLS • improper potentials: harmonic, cvff, class 2 (COMPASS) • polymer potentials: all-atom, united-atom, bead-spring, breakable • water potentials: TIP3P, TIP4P, SPC • implicit solvent potentials: hydrodynamic lubrication, Debye • long-range Coulombics and dispersion: Ewald, PPPM (similar to particlemesh Ewald), Ewald/N for long-range Lennard-Jones • force-field compatibility with common CHARMM, AMBER, OPLS, GROMACS options • handful of GPU-enabled pair styles hybrid potentials: multiple pair, bond, angle, dihedral, improper potentials can be used in one simulation overlaid potentials: superposition of multiple pair potentials Atom creation: • read in atom coords from files 40
• create atoms on one or more lattices (e.g. grain boundaries) • delete geometric or logical groups of atoms (e.g. voids) • replicate existing atoms multiple times • displace atoms Ensembles, constraints, and boundary conditions: • 2d or 3d systems • orthogonal or non-orthogonal (triclinic symmetry) simulation domains • constant NVE, NVT, NPT, NPH, Parinello/Rahman integrators • thermostatting options for groups and geometric regions of atoms • pressure control via Nose/Hoover or Berendsen barostatting in 1 to 3 dimensions • simulation box deformation (tensile and shear) • harmonic (umbrella) constraint forces • rigid body constraints • SHAKE bond and angle constraints • bond breaking, formation, swapping • walls of various kinds • non-equilibrium molecular dynamics (NEMD) • variety of additional boundary conditions and constraints Integrators: • velocity-Verlet integrator • Brownian dynamics • rigid body integration • energy minimization via conjugate gradient or steepest descent relaxation
41
• rRESPA hierarchical timestepping Output: • log file of thermodynamic info • text dump files of atom coords, velocities, other per-atom quantities • binary restart files • per-atom quantities (energy, stress, centro-symmetry parameter, CNA, etc) • user-defined system-wide (log file) or per-atom (dump file) calculations • spatial and time averaging of per-atom quantities • time averaging of system-wide quantities • atom snapshots in native, XYZ, XTC, DCD, CFG formats Pre- and post-processing : • Various pre- and post-processing serial tools are packaged with LAMMPS • separate toolkit called Pizza.py which provides tools for doing setup, analysis, plotting, and visualization for LAMMPS simulations. Pizza.py is written in Python and is available for download from the Pizza.py WWW site. Specialized features (these are LAMMPS capabilities which you may not think of as typical molecular dynamics options): • real-time visualization and interactive MD • atom-to-continuum coupling with finite elements • coupled rigid body integration via the POEMS library • parallel tempering • parallel replica dynamics • Direct Simulation Monte Carlo for low-density fluids • Peridynamics mesoscale modeling • targeted and steered molecular dynamics • two-temperature electron model 42
Implementation We now detail how Eq. 1.1 was implemented in LAMMPS. First, a set of consistent units suitable for macroscopic simulations was needed. LAMMPS allows the user to choose units convenient for their particular simulations, each of which is implemented as a handful of conversion factors. These are used when, for ex- ample, kinetic energy 1/2mv 2 is computed with velocities in Angstroms/fs and the result should be in Kcal/mole (for an atom- istic simulation). For peridynamic simulations, an “SI” option was added where energy = Joules, distance = meters, time = seconds, etc. LAMMPS operates in parallel in a spatial-decomposition mode, where each processor owns a sub-domain of the overall simulation box and the particles within the subdomain. To compute pairwise or bond forces, a processor communicates with its neighboring processors via distributed-memory message passing (MPI) to acquire information about nearby “ghost ”particles owned by those processors. To improve the efficiency of the computation of pairwise forces, LAMMPS uses Verlet neighbor lists that are recomputed every few timesteps via binning particles on a regu- lar grid and searching nearby bins to find neighbors. These are also the timesteps on which particles migrate to new processors as needed. For peridynamic simulations, Verlet neighbor lists are used only when computing the family of short-range interacting particles, FiS . Each PD particle stores 4 quantities in addition to the usual coordinates, velocities, forces, etc. These are the critical strain s0 , the particle volume V , the particle density ρ , and the initial position of the particle x0 . V is precomputed for each particle based on the initial problem geometry. These 4 quantities migrate with particles as they move from processor to processor. Additionally, the s0 value for each ghost particle is communicated every timestep, since it is a dynamic quantity and the breakage criterion for bond i , j depends on the s0 value of both particles. Similarly, the V value for ghost particles is needed to scale the bond force between particles i , j. Since V is a static quantity, it is only communicated on timesteps when reneighboring is done. The short-range term is conceptually identical to standard short-range pair potential computations within an molecular dynamics code. In LAMMPS a “half ”neighbor list is used to efficiently find neighboring particles within a distance d S on a given timestep. By “half ”we mean that any interacting pair i , j is only stored once in the list, either by particle i or by particle j. The forces on particles i , j are computed for each pair in the list. At the end of the force computation, forces on ghost particles are communicated back to the owning processor. 43
Computation of the bond term requires knowledge of which bonds are already broken. Thus each particle stores a list of its bond partners, denoted by Fi, and flags them as they break. For each partner, the initial bond distance xj − xi is also stored, since it is used to calculate bond strength. The bond family of particles Fi for each particle is computed only once, on the first timestep of the simulation, based on the initial undeformed state of the material. The union of Fi over all particles is effectively a “full ”neighbor list where the geometric neighbors of particle i within a cutoff distance δ are stored. By “full ”we mean the interacting pair i , j is stored twice, once by particle i and once by particle j. With this information, the bond term can be computed by looping over the particles in Fi for each particle i. The bond is skipped if it previously broken. If both particles i , j are owned by the processor, the bond is also skipped if i ¿ j, taking advantage of Newton’s 3rd law. For active bonds, the bond is flagged as “broken ”if the bond strain s(t, η, ξ) exceeds the current strain criterion s0 for either particle i or j. While bond forces are computed, a new strain criterion s0 is also calculated, which will be used to break bonds on the next timestep. Note that the min operation implies a loop over all bonds of particle i. As each bond is calculated, the contribution to the new s0 of both the i and j particles is accumulated. Thus at the end of the bond loop, each particle has a new s0 , valid for the next timestep. Since the bond partner list stores all bonds for each particle owned by a processor, no extra communication is necessary to generate s0 , e.g., due to bonds with ghost particles. The bond partner list is another particle property that must migrate with particles as they move to new processors. This is done by packing and unpacking the Fi set of bond neighbors and distances into a message-passing buffer each time a particle migrates. During this operation, broken bonds are pruned from the list, so that a minimal amount of information is communicated. Once the short-range and bond terms have been computed, a final scaling by the volume V i of each particle is applied. This results in an effective force on each particle that can be used by a standard MD time integrator (velocity Verlet in the case of LAMMPS) in the usual way to update particle velocities and coordinates. For PD models we use a constant NVE integrator. Thermostatting is not used, since temperature is an ill-defined quantity for macroscopic PD particles. The “pressure ”due to PD interactions can be computed via the virial in the usual MD manner, except that the kinetic energy term contributing to the pressure is ignored. Alternatively, a precise notion of PD stress has been for- mulated and may be used.
44
4.1.2
EMU
Emu is the first code based on the peridynamic theory of solid mechanics. It is currently in the research and development stage and is being licensed to a few “alpha” users, and there has never been any information about a public release of the code, both from a source and binary side. It is therefore impossible to describe the software itself; by the way, a qualitative comparison can be performed against its published results.
45
4.2
Simulations
Peridynamics has been previously applied to a bunch of problems: some of them, taken from [8], are below reported:
4.2.1
Perforation of ductile targets
4.2.2
Perforation of reinforced concrete
46
4.2.3
Perforation of a fabric
4.2.4
Ball explosion
47
4.2.5
Time rate influence
48
Capitolo 5 Simulations 5.1
Membrane
The first performed simulation is an obscillating membrane: every one of its 4 sides is constrained, and the center of the dominium is assigned an initial velocity to produce oscillations. The first figures describe the behaviour of the numerical method, showing the time step values and the number of times the cycle has been repeated internally before proceeding. The last ones describe the actual response of the structure, with the oscillations and damping that brings the dominium back to its rest position.
Figura 5.1: Model with constraints and initial conditions
49
(a)
(b)
Figura 5.2: ∆x = 0.5
50
(a)
(b)
Figura 5.3: ∆x = 0.5
51
(a)
(b)
Figura 5.4: ∆x = 0.5
52
(a)
(b)
Figura 5.5: ∆x = 0.5
53
(a)
(b)
Figura 5.6: ∆x = 2
54
(a)
(b)
Figura 5.7: ∆x = 2
55
(a)
(b)
Figura 5.8: ∆x = 2
56
(a)
(b)
Figura 5.9: ∆x = 2
57
(a)
(b)
(c)
58
(d)
(e)
(f)
59
(g)
(h)
(i)
60
5.2
Plate with hole
A plate with hole is launched to show how does peridynamic deal with fracture; the response is the one expected, confirming the value of the method. There is no simmetry in the response, as we would expect, just because the mesher did not place the nodes symmetrically.
Figura 5.10: Model with constraints and initial conditions
61
(a)
(b)
(c)
62
(d)
(e)
(f)
63
(g)
(h)
(i)
64
(j)
(k)
(l)
65
5.3
Traction test
A traction test is here described. The method is again demonstrated to be stable, but we especially get the fundamental result of the crack path, which is perfectly coherent with brittle fracture theory. Peridynamics is then a proven method for fracture mechanics.
Figura 5.11: Model with constraints and initial conditions
66
(a)
(b)
(c)
(d)
67
(e)
(f)
(g)
(h)
68
5.4
Impact
An overview of peridynamics impact calculations capabilities is here presented. This test case is heavily different from the previous ones, as it involves repulsion forces. This implies not only the development of a specific force model, but especially a new p ¨roximity check¨. When dealing with contact, the starting positions are not taken into account anymore, and every particle could interact with each other at every time step. This must be checked, and short range forces should be created throughout all the calculation process. Besides the correctness of the concept, this is computationally unacceptable: the new forces are calculated whenever the previously defined error ωK passes a defined threshold, giving the information that particles are moving enough to create the opportunity of a contact. This is a good power-saving idea, but the check is computationally expensive even in this latter case, so a simple two dimensional calculation has been performed.
69
(i)
(j)
(k)
70
(l)
(m)
(n)
71
(o)
(p)
(q)
72
Capitolo 6 Conclusions The brand new peridynamic method for continuum mechanics has been theoretically analyzed and implemented from scratch in a c++ code. This bond-based software had to lead with stability problems, and computational inefficiencies due to the linearized model which was used to determine the integration time step size. A modified numerical method has been introduced to fix this issue: it carries on the explicit integration, but for each time step it gives a δt which is the best to avoid instability and to have computational efficiency as well. This is performed by an implicit integration within the explicit time-step. The stable code has subsequently been used to test peridynamic capabilities of dealing with dynamic phenomena like cracks and impacts: although only simple calculations have been carried on, due to the low available computational power, qualitative results confirm that this approach is both powerful and effective.
73
Bibliografia [1] Stewart A. Silling, Reformulation of elasticity theory for discontinuities and long-range forces, Journal of the Mechanics and Physics of Solids 48 (2000) 175-h209 [2] Michael L. Parks, Richard B. Lehoucq , Steven J. Plimptonb, Stewart A. Silling Implementing peridynamics within a molecular dynamics code, Computer Physics Communications 179 (2008) 777–783 [3] R.B. Lehoucq, S.A. Silling, Force flux and the peridynamic stress tensor, Journal of the Mechanics and Physics of Solids 56 (2008) 1566–1577 [4] Paul Demmie and S.A. Silling, An approach to modeling extreme loading of structures using peridynamics, Journal of mechanics of materials and structures, Vol.2 No.10, 2007 [5] S.A. Silling and E. Askari, A meshfree method based on the peridynamic model of solid mechanics, Comput. Struct. 83 (2005) 1526-1535 [6] Stewart Silling, Peridynamic states, Minisymposium on multiscale material modeling with peridynamics (presentation) [7] Markus Zimmermann, A continuum theory with long-range forces for solids, phd thesis at MIT, Department of Mechanical engineering (2005) [8] Stewart Silling, Peridynamic Modeling of the Failure of Heterogeneous Solids ARO Workshop on Analysis and Design of New Engineered Materials and Systems with Applications [9] T.I. Zohdi, A computational framework for network modeling of fibrous biological tissue deformation and rupture, Comput. Methods Appl. Mech. Engrg. 196 (2007) 2972–2980
74
[10] C. Geuzaine and J.-F. Remacle, Gmsh: a three-dimensional finite element mesh generator with built-in pre- and post-processing facilities, International Journal for Numerical Methods in Engineering, Volume 79, Issue 11, pages 1309-1331, 2009. [11] G.P.Nikishkov, Programming Finite Elements in Java, Springer, 2010, 402 pp. [12] S. J. Plimpton, Fast Parallel Algorithms for Short-Range Molecular Dynamics, J Comp Phys, 117, 1-19 (1995)
75
