Cap 4 Metodo SPH [PDF]

Zitiervorschau

4 METODO SPH

Il modello SPH (Smoothed Particle Hydrodynamics) è un metodo particellare meshfree basato sulla formulazione lagrangiana. Il suo utilizzo spazia diversi campi, dalle applicazioni astrofisiche ai problemi idrodinamici. Nell’ingegneria costiera i problemi analizzati riguardano la propagazione di onde in avvicinamento alla linea di costa. I principali vantaggi dell’SPH sono legati alla sua formulazione di tipo lagrangiano, grazie alla quale non abbiamo imposizioni restrittive sulla geometria del sistema, né su quanto ci possiamo allontanare dalle condizioni iniziali, essendo un metodo meshfree abbiamo la possibilità di studiare grandi deformazioni e di utilizzale materiali con legami costitutivi complessi. (Crespo, 2008).

4.1 Numerical Simulation. Formulazione essenziale dell’SPH L’idea di base sta nel trasportare gli aspetti di un problema fisico in un modello numerico, evitando di dover simulare l’esperienza reale con esperimenti complicati e costosi Il metodo SPH fu sviluppato per problemi idrodinamici nella forma di equazioni differenziali parziali (PDE) delle variabili di campo densità, velocità, energia, ecc. Ottenere soluzioni analitiche delle PDE è difficile, eccetto che per casi semplici. Si ha la necessità quindi di discretizzare prima il dominio del problema in cui sono definite le PDE. In seguito, è necessario un metodo per fornire un’approssimazione per i valori delle funzioni di campo e loro derivate, in qualsiasi punto. Alle PDE è quindi applicata una funzione di approssimazione per produrre un insieme di equazioni differenziali ordinarie (ODE) in forma discretizzata. Questo insieme di ODE discretizzate può essere risolto utilizzando una delle routine di integrazione standard del metodo convenzionale delle differenze finite. Nel metodo SPH, per raggiungere l’obiettivo di cui sopra, sono utilizzate le seguenti idee chiave:

1) Il dominio del problema viene rappresentato da un set di particelle arbitrariamente distribuite, se il dominio non è già in forma particellare. Non si necessita di nessuna connessione tra le particelle (meshfree); 2) Il metodo di rappresentazione integrale è usato per l’approssimazione delle funzioni di campo. Nell’SPH questo metodo prende il termine di approssimazione kernel (funzione di rappresentazione integrale); 3) L’approssimazione kernel è poi ulteriormente approssimata usando le particelle. In SPH è detta approssimazione particellare. Si opera sostituendo l’integrazione nella rappresentazione integrale della funzione di campo e sue derivate con sommatorie su tutti i valori corrispondenti alle particelle confinanti in un dominio locale chiamato dominio d’interazione. 4) L’approssimazione particellare è effettuata ad ogni istante, e quindi l’uso delle particelle dipende dalla sua distribuzione locale corrente (adattabilità);

5) L’approssimazione particellare viene eseguita per tutti i termini relativi alle funzioni di campo nelle PDE per produrre un insieme di equazioni differenziali in forma discretizzata rispetto al solo tempo (Lagrangiana); 6) Le ODE sono risolte utilizzando un algoritmo di integrazione esplicito in modo da avere una storia temporale di tutte le variabili di campo delle particelle (Dinamico). (GR Liu, 2003)

4.2 I Metodi Grid-Based tradizionale La simulazione attraverso i computer è diventata sempre più uno strumento importante per la risoluzione di problemi di ingegneria. Essa svolge un ruolo fondamentale, fornendo test ed esami necessari alla teoria, offrendo spunti per problemi fisici complessi, assistendo l’interpretazione e spesso la scoperta di nuovi fenomeni. I metodi numerici grid-based, come il metodo alle differenze finite (FDM), il metodo ai volumi finiti (FVM) ed il metodo agli elementi finiti (FEM), sono stati ampiamente usati in varie aree di fluidodinamica computazionale (CFD) e di meccanica dei solidi computazionale (CSM), per risolvere equazioni differenziali ed equazioni differenziali parziali (PDEs). Una caratteristica del metodo grid-based è quella di dividere il dominio continuo in piccoli subdomini discreti, grazie ad un processo denominato discretizzazione (meshing). Ogni nodo o punto della griglia è collegato agli altri in una maniera predefinita da una mappa topologica. Un sistema a griglia che consiste il nodi, celle o elementi, deve essere definito in modo da fornire una relazione tra i nodi prima del processo di approssiamzione per le equazioni differenziali. Le equazioni governative, basate opportunamente su una mesh predefinita, possono essere convertite in una serie di equazioni algebriche con incognite nodali per le variabili di campo. Nonostante il grande successo, questi metodi presentano dei problemi in molti aspetti che si ripercuotono nel poco efficace utilizzo per problemi complessi. Le maggiori difficoltà derivano dall’utilizzo delle mesh, che dovrebbero sempre assicurare che le condizioni di compatibilità numerica rispecchino quelle di compatibilità fisica per in continuo. L’utilizzo della griglia può portare ad avere difficoltà quando si affrontano problemi con superfici libere, contorni deformabili, interfacce mobili e grandi deformazioni. La descrizione lagrangiana, alla quale il nostro metodo fa riferimento, è tipicamente rappresentata dal metodo agli elementi finiti (FEM) ed è detta di tipo materiale. Infatti la mesh è fissa o assegnata al materiale per tutto il processo di calcolo. Dal momento che ogni nodo della griglia segue il percorso del materiale, il moto relativo dei nodi di collegamento può causare espansione, compressione e deformazione di una cella. Massa, momento ed energia sono trasportati con il moto delle celle della maglia. Poiché all’interno di ogni cella la massa resta invariata, non c’è flusso di massa da una cella all’altra ovvero attraverso i confini della mesh. Le caratteristiche di tale rappresentazione sono le seguenti:



 

 

Il termine convettivo nelle relative equazioni differenziali parziali non è presente, il codice è concettualmente semplice e dovrebbe essere più veloce data la mancanza degli oneri computazionali a risolvere i questi termini Essendo la griglia fissata sul materiale in movimento, l’intera storia temporale di tutte le variabili di campo in un punto materiale può essere facilmente tracciata Nel calcolo lagrangiano alcuni nodi della maglia possono essere posti lungo i confini e le interfacce dei materiali. Le condizioni al contorno e sulle superfici libere, i confini in movimento e le interfacce dei materiali sono imposte automaticamente, grazie al semplice movimento dei nodi della griglia Geometrie irregolari o complesse possono essere trattate convenientemente usando una maglia irregolare Il passo temporale, controllato dalla grandezza degli elementi più piccoli, può diventare troppo basso per essere efficiente e può anche portare al collasso del calcolo Una possibile opzione per migliorare il calcolo Lagrangiano è ri-zonare la mesh o rimeshare il dominio di calcolo. La rizonazione della mesh produce un sovrastato di una nuova, non distorta mesh sulla vecchia, distorta mesh, così che il calcolo successivo possa avvenire sulla nuova mesh invece che su quella vecchia (distorta). Le proprietà fisiche nelle nuove celle della mesh sono approssimate dalle celle della vecchia mesh attraverso il calcolo della massa, momento e trasporto di energia in una descrizione Euleriana. Le rizonazioni sono spesso utilizzate per simulare esplosioni, penetrazioni, impatti, turbolenze nel fluidi. La rizonazione nei metodi Lagrangiani comporta una grande spesa di tempo. Oltretutto con la rizonazione si perde anche la storia del materiale. Inoltre, i codici Lagrangiani sotto frequente ri-mesh assomigliano ad un codice Euleriano in senso generale. Pertanto, anche se ci sono alcuni vantaggi molto buoni nei metodi griglia-basati Lagrangiani, gli svantaggi possono provocare difficoltà numeriche durante la simulazione di eventi caratterizzati da grandi deformazioni. (GR Liu, 2003) La descrizione Euleriana è tipicamente rappresentata dal metodo ai volumi finiti (FVM) dove la mesh è fissa nello spazio ed il materiale può muoversi attraverso la griglia Una comparazione tra i due metodi è illustrata in tabella 2 (GR Liu, 2003), mentre in figura 2 sono rappresentate due griglie, una Lagrangiana, l’altra Euleriana

Tabella 2

Fig. 2 Mesh Lagrangiana

Fig. 2 Mesh Euleriana

4.3 Metodo Meshfree I metodi numerici grid-based presentano difficoltà in molti aspetti, che limitano l’applicazione di questi metodi in un gran numero di problemi complessi. Una delle più grandi limitazioni è la generazione della griglia, che non è sempre un processo semplice e può portare ad un lavoro oneroso sia in termini di tempo sia in termini di complessità matematica. I metodi meshfree agevolano la simulazione di problemi che richiedono la capacità di lavorare con grandi deformazioni, materiali avanzati, geometrie complesse, comportamento non lineare dei materiali , discontinuità e singolarità. In tabella 2 sono elencati, in ordine cronologico, i principali metodi meshfree

Tabella 2

4.3.2 Smoothed Particles Hydrodynamics I metodi particellari meshfree (MPM) trattano i sistemi come gli insiemi di particelle, che rappresentano un oggetto fisico o una porzione di un dominio. Per i problemi di Fluido Dinamica Computazionale (CFD) le variabili quali la massa, la quantità di moto, l’energia, la posizione, etc, sono calcolato per ogni particella. Alcuni esempi di questi metodi sono illustrati in tabella 3 (Tabella 2.4 in (GR Liu, 2003)).

Tabella 3

4.4 Smoothed Particle Hydrodinamics L’SPH è un autentico metodo particellare meshfree originariamente utilizzato per applicazioni nel continuo e può essere considerato come il più vecchio tra i moderni metodi meshfree. Fu inventato inizialmente per risolvere problemi astrofisici nello spazio tridimensionale, poiché il movimento collettivo delle particelle si avvicina molto al movimento di un liquido o di un gas e può essere modellato dalle equazioni dell’idrodinamica Newtoniana. (M.B. Liu, 2009) Nell’SPH lo stato del sistema è presentato da un set di particelle, le quali possiedono le proprietà dei materiali ed interagiscono tra loro in un campo controllato da una funzione peso (Smoothing Function). la pressione del fluido è calcolata dalla densità utilizzando un equazione di stato , l’accelerazione delle particelle è quindi ricavata dal gradiente della pressione e dalla densità. Per fluidi viscosi può essere incluso l’effetto della viscosità fisica nel calcolo dell’accelerazione delle particelle. Allo stesso modo di un metodo particellare Lagrangiano, l’SPH conserva esattamente il valore della massa. Vediamo quali sono i principali vantaggi di questo metodo rispetto un metodo numerico tradizionale grid-based 





SPH è un metodo particellare Lagrangiano, utilizza un algoritmo definito : Galilean Invariant. È possibile ottenete la storia temporale delle particelle materiali, l’avvezione ed il trasporto del sistema può essere quindi calcolato. Grazie alla corretta disposizione particellare in una specifica posizione ad un istanti iniziale prima dell’analisi, la superficie libera, le interfacce dei materiali ed i contorni mobili possono essere tracciate nel processo di simulazione indipendentemente dal movimento delle particelle. L’SPH è un metodo che non utilizza una mesh/griglia. Questo consente un trattamento corretto delle grandi deformazioni, dal momento che la connettività è generata come parte del calcolo e può cambiare nel tempo. Per questi motivi l’SPH viene utilizzato nei casi di fenomeni ad alta energia quali esplosioni, esplosioni sottomarine, impatti ad alta velocità, penetrazioni.







Nel metodo SPH , una particella rappresenta un volume finito alla scala del continuo; siamo abbastanza vicini al classico metodo della dinamica molecolare che utilizza una particella per rappresentare un atomo o una molecola, ed al metodo della dinamica delle particelle dissipative, che utilizza una particella per rappresentare un piccolo gruppo di molecole. L’SPH è adatto in quei casi dove l’oggetto in considerazione non è un continuo. Questo è vero specialmente nei casi di bio e nano ingegneria alla scala micro e nano, e di ingegneria astrofisica alla scala astronomica. L’SPH è più semplice per l’implementazione numerica ed è più naturale per lo sviluppo di modelli numerici tridimensionali dei metodi Grib-based Il primo algoritmo SPH era derivato dalla teoria della probabilità, questo algoritmo non conservava la quantità di moto ed il momento della quantità di moto. Con lo sviluppo del metodo e la sua grande applicazione molte caratteristiche interessanti sono state messe in luce, cosi come molte incongruenze sono state identificate, per questo sono state proposte varianti o modifiche. Per esempio Gingold e Monaghan introdussero un algoritmo capace di conservare sia la quantità di moto sia il momento della quantità di moto (Gingold RA, 1982). Hu e Adams presentarono anch’essi un algoritmo con le stesse capacità conservative per i flussi viscosi incompressibili. (Hu XY, 2006). Molti ricercatori hanno condotto studi su questo metodo riguardo l’accuratezza degli aspetti numerici, sulla stabilità, la convergenza e l’efficienza. Swegle identificò il problema della istabilità di tensione (Swegle JW, 1995). Morris notò il problema dell’inconsistenza delle particelle che era condotto con poca accuratezza nella soluzione dell’SPH (JP, 1996; Monaghan J. , 1994). Monaghan propose una simmetrizzazione della formulazione che ha mostrato ottimi risultati (Monaghan, Why particle methods work, 1982) (Monaghan J. , Smooth particle hydrodynamics, 1992). Randles e Libersky derivarono un formulazione normalizzata per il calcolo della densità approssimata ed una normalizzazione per la divergenza dei tensori degli sforzi (Randles PW, 1996). Chen et al. Proposero un metodo SPH corretto (CSPM) che sviluppava una accuratezza sia nei problemi del dominio sia intorno alle aree di confine. (Chen JK, 1999) (Chen JK B. J., 2000). (M.B. Liu, 2009)

4.5 Interpolazione Integrale L’SPH è basato sull’interpolazione integrale. Il principio fondamentale è quello di approssimare ogni funzione A(r) come (kernel appoximation):

4.1

dove r è il vettore posizione; W è la funzione peso o kernel; h è chiamata “smoothing length” e controlla il dominio di influenza Ω (figura 2)

Figura 2

Solitamente il valore di h deve essere più grande della distanza che separa inizialmente le particelle. L’approssimazione (4.1), in forma discreta, ci conduce all’approssimazione particellare:

4.2

dove la sommatoria è su tutte le particelle nella regione del supporto compatto della funzione kernel. La massa e la densità sono denominate rispettivamente e ; = W( − , ℎ) è la funzione peso o kernel. Uno dei vantaggi di questo approccio è che la derivata di una funzione viene calcolata analiticamente diversamente da come succede, ad esempio, nel metodo alle differenze finite. Le derivate in questo metodo possono essere ottenute grazie ad una derivazione ordinaria e non è necessario né un metodo di differenze finite né una mesh.

4.3

4.6 Il Kernel “Smoothing” La prestazione di un modello SPH dipende dalla scelta della funzione peso. Queste dovrebbero soddisfare condizioni restrittive quali positività, supporto compatto e normalizzazione; in più deve essere monotona decrescente con l’aumento della distanza tra le particelle a e deve comportarsi come una funzione delta quando la lunghezza h tende a zero:

Positività: nel dominio Ω 4.4 Supporto Compatto: fuori dal dominio Ω 4.5 Normalizzazione: 4.6

Comportamento funzione delta: 4.7

Comportamento monotono decrescente di

4.8

Il kernel dipende dalla lunghezza h e dalla distanza adimensionalizzata tra le particelle q=r/h, essendo r la distanza tra le particelle a e b. il parametro h controlla la grandezza dell’area intorno alla particella a all’interno della quale il contributo delle altre particelle non può essere trascurato. (Crespo, 2008)

4.7 Equazioni Governanti Le equazioni che stanno alla base della fluido dinamica si seguono se seguenti tre leggi di conservazione:   

Conservazione della massa Conservazione della quantità di moto Conservazione dell’energia Le equazioni del moto in SPH sono derivate da queste equazioni nella forma Lagrangiana

4.7.1 Conservazione della Quantità di Moto L’equazione di conservazione della quantità di moto in un campo continuo è:

2.9

Dove v è la velocità, P e ρ sono la pressione e la densità; g= (0.0.-9.82) ms² è l’accelerazione gravitazionale; Θ fa riferimento ai termini diffusivi. Essendoci diverse formulazioni riguardo i termini diffusivi, per ognuna di esse può esserci un approccio diverso nel descrivere le equazioni di conservazione della quantità di moto. Abbiamo tre possibili opzioni:   

Artificial viscosity Laminar viscosity Turbulence modeling (laminar viscosity + Sub Particle Scale (SPS) Turbolence)

Viscosità artificiale La viscosità artificiale proposta da (Monaghan J. , Smooth particle hydrodynamics, 1992) è stata molto spesso utilizzata a causa della sua semplicità. Nelle notazioni SPH la (2.9) può essere scritta come: 4.20

Il gradiente della pressione in forma simmetrica è espresso in notazione SPH come:

4.21

Πab indica la viscosità

4.22

con

dove

=

(

+

),

= (

+

), η²= 0.02h², α è un parametro libero che può

cambiare per ogni problema.

Viscosità laminare L’equazione di conservazione della quantità di moto con viscosità laminare è data da:

4.23

dove il termine dello sforzo laminare si semplifica:

2.24

Il termine v0 è la viscosità cinetica di un flusso laminare (0.893*20

m²)

Quindi in notazione SPH la (2.23) può essere riscritta come:

4.25

Viscosità laminare e SPS (Sub-Particle Scale) Per rappresentare adeguatamente la viscosità di un fluido ed un moto turbolento viene utilizzato il Large Eddy Simulation (LES), una modello matematico usato nella fluidodinamica computazionale (CFD) per lo studio di fenomeni turbolenti. L’equazione di conservazione della quantità di moto è:

4.26

dove il termine laminare può essere trattato seguendo l’equazione (2.24) e ̅ il tensore degli sforzi.

Le assunzioni di questa teoria (sotto le ipotesi di Boussinesq) sono spesso utilizzate per modellare il tensore degli sforzi SPS utilizzando una media alla Favre (per un fluido comprimibile):

4.27

dove è il tensore degli sforzi, = [min(Cs,Δl)*|S| (turbolence eddy viscosity), k la SPS energia cinetica della turbolenza, Cs la costante di Smagorinsky (0.22) CI =0.0066, Δl la distanza particella-particella, |S|= (2SijSij)^2/2. Quindi seguendo (Dalrymple R. A, 2006) l’equazione (2.26) può essere riscritta nella forma:

4.28 4.7.2 Equazione di Continuità Il fluido nella trattazione SOH standard è trattato come comprimibile, il che permette di utilizzare le equazioni di stato per determinare la pressione del fluido piuttosto che risolvere altre equazioni differenziali. In ogni modo la comprimibilità è regolata per rallentare la velocità del suono affinché il tempo di calcolo sia ragionabile. Le variazioni di densità sono calcolate a partire dalla:

4.29

invece di utilizzare una sommatoria pesata dei termini di massa (Monaghan,2992), poiché è nota per provocare un decremento della densità artificiale vicino alle interfacce fluide.

4.7.3 Equazioni di Stato Basandoci su quanto scritto da (Monaghan J. , 1994) e (Batchelor, 1974) la relazione tra pressione e densità è data dalla seguente espressione, nota come equazione di stato Tait:

4.30 si vede come una piccola oscillazione della densità produce una grande variazione della pressione. Questo fluido comprimibile ha un valore della velocità del suono, c, che è data dalla:

4.31 4.32

dove c0 è la velocità del suono ad una densità di riferimento (sulla superficie libera); la costante B=c0²ρ0²/γ pone un limite alla variazione della densità la scelta di B gioca un ruolo fondamentale dal momento che essa determina la velocità del suono. Usando un valore corrispondente al valore reale della velocità del suono in acqua, deve essere stabilito un passo in termini di tempo molto piccolo, che si basa sulla Courant-Fredrich-Levy condition. Monaghan mostrò che la velocità del suono deve essere rallentata artificialmente, tuttavia (Monaghan, Simulating free surface flows with SPH, 1994) suggerì che il minimo valore della velocità del suono dovesse essere circa dieci volte più grande che la massima velocità del fluido attesa.

4.7.4 Movimento delle Particelle Le particelle sono mosse utilizzando la variante XSPH (Monaghan, On the Problem of Penetration in Particle Methods, 1989)

4.33

Dove ε è una costante che varia tra 0 e 2, solitamente viene utilizzato ε=0.5.

4.7.5 Conservazione ell’Energia Durante la simulazione viene calcolata l’energia cinetica, potenziale e termica. L’energia termica associata ad ogni particella usando la viscosità artificiale è calcolata attraverso l’espressione data da (Monaghan, Smooth particle hydrodynamics, 1992)

4.34

L’energia totale del sistema è calcolata come la somma dell’energia cinetica, potenziale e termica. A titolo di esempio riportiamo in figura 3 le energie calcolate nel caso di collasso di una colonna d’acqua

Figura 3

Analizzando la figura, si osserva come l’energia cinematica per contorni stazionari sia uguale a zero. L’energia potenziale iniziale è zero. L’energia termica è calcolata basandosi sulla (2.24). Nell’ultima figura viene rappresentata l’energia totale del sistema; in accordo con (Monaghan, Smooth particle hydrodynamics, 1992), l’energia è conservata con un limite dello 0.5% in 400 step. In questa simulazione l’energia cresce

dello 0.3% in 500 step, quindi la variazione di energia è dentro i limiti proposti da Monaghan. (Crespo, 2008)

4.8 Scelta del Kernel Le approssimazioni kernel sono descritte ampiamente in (Monaghan, Smooth particle hydrodynamics, 1992), (GR Liu, 2003) e (Monaghan, Smoothed Particle Hydrodynamics., 2005). Nel nostro caso (DualSPHysics) è possibile scegliere tra i seguenti kernel:    

Gaussian Quadratic Cubic spline Quintic

Gaussiano 4.35

dove q=r/h, r è la distanza tra le particelle a e b e

= 2/(πh²) in 2D e 2/(

Figura 4

La figura 4 mostra il valore del kernel Gaussiano e le sue derivate.

ℎ ) in 3D

Qudratic

4.36

dove

= 2/(πh²) in 2D e 5/(4 ℎ ) in 3D

Figura 5

(GR Johnson, 1996) ha utilizzato questo kernel per simulare i problemi di impatti ad alta velocità. Questa funzione impedisce il raggruppamento delle particelle nei problemi di compressione.

Cubic Spline

4.37

dove

= 20/(7πh²) in 2D e 2/( ℎ ) in 3D

Figura 6

Questa funzione è stata la più usata in letteratura data la somiglianza con la funzione Gaussiana pur avendo un supporto compatto stretto. Uno dei vantaggi nell’utilizzo di questo kernel al posto di quello Gaussiano è che questo ha un supporto compatto e l’onere computazionale numerico è ridotto.

Quintic (Wendland, 1995)

4.38 dove

= 7/(4πh²) in 2D e 7/(8 ℎ ) in 3D

Figura 7

I risultati mostrano come il miglior compromesso tra l’accuratezza e onere computazionale del tempo è raggiunto dal kernel Wenland.

4.9 Consistenza e Stabilità In molti modelli numerici esiste un dilemma: consistenza o stabilità? Per settare un modello numerico come il metodo particellare, dobbiamo sceglierne una rispetto all’altra. L’SPH originale ha chiaramente preferito la stabilità (ed anche la flessibilità) sulla consistenza, che consente a questo metodo di lavorare bene con un gran numero di problemi complicati, senza però considerare troppo l’accuratezza. Questa sembra essere una scelta pratica per i problemi, pratici appunto, che riguardano l’ingegneria. Questa scelta dovrebbe essere considerata un vantaggio del metodo SPH. I tentativi di migliorare l’accuratezza possono essere di aiuto, a condizione però che la stabilità e l’efficienza non siano troppo compromesse. La questione è se si possano avere contemporaneamente accuratezza e stabilità. La risposta è si, a patto di cambiare il settaggio ed essere disposti a pagarne il prezzo. Il recente metodo proposto GSM è un tipico esempio in questa direzione (GR Liu Z. J., 2008) (XG Xu, 2009). Tuttavia questo metodo non è più un metodo particellare. Il GSM richiede precise stime degli integrali per scegliete attentamente il tipo di dominio “smoothing” da utilizzare, ed in un certo senso lavora come un FVM.