Appunti di Modellistica Numerica Avanzata nella Progettazione Meccanica

Modellistica numerica

avanzata nella progettazione

meccanica

Algoritmi di time-stepping

Un problema dipendente dal tempo ha bisogno di uno schema che converta equazioni differenziali

in equazioni algebriche. È necessario quindi discretizzare una soluzione del tipo u(x, t) che quindi

dipende dallo spazio e dal tempo. Una possibile discretizzazione della soluzione è la seguente dove

per ora consideriamo un solo elemento:

L’approssimazione qui è realizzata con i polinomi interpolanti a cui abbiamo aggiunto la

dipendenza dal tempo (n nella sommatoria è il numero di nodi nell’elemento).

Può essere fatta anche la seguente formulazione:

Con questa formulazione si tengono le funzioni di forma che avevamo e si associa agli spostamenti

nodali la variabilità nel tempo. Questa formulazione è più semplice ma non è detto che una

combinazione lineare di funzioni rispettivamente nello spazio e nel tempo sia sempre adatta per

descrivere la variabilità della variabile primaria nella realtà.

Perciò una formulazione generale nello spazio e nel tempo della realtà tradotta nel linguaggio ad

elementi finiti è data dalla prima formulazione, ma per ∆t tendenti a 0, cioè per time-step molto

piccoli nella realtà, anche se la separazione delle variabili spaziali e temporali è scorretta l’errore

introdotto è tollerabile. Di conseguenza nei codici di time-stepping è utilizzata la seconda

approssimazione.

Gli algoritmi di time-stepping sono dunque degli algoritmi che si usano per le analisi di dinamica

veloce.

Discretizzazione del tempo

Per la discretizzazione del tempo è necessario:

• definire un passo temporale ∆t;

• trovare delle relazioni di ricorrenza che legano la soluzione al tempo n con la soluzione al

tempo n + 1.

Prendiamo la solita equazione differenziale che abbiamo visto fin dall’inizio del precedente corso

(un po’ semplificata) e abbiamo adesso un carico applicato nel tempo e anche le nostre incognite

varieranno nel tempo; si parte inoltre da una condizione iniziale, cioè al tempo t = 0, che è nota.

Si possono trovare le equazioni di ricorrenza con vari metodi e noi vedremo la famiglia-α, in cui α

non è altro che una sorta di peso che definisce il comportamento dell’algoritmo. Con questo

parametro si introduce un’equazione secondo cui la derivata prima della variabile è uguale ad una

somma pesata della derivata al tempo s e della derivata al tempo s + 1.

Abbiamo una somma pesata perché α è

compreso tra 0 e 1 e se aumenta α

diminuisce 1 – α e quindi viene dato più peso

alla derivata futura rispetto alla derivata al

tempo attuale. Nel caso in cui abbiamo α = 0

e α = 1 abbiamo rispettivamente le formule

di forward difference e backward difference:

−

+1

̇ ~

−

−1

̇ ~

−1

L’equazione vista viene riscritta in modo più sintetico, dove il primo membro dell’equazione viene

̇

scritto come :

+

In questo modo abbiamo già trovato un qualcosa che esprime la soluzione al tempo s + 1 come

una somma della soluzione al tempo s e di ∆t per una derivata della soluzione.

Riscriviamo l’equazione differenziale di prima all’istante di tempo s e all’istante di tempo s + 1 e

poi sostituiamo queste due scritture nella definizione della derivata vista prima:

Siamo arrivati a una scrittura in cui la soluzione al tempo s + 1 dipende dalla soluzione al tempo s

moltiplicata per un coefficiente dipendente da a e b, che sono informazioni specifiche del

problema, e da α e ∆t, più un termine in cui compaiono di nuovo il passo temporale, a e b, α e in

cui ci sono anche le informazioni sui carichi ai tempi s e s + 1. Tutte le informazioni a secondo

membro sono note perché conosciamo la soluzione al tempo s (quando si parte con l’analisi

questa sarà uguale alla soluzione iniziale), α che è definito prima di partire con l’analisi, ∆t, a e b, e

i carichi.

Bisogna tenere presente che le soluzioni u e u sono soluzioni nodali e poi se vogliamo la

s s+1

soluzione del modello ad un determinato istante di tempo si interpola con le funzioni di forma.

Stabilità

Riscriviamo innanzitutto l’equazione appena vista in termini più sintetici:

Se nel metodo di Newton-Raphson avevamo il controllo sul residuo, cioè finché non si arrivava a

convergenza non si incrementava il carico, in questo metodo non c’è un controllo, per cui si passa

da un istante temporale al successivo sapendo che la discretizzazione nel tempo porta a degli

errori perché non ci interessa il contenimento dell’errore al di sotto di una soglia stabilita

dall’utente, ma si vuole solo verificare che l’errore non aumenti nel tempo in modo esplosivo.

Vogliamo che l’errore sia stabile quindi anche se non si può sapere quanto è, e perché ciò avvenga

il parametro A, detto fattore di amplificazione dell’errore, deve essere minore di 1 in valore

assoluto.

Se la condizione non è rispettata l’algoritmo non è in grado di darci una soluzione. Questa

condizione dipende dalla fisica del problema perché compaiono a e b, dal time-step, e dal tipo di

algoritmo perché compare α.

Se la condizione è soddisfatta per ogni time-step allora l’algoritmo è detto incondizionatamente

stabile o più semplicemente stabile, mentre se è soddisfatta solo per alcuni time-step è detto

condizionatamente stabile.

La scelta del ∆t influenza l’accuratezza al di là della discretizzazione spaziale (la mesh) che

influenza la soluzione ad un dato istante di tempo, e inoltre la grandezza del passo con cui ci

muoviamo influenza anche la soluzione e come si accumula l’errore.

Possiamo tracciare il seguente diagramma:

Sulle ordinate si riporta A e sulle ascisse si riporta λ∆t, dove ∆t è il passo temporale con cui si

risolve e λ è un autovalore del sistema, cioè al suo interno ci sono a e b, cioè le informazioni

relative al problema. Possiamo tracciare varie curve al variare di α e vediamo che le curve con α =

1 e α = 0.878 stanno nel semipiano con A positivo, mentre le altre curve ad un certo punto vanno

al di sotto della linea dello 0 ma non intersecano mai la retta A = -1. La curva con α = 0 invece

interseca la retta A = -1 per un valore di ascissa 2. Fra α = 1/2 e α = 0 ci saranno delle curve che

prima o poi intersecheranno la retta A = -1. La retta con A = -1 definisce il semipiano con la

condizione di instabilità, cioè per A < -1 si sta violando la condizione di stabilità. Da un punto di

vista numerico c’è una differenza se l’algoritmo lavora con A tra -1 e 0 o con A maggiore di 0: se A

è maggiore di 0 la soluzione non è affetta da errore, mentre se A è compreso tra -1 e 0 c’è un certo

errore sovrapposto alla soluzione. Questo è il motivo per cui ad esempio nelle analisi di crash

vengono usati dei filtri passa basso sulla soluzione numerica perché certe volte questo errore crea

problemi (si lavora infatti nella zona con A compreso tra -1 e 0 in questo tipo di analisi).

Esempio: raffreddamento di una barra

Abbiamo una barra a temperatura iniziale uniforme soggetta ad una temperatura di 0° alle

estremità. Dovremmo avere un decadimento esponenziale della temperatura; la soluzione esatta

è rappresentata dalla curva exact tracciata con linea continua.

Queste soluzioni sono state ottenute con un ∆t di 1.5 s e con i vari algoritmi enunciati prima

(grafico a sinistra).

Il metodo della backward difference è incondizionatamente stabile come quello di Galerkin e di

Crank-Nicolson, e essi non vanno mai sotto lo 0, perciò non entrano mai in una zona in cui sono

affetti da rumore. Questi tre metodi riproducono con una certa approssimazione l’andamento del

raffreddamento, mentre il metodo di Eulero converge a temperatura 0° ma oscillando e andando

addirittura al di sotto della temperatura finale. Se rifacciamo lo stesso calcolo aumentando il passo

temporale fino a 2.5 s (grafico a destra) vediamo che i metodi di backward difference e di Galerkin

hanno lo stesso comportamento, anche se arriviamo alla soluzione con due passi di calcolo anzi

che tre (prima ci arrivavamo a 4.5 s e adesso a 5 s), mentre il metodo di Crank-Nicolson è sempre

incondizionatamente stabile ma è sceso al di sotto della linea dello 0, quindi presenta un

comportamento non fisico (c’è del rumore) secondo cui la barra si raffredda eccessivamente e poi

si riscalda di nuovo fino alla temperatura finale. Il metodo di Eulero invece oscilla con un’ampiezza

non più smorzata ma che aumenta, quindi l’algoritmo è andato a lavorare nella zona di instabilità.

Il rumore nel metodo di Crank-Nicolson potrebbe essere eliminato con un filtraggio mentre sul

metodo di Eulero non ci si può fare nulla. RADIOSS o LsDyna

usano lo schema di

risoluzione di Eulero

ma lo fanno perché

con α = 0 si hanno

delle semplificazioni

benefiche visto che i

calcoli sono molto

pesanti; lavorano nella

zona al limite

dell’instabilità

(cerchiata), e quindi la

soluzione sarà sempre

affetta da errore. Per ottimizzare ∆t, dato che c’è la possibilità di cambiarlo ad ogni passo

temporale, questo viene ricalcolato dopo ogni time-step.

Per i metodi condizionatamente stabili abbiamo una condizione sul time-step critico in

corrispondenza del quale ci troviamo al limite dell’instabilità, per cui a parità di problema fisico,

cioè a parità di λ, abbiamo un time-step critico e con la soluzione si deve stare al di sotto di questo,

pena il comportamento visto nel metodo di Eulero di prima.

λ è il più grande degli autovalori del modello definiti così:

Apparentemente sorge un problema perché per trovare il ∆t ad ogni passo temporale dovremmo

cri

fare un’analisi agli autovalori ad ogni passo temporale e individuare tra questi il massimo, il che

sarebbe dispendioso. In realtà ci viene in aiuto il teorema di Irons e Treharne che afferma che gli

autovalori del problema sono limitati superiormente e inferiormente dagli autovalori degli

elementi: il massimo degli autovalori del sistema è minore o uguale del massimo degli autovalori

di un elemento e il minimo degli autovalori del sistema è maggiore o uguale del minimo degli

autovalori di un elemento.

Quindi possiamo adesso ragionare in termini cautelativi, cioè se il più grande degli autovalori del

sistema ci dà un certo time-step critico il più grande degli autovalori dell’elemento dà un time-step

critico più piccolo perché il massimo degli autovalori dell’elemento è maggiore o uguale del

massimo degli autovalori del sistema, quindi basta assicurare che il time-step sia più piccolo del

time-step critico calcolato col massimo degli autovalori degli elementi e automaticamente esso

sarà minore del time-step critico calcolato col massimo degli autovalori del sistema. Quindi non c’è

bisogno di fare un’analisi agli autovalori del sistema ma basta trovare gli autovalori degli elementi.

Nel caso di problemi differenziali del secondo ordine il time-step critico è proporzionale quindi alla

dimensione dell’elemento. Esso può variare durante l’analisi e dipende dalla dimensione

dell’elemento e dalla velocità di propagazione delle onde longitudinali nel materiale. Se

consideriamo la formula del time-step critico di prima (considerando α = 0) per un elemento rod

lineare vediamo che il time-step critico dipende da h, che rappresenta la dimensione

dell’elemento, e da c, che rappresenta la velocità di propagazione delle onde longitudinali nel

materiale (a seconda che la formulazione della matrice di massa sia lumped o consistent avremo

due formule diverse): 2 ℎ ℎ

≤ = ≤

ℎ

√3

max

Per un elemento rod quadratico con matrice di massa lumped abbiamo

ℎ

≤ √6

Per un elemento beam lineare abbiamo 1

−

2

̂ 2

ℎ ℎ ℎ

( ) ( )

≤ min { , [1 + ] }

2

Per elementi quadrangolari abbiamo h e h che rappresentano le dimensioni dei lati e c che

1 2 d

rappresenta sempre la velocità di propagazione delle onde longitudinali nel materiale.

2 4

4

= ∑ ∑

2

=1 =1

( ) ( ) ( )

− − − ( − )

2 4 3 1 4 2 1 3

[] = [ ]

( ) ( ) ( )

− ( − ) − −

4 2 1 3 2 4 3 1

2

≤ ⁄

1 2

1

≤ 1

1 1 2

( )

+

2 2

ℎ ℎ

1 2

Dall’ultima formula si vede che se una delle due dimensioni della piastra è molto piccola sarà

quella che influenzerà di più il time-step.

Il solutore sa che all’interno del modello ci sono dei rod, o dei beam o ecc. e quindi per ogni

elemento da tali formule trova il time-step critico, e infine va a vedere quale è l’elemento che dà

una condizione più restrittiva sul passo temporale critico. Supponiamo adesso che c sia costante

per tutti gli elementi; h è direttamente proporzionale con il time-step critico, quindi si deduce che

per un modello di dinamica veloce conviene fare una mesh con elementi più grandi possibili

compatibilmente con il grado di precisione che vogliamo in modo che il time-step critico più

restrittivo sia abbastanza grande così da avere tempi di simulazione rapidi. In generale potremmo

pensare di fare una mesh biased, cioè più piccola in una zona e più grande in un’altra in modo da

risparmiare elementi (e quindi nodi) addensandoli in una zona di interesse, ma adesso l’effetto del

bias è negativo perché gli elementi più piccoli daranno una condizione restrittiva sul time-step

critico aumentando il numero di passi temporali. Nelle analisi di dinamica veloce si usano perciò

mesh di tipo uniforme.

Trattiamo adesso c che in un materiale è la velocità di propagazione del suono:

=√

c dipende dal modulo di Young e dalla densità. Se abbiamo un modello fatto da un solo materiale

l’unica variabile che fa variare il time-step critico è h, mentre se abbiamo un modello fatto da più

materiali il time-step critico è dato dalla combinazione della dimensione dell’elemento e della

velocità del suono all’interno di esso. Nella realtà se c < c dovremo imporre h < h in modo da

2 1 2 1

ottenere una relazione del tipo ℎ ℎ

1 2

=

1 2

in modo da realizzare una mesh con elementi di dimensione diversa a seconda del materiale pur

rispettando questa equivalenza per tutti gli elementi.

Il controllo su ∆t viene eseguito ad ogni istante temporale quindi h che usiamo per calcolare

inizialmente il time-step critico è la dimensione degli elementi che creiamo all’istante iniziale, ma

nei passi successivi gli elementi cambiano dimensione perché si sono deformati e h sarà quello

ottenuto dopo ogni deformazione ad ogni passo temporale. Quindi se ci fossero delle

deformazioni molto localizzate che tendono quasi ad annullare la dimensione dell’elemento la

simulazione andrebbe avanti con i passi temporali ma aumenterebbe il tempo necessario per

terminarla perché il time-step diminuirebbe sempre di più. Per ovviare a questo problema si può

fare la mesh più fitta nella zona in cui è presente una deformazione localizzata in modo da

ripartirla sui vari elementi.

Gli elementi solidi hanno un ulteriore problema: supponiamo che quello

in figura sia un elemento solido e che nella deformazione il nodo in nero

vada a finire al di sotto della base dell’elemento; questo porta ad

ottenere un valore negativo del volume e la simulazione si ferma.

Mass scaling

Fin dall’inizio potremmo fare una mesh con elementi piccoli oppure con elementi troppo piccoli in

una determinata zona e ciò sarebbe un problema: se a causa di questi elementi il time-step si

riduce ad esempio di un fattore 100 il tempo di tutta la simulazione si moltiplica per 100 perché

l’elemento che prevede la condizione più restrittiva sul time-step critico la impone su tutto il

modello. Per ovviare a questo problema si fa un mass scaling: le modifiche di h si subiscono

passivamente perché derivano dalla deformazione del modello, perciò si deve agire su qualche

altro parametro per limitare la riduzione di ∆t quando abbiamo uno o pochi elementi che limitano

tutto il modello. Si può agire su c ma non si può agire su E perché così cambiano le caratteristiche

del materiale e quindi le tensioni e le deformazioni, e dunque si agisce su ρ: si specifica a priori al

software di quanto può aumentare la densità di questi elementi. Aumentando ρ diminuisce c e

possiamo avere un h più piccolo a parità di time-step, oppure possiamo aumentare il time-step

critico fissato h. Tuttavia così facendo si aumenta l’inerzia del sistema e stiamo ragionando in

un’analisi di dinamica dove l’inerzia ha molto impatto, perciò non si può aumentare la densità di

due o tre volte altrimenti è come se avessimo una massa concentrata che modifica l’inerzia di

tutto il sistema.

Queste analisi è bene farle dove ci sono effetti dinamici molto forti per