Universita degli Studi di Roma Tre

Dipartimento di Ingegneria Civile

Corso di Laurea Magistrale in Ingegneria Civile per la Protezione dai Rischi Naturali

Relazione di fine tirocinio

Metodi numerici di integrazione nel tempo di problemidinamici

Studente:

Gabriele PellicciaMatricola 430436

Tutor:

Prof. Giovanni Formica

Anno Accademico 2018-2019

Sommario

La seguente relazione descrive le attività svolte ai fini dello svolgimento della tesi di lau-rea, con particolare riferimento all’acquisizione di ulteriori conoscenze informatiche nelcampo della meccanica computazionale. Tali attività sono previste dall’art. 10, co. 5 let.d/e e considerate equivalenti al tirocinio; lo svolgimento delle attività di seguito illustratecorrisponde al riconoscimento di 6 CFU per un numero di ore complessivo pari a 150.Le attività si sono svolte nel periodo 04/03/2019 - 20/05/2019 presso il Laboratorio diModellazione e Simulazione (Dipartimento di Architettura dell’Università degli Studi diRoma Tre).L’obiettivo prevedeva l’approfondimento delle conoscenze relative all’applicazione del soft-ware MATLAB per lo studio delle proprietà di convergenza di solutori numerici di inte-grazione nel tempo.

Indice

1 Introduzione 3

2 analisi del problema dinamico 42.1 Equazione lineare del moto . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

3 Metodi di risoluzione ODE adottati 53.1 ODE45 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

3.1.1 Metodo di Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . 53.2 Metodo beta-Newmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

4 Confronto 10

1

Elenco delle figure

4.1 Spostamenti nel tempo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114.2 Velocità nel tempo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124.3 Fasi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134.4 Forza-spostamenti: non linearità . . . . . . . . . . . . . . . . . . . . . . . . 144.5 Spostamenti nel tempo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154.6 Spostamenti nel tempo tramite iterazione di Newton-Raphson . . . . . . . 16

2

Capitolo 1

Introduzione

Nel corso del presente Tirocinio è stato svolto uno studio di convergenza su problemidinamici semplici attraverso il confronto computazionale di diversi metodi: ODE45, beta-Newmark. Entrambi i metodi consistono nella risoluzione numerica di problemi con equa-zioni differenziali oridnarie al secondo ordine per poi essere ricondotte al primo. Il metodoODE45 che si basa sul metodo Runge Kutta e un linguaggio già implementato in MATLABdi default mentre il metodo Beta Newmark è stato implementato a mano.

3

Capitolo 2

analisi del problema dinamico

2.1 Equazione lineare del motoIl problema analizzato è del tipo

Mu+ Cu+Ku = f(t) (2.1)

dove M è la matrice delle masse, C la matrice di smorzamento viscoso, K la matricedelle rigidezze e f(t) il vettore delle forze esterne. Per u, u e u si intendono rispettiva-mente le evoluzioni temporali di accelerazioni, velocità e spostamenti. Si considerano lematrici M ,C e K simmetriche. I valori iniziali del problema consistono nel cercare unospostamento u = u(t) che soddisfi l’equazione del moto sopra indicata e i valori iziali:{

u(0) = u0

u(0) = v0

Dall’equazione omogenea associata è possibile, risolvendo un problema agli autovalori,individuare le N frequenze proprie, gli smorzamenti e le N forme modali (modi di vibrare)del sistema. Questo elenco costituisce ciò che si definisce come modello modale, cheidentifica una struttura dal punto di vista dinamico.Un ulteriore passo che si può fare è derivare il cosiddetto modello di risposta, che individuail comportamento dinamico della struttura in termini di risposta dinamica (spostamenti)per unità di forza esterna, nel dominio delle frequenze. Si ottiene così una funzione dettadi risposta in frequenza (FRF) hij(ω), ordinata in una matrice N x N, le cui entratecostituiscono la risposta dinamica nel grado di libertà i-esimo sotto una forzante unitarianel grado di libertà j-esimo.

4

Capitolo 3

Metodi di risoluzione ODE adottati

3.1 ODE45In matematica, un’equazione differenziale ordinaria (ODE) è un’equazione differenzialeche contiene una o più funzioni di una variabile indipendente e le derivate di tali funzioni.Il termine ordinario è usato in contrasto con il termine equazione differenziale parziale chepuò essere rispetto a più di una variabile indipendente.

Il software MATLAB dispone di diverse functions per risolvere sistemi di equazionidifferenziali che si differenziano per metodi impliciti e metodi espliciti, come in questocaso per quest’ultimo caso il comando ODE45, il più usato. Tutte le functions per ODEdi MATLAB sono basate sulla stessa sintassi, ne segue un esempio:

[t, y] = ODE45(fun, tspan, y0) (3.1)

dove

• fun è una stringa di caratteri che contiene il nome della funzione f(t, y) che definisceil sistema differenziale

• tspan è un vettore che contiene l’istante iniziale e lìistente iniziale:

tspan = [tf , t0] (3.2)

• y0 contiene le condizioni iniziali. È uno scalare se è presente un’unica equazionedifferenziale mentre è un vettore con m componenti se c’è da integrare un sistemacon m equazioni.

Il metodo esplicito ODE45 è basato su schemi Runge-Kutta.

3.1.1 Metodo di Runge-Kutta

In analisi numerica i metodi di Runge-Kutta sono un’importante famiglia di metodi itera-tivi impliciti ed espliciti per l’approssimazione delle soluzioni delle equazioni differenzialiordinarie. Queste tecniche furono sviluppate intorno al 1900 dai matematici tedeschi CarlRunge e Martin Wilhelm Kutta.

I metodi di Runge-Kutta (spesso abbreviati con "RK") sono una famiglia di metodiiterativi discreti utilizzati nell’approssimazione numerica di soluzioni di equazioni diffe-renziali ordinarie (ODE), e più specificatamente per problemi ai valori iniziali. Fannoparte della famiglia più generale di metodi discreti per le equazioni differenziali ordinarie,

5

ovvero di quella classe di metodi numerici che fornisce un’approssimazione della soluzionedi un’equazione differenziale (o più precisamente di un problema di Cauchy) in un insiemediscreto di punti.

Per trovare un un’approssimazione della funzione y(t) : R → Rd che verifichi il gene-rico problema di Cauchy:

{y′(t) = f(t, y(t))

y(t0) = y0

in un insieme discreto di punti in cui si considera il problema (solitamente nell’inter-vallo [t0, tf ], si considera un campionamento dell’intervallo ∆ in un insieme di punti{ti|i = 0...n}, dove ti = t0 + ih e h = (tf − t0)/n. Il metodo numerico fornisce allo-ra l’approssimazione dei valori y(tj), e per ottenere una ricostruzione abbastanza fedeledella funzione il numero n deve essere sufficientemente elevato.

Metodo esplicito

Dato il problema ai valori iniziali:

y = f(t, y) y(t0) = y0 (3.3)

dove i valori di t0 e y0 sono noti, si consideri un intervallo sufficientemente piccolo h > 0e si definiscano:

yn+1 = yn +h

6(k1 + 2 k2 + 2 k3 + k4)

tn+1 = tn + h(3.4)

per n = 1, 2, 3, .... In questo modo y(tn+1) è approssimato con yn+1, e yn+1 èdeterminato da yn più la media pesata di quattro incrementi k1, k2, k3, k4:

k1 = f(tn, yn)

k2 = f

(tn +

h

2, yn +

h

2k1

)k3 = f

(tn +

h

2, yn +

h

2k2

)k4 = f (tn + h, yn + h k3)

(3.5)

dove ogni incremento è il prodotto di h e una stima della pendenza di f . Nellospecifico:

• k1 è l’incremento basato sulla pendenza all’inizio dell’intervallo, utilizzando y(metodo di Eulero)

• k2 è l’incremento basato sulla pendenza alla metà dell’intervallo, utilizzandoy + h

2k1

• k3 è un altro incremento basato sulla pendenza alla metà dell’intervallo, utiliz-zando y + h

2k2

6

• k4 è l’incremento basato sulla pendenza alla fine dell’intervallo, utilizzando y+hk3

Nel fare la media, gli incrementi valutati in un punto intermedio dell’intervallo hannopeso maggiore, ed i coefficienti sono scelti in modo che se f è indipendente da y, sicchél’equazione dipende da un semplice integrale, allora il metodo RK coincide con la regoladi Cavalieri-Simpson.

3.2 Metodo beta-NewmarkL’algoritmo più vecchio e probabilmente più ampiamente utilizzato per l’integrazione delleequazioni delle dinamiche strutturali è quello di Newmark (1959). In questo caso talealgoritmo viene discusso in relazione all’equazione lineare del moto (3.2).

L’obiettivo della procedura di integrazione Newmark è il calcolo del vettore di stato(un+1, un+1) al tempo tn+1 = tn + h, dato il vettore di stato (un, un) allo step precedentetn e il vettore di carico in entrambe le volte, fn e fn+1. La procedura viene stabilitain due passaggi: prima gli incrementi di u(t) e u(t) sono espressi in termini di integralidell’accelerazione u nell’intervallo di tempo [tn, tn+1], e quindi l’accelerazione viene valu-tata dall’uso delle equazioni del moto. Pertanto, l’algoritmo di Newmark si basa su dueprincipi: l’espressione di spostamento e incrementi di velocità in termini di accelerazionecon una rappresentazione approssimata – asintoticamente valida per piccoli incrementi ditempo – e collocazione dell’equazione del moto al successivo step temporale tn+1. Appro-priati parametri adottati nelle rappresentazioni di spostamenti e velocità sono determinaticonsiderando il problema lineare, e il metodo viene quindi esteso a problemi non linearirisolvendo l’equazione non lineare del movimento per iterazione.

Segue immediatamente dal calcolo integrale che il valore del vettore u(tn+1) al tempo(tn+1) può essere espresso come il suo valore u(tn) al tempo tn più un incremento valutatocome integrale della derivata temporale u(τ) nell’intervallo di tempo:

u(tn + h) = u(tn) +

∫ tn+1

tn

u(τ) d(τ) (3.6)

Ora il fattore 1 può essere introdotto sotto il segno dell’integrale e lo stesso può essererisolto tramite integrazione per parti. Quando l’integrale del fattore 1 è selezionato comeτ − tn+1, si ha:

(??)u(tn+1) = u(tn)− [(tn+1 − τ)u(τ)]tn+1

tn +

∫ tn+1

tn

(tn+1 − τ)u(τ) d(τ) (3.7)

Si noti che il limite superiore tra parentesi quadre non contribuisce a causa del fattoretn+1 − τ .

Nella procedura Newmark un punto centrale è la rappresentazione degli incrementidi velocità e spostamento in termini di integrali dell’accelerazione. La prima formula èottenuta da (3.6), quando si sostituiscono u e u con u e u, rispettivamente. La secondaformula deriva da (??) dalla valutazione del termine tra parentesi quadre:

un+1 = un +

∫ tn+1

tn

u(τ) d(τ)

un+1 = un + hun +

∫ tn+1

tn

(tn+1 − τ)u(τ) d(τ)

(3.8)

dove h = tn+1 − tn è la lunghezza dell’intervallo di tempo e un = u(tn).

7

Nel presente uso della formula di integrazione l’accelerazione non sarà nota nell’interointervallo e gli integrali dovranno essere valutati approssimativamente dal valore del vet-tore di accelerazione agli estremi dell’intervallo. Le formule appropriate sono date dallemedie ponderate dell’accelerazione nella forma∫ tn+1

tn

u(τ) d(τ) ' (1− γ)hun + γhun+1 (3.9a)∫ tn+1

tn

(tn+1 − τ)u(τ) d(τ) ' (1

2− β)h2un + βh2un+1 (3.9b)

I parametri 0 < γ < 1 e 0 < β < 12determinano il grado di avanzamento della ponderazio-

ne, con γ = 0,β = 0 che corrisponde alla ponderazione totale all’indietro, e γ = 1,β = 12

che corrisponde alla piena ponderazione in avanti. Si noti che in (3.9a) la somma dei pesidi integrazione è h, mentre in (3.9b) la somma è 1

2h a causa del fattore lineare nell’inte-

grando. Queste proprietà sono necessarie per garantire convergenza al diminuire di h, ilche implica che le formule debbano essere esatte per un’accelerazione costante.

Quando le formule di integrazione approssimata vengono sostituite nelle si ottengonole seguenti relazioni discrete tra i vettori spostamento, velocità e accelerazione ai due puntifinali dell’intervallo:

un+1 = un + (1− γ)hun + γhun+1, (3.10a)

un+1 = un + hun +

(1

2− β

)h2un + βh2un+1. (3.10b)

L’algoritmo di integrazione Newmark si ottiene soddisfacendo le equazioni del moto altempo tn+1, usando la rappresentazione e per il vettore di velocità e spostamento, ri-spettivamente. La sostituzione di queste rappresentazioni nelle equazioni lineari del motofornisce la seguente equazione per il vettore accelerazione un+1:(M + γhC + βh2K

)un+1 = fn+1 − C

(un + (1− γ)hun

)−K

(un + hun + (

1

2− β)h2un

)(3.11)

Questa equazione consente di calcolare il vettore accelerazione un+1 allo step tn+ 1, ela velocità e lo spostamento possono essere calcolati dalla .

È conveniente organizzare il calcolo di un passo temporale nell’algoritmo di Newmarksotto forma di una fase di predizione seguita da una fase di correzione. Nel problemalineare la definizione dei predittivi è semplicemente una questione di convenienza, mentrenei problemi non lineari il predittivo funge da punto di partenza per le iterazioni. Nellafase di predizione valori preliminari u∗n+1 e u∗n+1di velocità e spostamento sono valutatidalla , senza l’ultimo termine contenente l’accelerazione ancora sconosciuta un+1:

u∗n+1 = un + (1− γ)hun, (3.12a)

u∗n+1 = un + hun +

(1

2− β

)h2un. (3.12b)

Questi valori predittivi figurano direttamente alla destra dell’equazione per l’accelerazioneun+1. Con la notazione

M∗ = M + γhC + βh2K (3.13)

per la matrice delle masse modificata, l’equazione del moto prende la forma semplificata

M∗un+1 = fn+1 − Cu∗n+1 −Ku∗n+1 (3.14)

8

dove quest’ultima equazione è risolta per l’accelerazione un+1.La correzione consiste nell’aggiungere l’ultimo termine nella al predittivo vettore di velo-cità e spostamento:

un+1 = u∗n+1 + γhun+1, (3.15a)un+1 = u∗n+1 + βh2un+1. (3.15b)

completando così la fase temporale.Per le strutture non smorzate è possibile formulare l’algoritmo in forma esplicita se

la matrice delle masse è diagonale e β = 0. In tal caso la può essere risolta diretta-mente per un+1, senza l’inversione della matrice. Questo particolare caso è stabile solocondizionalmente e richiede l’incremento temporale h entro un limite superiore.

9

Capitolo 4

Confronto

Nell’analisi computazionale tramite MATLAB sono stati analizzati sia il caso lineare siaquello non lineare.I parametri presi in considerazione sono i seguenti:

• m = 1.0 massa

• k = 15.0 rigidezza lineare

• ζ variabile tra 0.001 e 0.05

• f(t) = Acos(ωt) forzante esterna, con ampiezza A = 0.05 e frequenza ω = 0.9

Per il caso non lineare è stato preso un rapporto tra le rigidezze variabile

a =knlk

(4.1)

dove per knl si intende proprio l’apporto di rigidezza non lineare.

Di seguito gli output del comando ODE45 e gli output dell’algoritmo di Newmark

10

0 50 100 150 200 250 300 350 400

time

-0.025

-0.02

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

0.02

0.025

deflection

Figura 4.1: Spostamenti nel tempo

11

0 50 100 150 200 250 300 350 400

time

-0.1

-0.08

-0.06

-0.04

-0.02

0

0.02

0.04

0.06

0.08

0.1

velocity

Figura 4.2: Velocità nel tempo

12

Figura 4.3: Fasi

13

Figura 4.4: Forza-spostamenti: non linearità

14

0 50 100 150 200 250 300 350 400

time

-0.025

-0.02

-0.015

-0.01

-0.005

0

0.005

0.01

0.015

0.02

0.025

deflection

Figura 4.5: Spostamenti nel tempo

15

Figura 4.6: Spostamenti nel tempo tramite iterazione di Newton-Raphson

16

Conclusioni

Lo scopo per il quale si è deciso di usare il software MATLAB è quello di potersi ricavarein forma numerica la convergenza del problema fin qui delineato tramite iterazione. Inquesto modo si è ragionato su come sarebbe potuta procedere l’analisi per il prosieguodella tesi su problemi di carattere isteretico.Tale processo porterà all’implemetenzaione di un nuovo algoritmo che, servendosi dell’ap-profondito studio fatto del metodo beta-Newmark, richiederà lo studio per l’aggiunta diun nuovo parametro che caratterizzerà, tramite discretizzazione temporale in un intervalloben definito, la forza di richiamo in campo plastico.

17