POLITECNICO DI MILANO

Corso di Studi in Ingegneria Matematica

Progetto di programmazione avanzata

Problemi di contatto in elasticita lineare

Umberto Villa Matr. 680022

ANNO ACCADEMICO 2006-2007

Indice

1 Il problema di Signorini 31.1 Statica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.1.1 Notazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31.1.2 Problema di Signorini . . . . . . . . . . . . . . . . . . . . . . . . . 41.1.3 Algoritmo risolutivo . . . . . . . . . . . . . . . . . . . . . . . . . . 6

1.2 Dinamica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.2.1 Notazioni . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81.2.2 Discretizzazione temporale . . . . . . . . . . . . . . . . . . . . . . 81.2.3 Problema variazionale . . . . . . . . . . . . . . . . . . . . . . . . . 101.2.4 Approssimazione numerica del problema variazionale . . . . . . . . 101.2.5 Sistema accoppiato di due corpi . . . . . . . . . . . . . . . . . . . . 11

2 Implementazione 132.1 Struttura delle classi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 132.2 Classi per la gestione di corpi deformabili e rigidi . . . . . . . . . . . . . . 14

2.2.1 La classe ElasticStructureSolver . . . . . . . . . . . . . . . . . . . . 142.2.2 La classe RigidBody . . . . . . . . . . . . . . . . . . . . . . . . . . 172.2.3 La classe DataNewmark . . . . . . . . . . . . . . . . . . . . . . . . 19

2.3 Classi per la gestione del contatto . . . . . . . . . . . . . . . . . . . . . . . 202.3.1 La classe Contact . . . . . . . . . . . . . . . . . . . . . . . . . . . . 202.3.2 La classe ContactFriction . . . . . . . . . . . . . . . . . . . . . . . 242.3.3 La classe SlipFriction . . . . . . . . . . . . . . . . . . . . . . . . . . 252.3.4 La classe GeneralContact . . . . . . . . . . . . . . . . . . . . . . . 262.3.5 La classe VincoloHD . . . . . . . . . . . . . . . . . . . . . . . . . . 272.3.6 La classe Vincolo . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

2.4 Classi di supporto . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 292.4.1 Class AnaliticalFunction . . . . . . . . . . . . . . . . . . . . . . . . 29

3 Simulazioni 313.1 Sfera su basamento fisso . . . . . . . . . . . . . . . . . . . . . . . . . . . . 313.2 Sfera su basamento mobile . . . . . . . . . . . . . . . . . . . . . . . . . . . 333.3 Sfera in una scatola . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 343.4 Attrito disco-pastiglia . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

1

Introduzione

Nel presente progetto viene proposto ed implementato un approccio variazionale per lasoluzione del problema del contatto, anche in presenza di attrito, tra un corpo elasticoed un corpo rigido, la cui geometria e descritta da un piano.

In particolare si utilizzano i moltiplicatori di Lagrange per la minimizzazione delfunzionale energia di deformazione del corpo elastico sotto il vincolo di non compenetra-zione. Tali moltiplicatori assumono un chiaro significato fisico, rappresentando gli sforzinormali scambiati tra le due superfici a contatto. La tecnica degli active set riducel’onere computazionale, risolvendo tale problema di ottimizzazione senza aumentare ladimensione del sistema lineare.

Per quanto riguarda l’avanzamento temporale si utilizza il metodo di Newmark aver-age acceleration, accurato al secondo ordine, ma poco adatto all’integrazione di proble-mi stiff, caratterizzati dalla presenza di oscillazioni numeriche ad alte frequenze indottedalla discretizzazione spaziale piu che dalla fisica del fenomeno. Prossimamente verraimplementato un metodo generalized α, in grado di dissipare correttamente tali modispuri.

La presente relazione e suddivisa in tre parti. La prima fornisce l’ambientazioneteorica del problema e mostra l’algoritmo risolutivo: per maggiore chiarezza si analizzadapprima la statica e successivamente la dinamica. La seconda spiega come sono im-plementate le diverse classi, con particolare attenzione alle loro dipendenze reciproche(ereditarieta/inclusione), alle interfacce utente e strutture dati. La terza presenta irisultati delle simulazioni: il primo caso test controlla la solidita dell’algoritmo e la suaaderenza ai principi fisici di base (conservazione-dissipazione energia cinetica) per unproblema di una sfera che rimbalza su un fondamento rigido, il secondo mostra l’accop-piamento tramite contatto nella dinamica di due corpi, il terzo risolve un problema dallegeometrie piu complesse. Infine si effettua una simulazione di una prova al banco perun impianto frenante, in cui si calcola il contatto e l’attrito tra pastiglia e disco. [5] [2][6] [4] [3][1][7][8]

2

Capitolo 1

Il problema di Signorini

1.1 Statica

1.1.1 Notazioni

Il tensore delle piccole deformazioni εij e definito in funzione del vettore spostamentou = (u1, u2, u3) dalla

εij(u) =12

(∂ui

∂xj+

∂uj

∂xi

)(1.1)

Per le ipotesi di elasticita lineare, il legame tra il tensore degli sforzi e quello dideformazioni e il seguente:

σij = Eijklεkl (1.2)

con Eijkl = Ejikl, Eijkl = Eijlk, Eijkl = Elkij . Nel caso di legame elastico lineareisotropo, la relazione tra σ e ε e

σij = λδij div u + 2µεij (1.3)

dove λ e µ sono dette costanti di Lame.

Sia Ω la regione di spazio occupata dal corpo elastico, ΓD ⊂ ∂Ω la superficie vincolatadel corpo, ΓN ⊂ ∂Ω la superficie caricata, consideriamo il seguente problema:

σij,j = fi in Ωσijnj = ti su ΓN

ui = 0 su ΓD

dove fi rappresentano le forze di volume, ti le forze di superficie (positive se ditrazione, negative se di compressione).

Definiamo lo spazio V :=v ∈ [H1(Ω)]3 : v = 0 su ΓD

, e l’operatore γ(·) che

associata ad ogni funzione di V la sua traccia in [H1/2(∂Ω)]3.

3

Mediante il principio dei lavori virtuali, il problema elastico non vincolato puo esserecosı formalizzato:

Problema 1.1 Trovare u ∈ V tale che∫

ΩEijklεij(u)εkl(v) =

∫

Ωfivi +

∫

ΓN

tiγ(vi) ∀v ∈ V

con f ∈ [L2(Ω)]3 e t ∈ [H−1/2(∂Ω)]3.

Introduciamo, infine, l’energia elastica di deformazione w

w =12σijεij =

12Eijklεijεkl (1.4)

e il funzionale energia J che tiene conto anche delle forze esterne

J(v) =∫

Ωw(v)−

∫

Ωf · v −

∫

ΓN

t · γ(v) (1.5)

1.1.2 Problema di Signorini

Detta ΓC la superficie di contatto sottoposta al vincolo cinematico di non compenetra-zione, il problema di Signorini si scrive nella seguente forma:

σij,j = fi in Ω (1.6)ui = 0 su ΓD (1.7)

σijnj = ti su ΓN (1.8)un − g ≤ 0 su ΓC (1.9)

σn ≤ 0 su ΓC (1.10)(un − g)σn = 0 su ΓC (1.11)

dove un = uini rappresenta lo spostamento in direzione normale alla superficie, σn =σijnjni lo sforzo normale. In particolare le condizioni imposte su ΓC sono le seguenti:

• vincolo cinematico di non compenetrazione (1.9)(per la sua derivazione vedi [4])

• contatto unilatero: sono ammissibili solo forze di compressione (1.10)

• scarti complementari: se non c’e contatto non c’e nemmeno sforzo normale (1.11).

Definta C := v ∈ V : γn(v) ≤ g1 la regione ammissibile per la soluzione u, scri-viamo il problema elastico in presenza di vincoli cinematici come un problema diottimizzazione:

1γn e l’operatore di traccia normale γn(v) = γ(vini).

4

Problema 1.2 Trovare u ∈ V tale che

u = argminv∈C

J(v) = argminv∈C

∫

Ωw(v)−

∫

Ωf · v −

∫

ΓN

t · γ(v)

Per risolvere il problema (1.2), utilizziamo un metodo iterativo che consiste nellarisoluzione successiva di problemi di ottimizzazione con vincoli di uguaglianza. A talescopo definiamo lo spazio X dei moltiplicatori di Lagrange λ associati al problema (1.2):

X = (H1/2(∂Ω))′ (1.12)

Il metodo che viene proposto e la controparte continua dei metodi di tipo active setutilizzati in ricerca operativa nella soluzione di problemi di programmazione quadratica.In virtu di tale interpretazione tratteremo le proprieta di convergenza direttamente peril problema discreto.

• porre Γ(0)AC = ΓC , essendo ΓAC ⊂ ΓC la superficie effettiva di contatto tra i due

corpi

• iterare i=0....

– trovare (u(i), λ(i)) punto sella della lagrangiana:

L(i)(v, µ) = J(v)−∫

ΓAC(i)

(µ(g − γn(v))) (1.13)

– aggiornare la superficie effettiva di contatto 2:

Γ(i+1)AC =

x ∈ ΓC : λ(i) + c(γn(u(i))− g) ≤ 0

(1.14)

L’interpretazione del moltiplicatore di Lagrange λ, come lo sforzo normale σn scam-biato tra il corpo elastico e la superficie di contatto, si ha scrivendo esplicitamente ilsistema lineare per la ricerca del punto-sella.

∂L

∂u=

∫

ΩEijkluk,lvi,j −

∫

Ωfivi −

∫

ΓN

tiγ(vi)−∫

ΓAC

λγn(v) = 0 ∀v ∈ V

∂L

∂λ=

∫

ΓAC

µ(γn(u)− g) = 0 ∀µ ∈ X

2essendo c una costante di penalizzazione per il vincolo di impenetrabilita

5

1.1.3 Algoritmo risolutivo

Per risolvere il problema (1.2) usiamo un approccio del tipo discretize then optimize.I vantaggi di tale approccio risiedono nella particolare struttura del problema discretocosı ottenuto, per il quale e possibile applicare algoritmi della ricerca operativa assaiefficienti.Applicando una discretizzazione tipo elementi finiti del problema (1.2) si ha il seguente

Problema 1.3 Trovare u ∈ RN tale che

u = arg minu

Jh(u) =12u′Au− f ′u

sotto il vincolo (un)i ≤ gi ∀nodo i di bordo.

u e il vettore dei coefficienti degli elementi di base nella soluzione, A la matrice distiffness dell’elasticita lineare, f la discretizzazione del funzionale lineare legato alle forzedi volume e di superficie, g il vettore dei coefficienti dell’interpolante elementi finiti delvincolo g, N il numero di gradi di liberta dello spazio Vh.

E dunque riconoscibile la struttura di problema di programmazione quadratica, nelquale si deve ottimizzare una funzione quadratica strettamente convessa sotto vincolilineari.

Per tale categoria di problemi risultano molto efficienti i metodi di tipo active set,nei quali si determina in modo iterativo l’insieme dei vincoli attivi (ossia quelli in cuivale l’uguaglianza e non la disuguaglianza stretta) tramite il controllo sull’ammissibilitadella soluzione e del segno del moltiplicatore associato al vincolo.

Nel dettaglio i passi dell’algoritmo sono i seguenti [2]:

Algoritmo 1.1 Metodo active set per vincoli di minore-uguale

1. Inizializzazione: porre k=0, scegliere λ0 ∈ RN , u0 ∈ RN .

2. Ciclo su k=0,....

(a) calcolo dell’active set3 Ik = i ∈ N : λk,i + c(gi − uk,i) ≤ 0(b) calcolo del suo complementare Jk = N \ Ik

(c) costruzione della matrice4 A e del vettore f :

Aij = Aij se i ∈ Jk, Aij = δij se i ∈ Ik

fi = fi se i ∈ Jk, fi = gi se i ∈ Ik

3c e una costante di penalizzazione che forza il rispetto del vincolo4Tale costruzione puo essere fatta in modo efficiente e veloce utilizzando un metodo di penalizzazione

analogo a quello per imporre le condizioni di Dirichlet sui nodi di bordo

6

(d) soluzione del sistema lineare Auk+1 = f(e) aggiornamento dei moltiplicatori: λk+1 = Auk+1 − f

In riferimento alla definizione dell’active set Ik notiamo che

• se il vincolo i-esimo non e soddisfatto alla k-esima iterazione, gi−uk,i ≤ 0 e dunque,considerato il valore elevato della costante di penalizzazione c, il vincolo i-esimoviene incluso in Ik;

• se il vincolo i-esimo e soddisfatto strettamente alla k-esima iterazione esso nonviene incluso in Ik a prescindere dal valore del moltiplicatore ad esso associato;

• se il vincolo i-esimo e soddisfatto con l’uguaglianza alla k-esima iterazione alloraesso viene incluso in Ik se e solo se il suo moltiplicatore ha segno negativo (seλi ≥ 0 la soluzione e ammissibile ma non ottimale).

E possibile dimostrare che tale scelta di aggiornamento dell’active set permette ilraggiungimento delle condizioni di ottimo. Indicata con u∗,λ∗,I∗,J∗ la soluzione a cuiconverge l’algoritmo, si ha infatti

• ammissibilita della soluzione: u∗i ≤ gi ∀i ∈ N ;

• segno dei moltiplicatori: λ∗i ≤ 0 ∀i ∈ N ;

• scarti complementari: λ∗i (u∗i − gi) = 0

Immediata e l’estesione dell’algoritmo a vincoli di maggiore-uguale del tipo u ≥ l:

Algoritmo 1.2 Metodo active set generale

1. Inizializzazione: porre k=0, scegliere λ0 ∈ RN , u0 ∈ RN .

2. Ciclo su k=0,....

(a) calcolo degli active sets

I−k = i ∈ N : λk,i + c(gi − uk,i) ≤ 0, I+k = i ∈ N : λk,i + c(li − uk,i) ≥ 0

(b) calcolo dell’insieme complementare Jk = N \ (I+k

⋃I−k )

(c) costruzione della matrice A e del vettore f :

Aij = Aij se i ∈ Jk, Aij = δij se i ∈ I+k

⋃I−k

fi = fi se i ∈ Jk, fi = gi se i ∈ I−k , fi = li se i ∈ I+k

(d) soluzione del sistema lineare Auk+1 = f .(e) aggiornamento dei moltiplicatori: λk+1 = Auk+1 − f

In questo secondo caso avremo

• ammissibilita della soluzione: li ≤ u∗i ≤ gi ∀i ∈ N ;

• segno dei moltiplicatori: λ∗i ≤ 0 ∀i ∈ I−∗, λ∗i ≥ 0 ∀i ∈ I+∗;

• scarti complementari: λ∗i (u∗i − gi)(u∗i − li) = 0

7

1.2 Dinamica

1.2.1 Notazioni

Il moto di un corpo lineare elastico in presenza di attrito viscoso puo essere descrittocon la seguente equazione alle derivate parziali:

utt − Cut −Ku = f in QT = Ω× [0, T ]u(x, 0) = g(x), ut(x, 0) = h(x) in Ω

u(x, t) = 0 su ST = ∂Ω× [0, T ](1.15)

dove l’operatore C e lo smorzamento di Rayleigh:

C = aI + bK

essendo I l’operatore identita, a e b parametri che determinano lo smorzamento ξ allediverse frequenze di vibrazione ω del sistema. In particolare il parametro a ha un pesomaggiore nello smorzamento delle basse frequenze, b delle alte frequenze:

ξ(ω) =a/ω + bω

2(1.16)

1.2.2 Discretizzazione temporale

Una semidiscretizzazione in spazio del problema (1.15) conduce alla risoluzione di unsistema di equazioni differenziali ordinarie della forma

Mutt + Cut + Ku(t) = M f(t) (1.17)

dove M e la matrice di massa, K la matrice di stiffness associata al problema elasticolineare stazionario, C = aM + bK la matrice di smorzamento viscoso.

Per integrare in tempo il sistema (1.17) useremo lo schema di Newmark che e unodei piu usati in meccanica strutturale. Detto ∆t il passo di avanzamento temporale,siano un, vn, an le approssimazioni della soluzione u(t) e delle sue derivate temporaliall’istante n∆t.

La famiglia dei metodi di Newmark e definita in funzione di due parametri β e γ chedeterminano le proprieta di stabilita ed accuratezza del metodo. Lo schema nella formagenerale consiste nelle seguenti equazioni:

Man+1 + Kun+1 = M fn+1 (1.18)

un+1 = un + ∆tvn +∆t2

2[(1− 2β)an + 2βan+1] (1.19)

vn+1 = vn + ∆t[(1− γ)an + γan+1] (1.20)

Un metodo efficiente per implementare tale metodo e dato dal seguente algoritmo5:5Per verificarne la consistenza e sufficiente sostituire le definizioni di spostamento (1.19) e velocita

(1.20) nell’equazione per l’accelerazione (1.18)

8

Algoritmo 1.3 Metodo di Newmark:

1. Inizializzazione t = 0, Ma0 = F0 − Cv0 −Ku0.

2. Ciclo su n=0,1....

(a) Fase predictor

un+1 = un + ∆tvn + (1− 2β)∆t2

2an (1.21)

vn+1 = vn + (1− γ)∆tan (1.22)

(b) Fase eval (calcolo dell’accelerazione):(M + γ∆tC + β∆t2K

)an+1 = Fn+1 − Cvn+1 −Kun+1 (1.23)

(c) Fase correctorun+1 = un+1 + β∆t2an+1 (1.24)

vn+1 = vn+1 + γ∆tan+1 (1.25)

Riportiamo,infine, un risultato di stabilita per la famiglia dei metodi di Newmark[3]:

• incondizionatamente asintoticamente stabile

2β ≥ γ ≥ 12

• condizionatamente asintoticamente stabile

γ ≥ 12

β <γ

2

con regione di stabilita:

∆t ≤ Ωcrit

ω

dove ω e la massima frequenza di oscillazione del sistema e Ωcrit e la frequenza dicampionamento critica:

Ωcrit =(γ

2− β

)− 12

Di particolare interesse sono il metodo average acceleration (β = 14 , γ = 1

2) ed ilmetodo differenze centrate(β = 0, γ = 1

2):

• average acceleration e un metodo implicito che approssima il moto in ogni inter-vallo di tempo ∆t con la legge oraria di un moto uniformemente accelerato conaccelerazione pari alla media tra l’istante n∆t e (n + 1)∆t. Garantisce la stabilitaincondizionata ed ordine di accuratezza O(∆t2) in velocita e spostamento.

9

• differenze centrate e un metodo esplicito, ossia non utilizza l’accelerazione al-l’istante attuale nel calcolo dello spostamento. E condizionatamente stabile(Ωcrit = 2) e con ordine di accuratezza O(∆t2).

Nel problema del contatto tra un corpo rigido ed un corpo viscoelastico che analiz-zeremo in seguito, useremo lo schema differenze centrate per il corpo rigido ed averageacceleration per il corpo deformabile. Con tale approccio, come vedremo, e infatti pos-sibile risolvere le equazioni del moto senza sottoiterazioni nel calcolo della posizione deicorpi a ciascun istante temporale.

1.2.3 Problema variazionale

Consideriamo il seguente problema di ottimizzazione vincolata nel quale per semplicitatrascuriamo il termine di smorzamento:

Problema 1.4 trovare u appartenente a

C =v ∈ L2(0, T ; [H1

0 (Ω)]3), γn(v(t)) ≤ g(t)∀t ∈ [0, T ]

, ut ∈ L2(0, T ; L2(Ω))

tale che

u(t) = arg minv(t)∈C

J(v) =12(vt,vt) +

12a(v(t),v(t))− < f(t),u(t) > ∀t ∈ [0, T ]

dove a(·, ·) e la forma bilineare simmetrica, continua, coerciva del problema dell’elasticitalineare; (·, ·) rappresenta il prodotto scalare in [L2(Ω)]3 e < ·, · > la dualita tra [H1

0 (Ω)]3

e [H−1(Ω)]3; f ∈ L2(0, T ; [H−1(Ω)]3) e g ∈ L2(0, T ; H10 (Ω)) sono funzioni assegnate.

1.2.4 Approssimazione numerica del problema variazionale

Data la grande difficolta teorica nell’analisi del problema di ottimizzazione vincolata nelcaso non stazionario, utilizzeremo l’approccio discretize-then-optimize.

Definiamo il funzionale discreto energia Jh(an) in modo tale che l’equazione (1.23)traduca ∇Jh = 0.

Jh(an) =12a′n(M + γ∆tC + β∆t2K)an − F′nan + (Cvn)′an + (Kun)′an (1.26)

dove vn, un sono definite in accordo a (1.22), (1.21).Infine il vincolo γn(u)(i)n ≤ g

(i)n ((i) indica l’i-esima componente del vettore) puo essere

riscritto in funzione di a(i)n grazie all’equazione (1.24) nel seguente modo:

γn(a)(i)n ≤ g(i)n g(i)

n =g(i)n − γn(u)(i)n

β∆t2(1.27)

Arriviamo dunque al seguente problema di programmazione quadratica:

10

Problema 1.5 Trovare an ∈ (R3)n tale che

an = argminJh(an)

sotto il vincolo γn(a)(i)n ≤ g(i)n ∀i.

Tale problema viene infine risolto mediante l’algoritmo active set (1.1).

1.2.5 Sistema accoppiato di due corpi

In questo paragrafo analizziamo come risolvere il problema dell’impatto tra un corpoelastico e uno rigido. Ad ogni istante temporale, note le forze agenti all’istante prece-dente, calcoliamo la posizione attuale del corpo rigido mediante il metodo delle differenzecentrate. Tale informazione viene utilizzata per calcolare il vincolo g che permette diottenere tramite la soluzione di (1.5), le accelerazioni del corpo deformabile e la forza dicontatto scambiata. Infine si integra il moto del corpo deformabile mediante il metodoaverage acceleration.

Nel seguito indicheremo con

• ud,vd,ad il vettore spostamento, velocita e accelerazione del corpo elastico

• ur,vr,ar spostamento, velocita e accelerazione del baricentro del corpo rigido

• g il profilo del corpo rigido

• Ik l’insieme dei nodi del corpo elastico a contatto con il corpo rigido

• σn le pressioni di contatto esercitate dal corpo rigido sul corpo elastico, e λ imoltiplicatori di Lagrange ad esse associati mediante la formula:

λ(i) =∫

∂Ω(σn)γn(v(i)) ∀i

• n la normale alla superficie di contatto

• Fc la forza di contatto esercitata dal corpo elastico sul corpo rigido

• m la massa del corpo rigido.

Presentiamo ora come viene implementato l’accoppiamento dei due corpi:

Algoritmo 1.4 Accoppiamento metodo average acceleration e differenze centrate:

1. Inizializzazione di ud0, vd

0, ad0, ur

0, vr0, ar

0;

2. Ciclo temporale n=0,1...

11

(a) Fase predictor

udn+1 = ud

n + ∆tvdn +

∆t2

4ad

n

vdn+1 = vd

n +∆t

2ad

n

urn+1 = ur

n + ∆tvrn +

∆t2

2ar

n

vrn+1 = vr

n +∆t2

2ar

n

(b) Fase eval

• Valutazione del vincolo

gn+1 =(g)− γn(ud

n+1 + x)∆t2/4

• Calcolo dell’accelerazione adn+1 del corpo deformabile, della regione di

contatto Ik e dei moltiplicatori di Lagrange discreti λ mediante lasoluzione di (1.5)

• Valutazione della forza esercitata dal corpo deformabile al corpo rigido 6

Fc =∫

∂Ω(σn)n ≈

∑

i∈Ik

λ(i)n+1n

(i) arn+1 =

Fc

mr

(c) Fase corrector

udn+1 = ud

n+1 +∆t2

4ad

n+1

vdn+1 = vd

n+1 +∆t

2ad

n+1

urn+1 = ur

n+1

vrn+1 = vr

n+1 +∆t2

2ar

n+1

6Nel caso in cui la discretizzazione EF rappresenti una partizione dell’unita vale l’equivalenza

12

Capitolo 2

Implementazione

Il codice e stato implementato sulla base delle librerie gia presenti in LifeV [1]. Inparticolare sono state utilizzate le classi seguenti:

• DataMesh, RegionMesh3D, RegionMesh3D ALE per la descrizione della mesh dicalcolo

• RefFe, CurrentFe, Dof, QuadRule, ElemMat, ElemVec, BCHandler per la de-scrizione degli elementi finiti e della routine di integrazione numerica perl’assemblaggio delle matrici ad elementi finiti

• MSRPatt, MSRMatr, ScalUnknown, PhysVectUnknown, KN per la descrizione dimatrici sparse e vettori

• GetPot per la lettura da file

• SolverAtzec per l’interfacciamento con il solutore per sistemi lineari Atzec [7].

Il contributo originale del lavoro consiste nei seguenti aspetti

• implementazione ex novo di un algoritmo di contatto tra un corpo rigido ed unoelastico mediante un approccio di tipo variazionale

• modellazione dell’attrito dinamico tra due corpi a contatto in funzione dellavelocita di slip tra di essi

• creazione di una classe per la descrizione del moto (anche vincolato) di un corporigido

• implementazione della famiglia di metodi di Newmark per l’avanzamento in tempo.

2.1 Struttura delle classi

Le classi implementate possono essere divise in tre diversi gruppi a seconda della lorofunzione:

13

• classi per la gestione di corpi deformabili e rigidi: DataNewmark, ElasticStructure-Solver, RigidBody

• classi per la modellazione del contatto tra corpi: Contact, ContactFriction,FrictionSlip,GeneralContact,VincoloHD,Vincolo

• classi di supporto, legate all’interpolazione di funzioni analitiche AnaliticalFunctione classi figlie.

Figura 2.1: Struttura delle classi implementate: linea continua ereditarieta, lineatratteggiata inclusione

2.2 Classi per la gestione di corpi deformabili e rigidi

Tali classi permettono l’integrazione delle equazioni del moto sia per un corpo rigido cheper un corpo elastico. Esse presentano un’interfaccia utente il piu possibile omogeneaper corpi di natura diversa. Le informazioni e gli algoritmi riguardanti l’avanzamento intempo sono contenute in una classe template virtuale DataNewmark capace di gestiresia corpi rigidi che deformabili.

2.2.1 La classe ElasticStructureSolver

Tale classe e le sue classi genitori (DataElasticStructure, ElasticStructureHandler,DataNewmark) contengono gli attributi e i metodi per la soluzione di problemi dielasticita lineare tempo dipendenti mediante il metodo di Newmark.

14

Figura 2.2: Diagramma di ereditarieta per il solutore delle strutture; le caselletratteggiate contengono classi virtuali

Parametri template

La classe contiene un parametro template di tipo Mesh.

Attributi

Elenchiamo i principali attributi della classe e delle sue classi genitori:

• classe DataElasticStruture:

– caratteristiche fisiche del corpo elastico: densita ( rho), modulo di Young( E ), coefficiente di Poisson ( nu), costanti di Lame ( lambda, mu)

Figura 2.3: Struttura delle classi implementate: linea blu ereditarieta, linea nerainclusione.

15

• classe ElasticStrutureHandler :

– oggetti per la discretizzazione ad elementi finiti: elemento finito di riferimen-to ( refFe), descrizione dei gradi di liberta del problema ( dof ), regole diquadratura ( Qr, bdQr), elemento finito corrente ( fe, feBd)

– vettori spostamento ( d), velocita ( w), accelerazione ( a)

• classe ElasticStrutureSolver :

– matrici sparse per la descrizione del problema elastico: pattern ( pattK ),matrice di stiffness ( K ), massa ( M ), damping ( D), matrice del sistemalineare di Newmark con condizioni al bordo ( A) e senza ( AConst)

– termine noto del sistema lineare di Newmark ( rhs)

– matrici e vettori locali per l’assemblaggio delle matrici globali: elmatK,elmatM, elmatD, elmatA, elvec

– linearSolver : solutore del sistema lineare di Newmark tramite la libreriaAtzec.

Metodi

• Costruttore: inizializza gli attributi tramite GetPot, assembla le matrici globalidel problema. Non e fornito ne il costruttore di default ne quello di copia.

• Setters&Getters per i principali attributi.

• void initialize(const Function & d0, const Function & w0): permette diinizializzare il vettore spostamento d e velocita w mediante una funzione f conla seguente firma: Real f(const Real & t, const Real & x, const Real & y,const Real & z, const ID & i).

• void initialize(const std::string & depName, const std::string & vel-Name): permette di inizializzare il vettore spostamento d e velocita w mediantelettura da bb-file.

• void timeAdvance(source type const & source, const Real & time):gestisce l’avanzamento di uno step temporale con il metodo di Newmark.source type e un funtore con la seguente firma: Real f(const Real & t, constReal & x, const Real & y, const Real & z, const ID & i). Ad ogni istantetemporale vengono svolte le seguenti operazioni:

– fase predictor : si invoca DataNewmark::pred(...)

– fase eval :

∗ calcolo del termine noto del sistema: valutazione delle forze divolume (computeSource()); aggiunta del contributo di Newmark(DataNewmark::computeRhs(...))

16

∗ aggiunta delle condizioni al bordo (eSSaddBC())∗ soluzione del sistema lineare (eval())

– fase corrector : si invoca DataNewmark::corr(...).

• void computeSource(): assembla il termine noto relativo alle forze di volumetramite un ciclo su ciascun elemento finito della mesh.

• void aSSaddBC(): gestisce le condizioni al bordo tramite la funzione bcMan-age(...).

• void eval(): risolve il sistema lineare associato al metodo di Newmark tramite lalibreria Atzec.

• funzioni per il calcolo di grandezze globali relative al baricentro:

– posizione → Real centerOfMass(UInt comp):

xcm =1iMi,j(xj + dj)

1iMi,j1j

– velocita → Real centerOfMassW(UInt comp):

wcm =1iMi,jwj

1iMi,j1j

– accelerazione → Real centerOfMassA(UInt comp):

acm =1iMi,jaj

1iMi,j1j

dove componente vale: 0 X, 1 Y, 2 Z.

2.2.2 La classe RigidBody

Tale classe contiene gli attributi ed i metodi per soluzione delle equazioni del moto perun corpo rigido vincolato tramite lo schema di Newmark. Tale classe eredita dalla classeDataNewmark.

Paramentri template

Tale classe ha un parametro template di tipo vettore. Per tale parametro templatedevono essere implementati i seguenti operatori:

• operatore di assegnazione: T & operator=(T &)

• somma: T& operator+=(T &)

• differenza: T& operator-=(T &)

17

• prodotto con scalare: T& operator*=(Real &)

• prodotto componente per componente: T& operator*=(T &)

• quoziente componente per componente: T& operator/=(T &).

Attributi

• d, w, a: vettori spostamento, velocita ed accelerazione del baricentro del corporigido. Le prime tre componenti (0,1,2) si riferiscono alle traslazioni rispetto agliassi (X,Y,Z rispettivamente), le seconde tre (3,4,5) alle rotazioni rispetto agli assi(X,Y,Z rispettivamente)

• inerzia: vettore contente la massa (componenti 0,1,2) e le inerzie relative ai treassi coordinati (componenti 3,4,5)

• newmarkSys: matrice del sistema di Newmark per il calcolo dell’accelerazione

• force: vettore delle forze e dei momenti agenti sul corpo rigido

• vincoli : variabile booleana che indica se il corrispettivo grado di liberta e vincolato(valore 0) oppure no (valore 1)

• time: istante attuale della simulazione.

Metodi

• Costruttore: inizializza gli attributi tramite Getpot. Non e fornito ne uncostruttore di default ne un costruttore di copia.

• Getters&Setters per la maggior parte degli attributi (per maggiori dettagli siveda la documentazione doxygen).

• void timeAdvance(const vettore & force, const Real & time): esegue unpasso di avanzamento temporale con lo schema di Newmark (tre fasi: predictor,eval, corrector).

• void predictor(): walkaround per DataNewmark::pred(vettore & d, vet-tore & w, const vettore & a).

• void eval(): risolve il sistema di Newmark per il calcolo dell’accelerazioneall’istante attuale.

• void corrector(): walkaround per DataNewmark::corr(vettore & d, vettore& w, const vettore & a).

18

2.2.3 La classe DataNewmark

La classe DataNewmark e una classe template virtuale che contiene al suo interno le strut-ture ed i metodi per l’implementazione del metodo di Newmark. Da tale classe possonoereditare sia un oggetto di tipo ElasticStructureSolver che uno di classe RigidBody.

Parametri Template

Tale classe ha tre parametri template che le garantiscono la compatibilita con la classefiglia: vettore, matrice, matloc. Ciascun parametro template T richiede che sianoimplementati i seguenti metodi:

• operatore di assegnazione: T & operator=(T &)

• somma: T& operator+=(T &)

• differenza: T& operator-=(T &)

• prodotto con scalare: T& operator*=(Real &).

Si richiede inoltre l’implementazione del prodotto matrice vettore. Tali parametri tem-plate verranno prossimamente spostati dalla definizione della classe a quella dei metodiche le utilizzano, per aumentare la flessibilita della classe.

Attributi

• dt : passo di avanzamento temporale

• endT : istante finale della simulazione

• beta, gamma: coefficienti dello schema di Newmark

• aRay, bRay : coefficienti di damping debole e forte (Rayleigh).

Metodi

• Costruttore: Inizializza gli attributi o tramite GetPot o direttamente. Non eimplementato il costruttore di default ne quello di copia.

• Getters&Setters: permettono di leggere e modificare il valore di tutti gli attributidella classe.

• void assembleLocalDamping(matloc &D, const matloc &M, const mat-loc &K): costruisce la matrice locale di damping alla Rayleigh utilizzando la for-mula D = aRay*M + bRay*K, essendo M e K le matrici locali rispettivamentedi massa e stiffness.

• void assembleLocalAcc(matloc &A, const matloc &M, const matloc&K, const matloc &D): costruisce la matrice locale del sistema di Newmarkutilizzando la formula A = M + gamma* dt*D+ beta* dt* dt*K.

19

• void computeRHS(vettore & F, const matrice & D, const matrice & K,const vettore & w, const vettore & d): calcola il termine noto del sistema diNewmark F = F - D*w - K*d. Si assume che F contenga il contributo delle forzedi volume e superficie.

• void pred(vettore & d, vettore & w, const vettore & a) calcola lo sposta-mento d e la velocita w predette all’istante attuale in funzione di spostamen-to d, velocita w ed accelerazione a all’istante precedente: d = d + dt*w +dt* dt*(.5- beta)*a; w = w + dt*(1- gamma)*a.

• virtual void eval(): metodo virtuale per la risoluzione del sistema di Newmarkper l’accelerazione. La classe figlia deve implementare tale metodo.

• void corr(vettore & d, vettore & w, const vettore &a) calcola lo sposta-mento d e la velocita w all’istante attuale in funzione di spostamento d e velocitaw predette ed accelerazione a all’istante attuale: d = d + dt* dt*( beta)*a; w =w + dt* gamma*a.

2.3 Classi per la gestione del contatto

Le classi per la gestione del contatto sono divise in due categorie.

• classi per la gestione del contatto tra un corpo elastico ed un corpo rigido: Con-tact, ContactFriction, FrictionSlip. Sono le prime ad essere state implementate,la geometria del contatto e altamente semplificata e ridotta ad un singolo pianoparallelo ad uno degli assi cartesiani. Sono state scritte pensando alla geometriaspecifica del contatto tra disco e pastiglia in un freno.

• classi per la gestione del contatto tra un corpo elastico e piu vincoli fissi: General-Contact, VincoloHD, Vincolo. Con tali classi si riesce ad affrontare problemi congeometrie piu complesse.

2.3.1 La classe Contact

Tale classe permette di integrare le equazioni del moto di un corpo elastico e di uncorpo rigido sotto il vincolo di non compenetrazione. Per semplicita nella descrizione del-la geometria del problema si assume che la superficie di contatto sia un piano ortogonalead uno degli assi cartesiani (X,Y,Z). E in via di sviluppo la possibilita di generaliz-zare il vincolo a piani generici introducendo una matrice di rotazione per costruire lecomponenti normali e tangenziali al vincolo.

Per massimizzare l’accuratezza, mantenere la consistenza, ridurre i costi com-putazionali e consigliabile utilizzare per il corpo elastico il metodo average acceleration(schema di Newmark con β = .25, γ = .5), mentre per il corpo rigido il metodo leap-frog(schema di Newmark con β = 0., γ = .5). Tale scelta garantisce un second’ordine diaccuratezza e buone proprieta di stabilita nell’avanzamento temporale.

20

Figura 2.4: Struttura delle classi implementate: linea blu ereditarieta, linea nerainclusione, linea rossa inclusione tramite puntatore

Un’ultima importante limitazione e che per calcolare in modo corretto la forza direazione applicata sul corpo rigido, occorre utilizzare una discretizzazione EF del corpoelastico che sia una partizione dell’unita (come ad esempio elementi finiti lineari sutetraedri). Tale limitazione e essenzialmente dovuta al fatto che, durante il calcolo, none direttamente disponibile la pressione di contatto tra i due corpi σn, bensı la forzalocale λi =

∫∂Ω σnγn(vi) essendo vi la funzione test associata al nodo di bordo i-esimo.

E immediato verificare che la forza scambiata tra i due corpi e pari a Σλi solo se lefunzioni test rappresentano una partizione dell’unita.

Parametri Template

Tale classe mantiene il parametro template di tipo Mesh della classe ElasticStructure-Solver.

Attributi

Riportiamo solo gli attributi necessari alla comprensione del funzionamento della classe;per la lista esaustiva si rimanda alla documentazione doxygen.

• rigid : puntatore ad un oggetto di classe RigidBody

• flexible: puntatore ad oggetto di classe ElasticStructureSolver

• offsetVC : direzione ortogonale al vincolo descritto dal corpo rigido (0 normale delvincolo orientata come X, 1 Y, 2 Z))

21

• g : vettore contenente il vincolo cinematico sul vettore spostamento

• vincolo: vettore contenete il vincolo cinematico sull’accelerazione

• cvSmooth: parametro di rilassamento per il vincolo di non penetrazione: 0 ilvincolo non viene imposto, 1 il vincolo e imposto in modo inviolabile

• lambda: vettore contenente i valori nodali della forza di contatto λi =∫∂Ω σnγn(vi)

• forceOnRigidBody : vettore risultante delle forze e momenti agenti sul corpo rigido

• boundaryDof : lista contenente le Id dei dof di bordo

• activeSet : vettore contenente le Id dei dof dove il vincolo di contatto e attivo. Eun sottoinsieme di boundaryDof

• A, APertubed : rispettivamente la matrice del sistema di Newmark per l’accele-razione libera del corpo elastico e quella modificata per tenere conto dei vincolicinematici

• rhs, rhsPertured : rispettivamente il termine noto del sistema di Newmark perl’accelerazione libera del corpo elastico e quello modificato per tenere conto deivincoli cinematici

• a: vettore soluzione contenente i valori nodali dell’accelerazione del corpo elastico.

Metodi

• Costruttore: riceve in ingresso una referenza ad un oggetto GetPot, ElasticStruc-tureSolver, RigidBody.

• void computedt(): calcola il passo temporale per lo step successivo in modo dadescrivere in modo accurato la fase di impatto e di allontanamento dei due corpi.Tale routine di adattazione temporale e particolarmente importante per garantirela conservazione dell’energia cinetica specialmente nel caso di corpi elastici adelevato modulo di Young. Il controllo avviene a livello dei predictor e assicurache la compenetrazione nella fase predictor sia inferiore ad una frazione prefissatadello spostamento totale, inoltre il tempo viene ulteriormente ridotto nella fasedi allontanamento dei due corpi. In tale modo e possibile descrivere in modomaggiormente accurato le forze impulsive che si generano tra i due corpi.

• void timeAdvance( source type const & source, const Real & time):gestisce l’avanzamento di uno step temporale per il sistema accoppiato dei duecorpi in contatto. Le fasi in cui esso avviene sono le seguenti:

– fase predictor : si chiamano i metodi predictor() degli oggetti rigid eflexible

– fase eval :

22

∗ preparazione della matrice e del termine noto per il corpo elastico nonvincolato

∗ valutazione del vincolo cinematico sull’accelerazione nodale del corpoelastico (computeVincolo())

∗ ricerca della configurazione che minimizza l’energia di deformazione elas-tica. Vengono effettuate le seguenti operazioni fino a convergenza (ilcriterio di stop e la mancata modifica tra due iterazioni dell’active set):· aggiunta dei vincoli e risoluzione del problema perturbato

(solvePerturbedProblem())· calcolo dei moltiplicatori di Lagrange lambda (computeLambda())· aggiornamento dell’active set in base all’ammissibilita e all’ottimalita

della soluzione all’iterazione corrente∗ calcolo della forza e del momento risultante sul corpo rigido

(computeForceOnRigidBody())

– fase corrector : si chiamano i metodi predictor() degli oggetti rigid eflexible.

• void computeVincolo(): calcola il vincolo cinematico sullo spostamento g esull’accelerazione vincolo mediante la relazione (1.27).

• void solvePerturbedProblem(): APerturbed e rhsPerturbed vengono calco-lati mediante il metodo MSRMatr::diagonalize(dofID, tvg, rhs, valore) cheannulla le componenti extra-diagonali della riga e colonna dofID, pone uguale atvg l’elemento diagonale (dofID,dofID) e modifica il termine noto rhs in modoche la soluzione nel nodo dofID valga valore. Infine si risolve il sistema algebricomediante il solutore Atzec.

• void computeLambda(): calcola il valore dei moltiplicatori di Lagrange lambdacome il residuo della soluzione vincolata a rispetto al problema non perturbato.lambda = A* a - rhs. I lambda assumono il significato fisico delle forze esercitatesu ciascun nodo di bordo, e per risalire alla pressione di contatto σn andrebbe risoltoun ulteriore sistema algebrico:

Mi,j(σn)j = λi dove Mi,j =∫

∂Ωγn(vj)γn(vi)

• void refreshActiveSet(): aggiorna l’active set tenendo conto dell’ammissibilita(rispetto dei vincoli) e dell’ottimalita (segno dei moltiplicatori) della soluzionecome spiegato nell’algoritmo (1.1) al punto a.

• void computeForceOnRigidBody(): calcola forze e momenti sul corpo rigido(principio di azione e reazione). Nel calcolo si sfrutta il fatto che gli EF utilizzatisiano una partizione dell’unita.

23

2.3.2 La classe ContactFriction

Questa classe generalizza la classe Contact (da cui eredita) aggiungendo le routine peril calcolo dell’attrito dinamico e quasi statico.

Parametri Template

Tale classe mantiene il parametro template di tipo Mesh della classe Contact da cuieredita.

Attributi

• mu: oggetto di classe FrictionSlip. Si comporta come un funtore e restituisce ilvalore del coefficiente d’attrito in funzione della velocita di slip tra i due corpi

• friction: vettore contente la forza d’attrito tra i due corpi

• radiusYZ : vettore contenente la distanza dall’asse X di ciascun grado di libertadel corpo deformabile, serve per calcolare i momenti agenti sul corpo rigido.

Metodi

• Costruttore: riceve in ingresso una referenza ad un oggetto GetPot, ElasticStruc-tureSolver, RigidBody.

• void timeAdvance( source type const & source, const Real & time):gestisce l’avanzamento di uno step temporale per il sistema accoppiato dei duecorpi in contatto. Le uniche modiche nel metodo rispetto alla classe base sta nellafase eval: si utilizza un metodo iterativo nel quale prima si calcola la soluzioneal problema del solo contatto utilizzando la forza d’attrito calcolata all’iterataprecedente come forzante e successivamente si calcola la forza d’attrito all’itera-ta corrente mediante lo sforzo di contatto appena calcolato. L’algoritmo termi-na quando viene soddisfatta una certa tolleranza sulla norma della differenza delvettore accelerazione a tra due iterate successive.

• void computeFriction(): lo sforzo di attrito dinamico e pseudo-statico ff vienecalcolato utilizzando una generalizzazione della legge di Coulomb:

ff = µ(‖vslip‖)σnvslip

‖vslip‖dove il coefficiente d’attrito µ dipende in modo non lineare dalla velocita di slidingtra i due corpi (cfr class SlipFriction). Poiche non si ha a disposizione direttamenteσn, ma solo la sua valutazione contro le funzioni test v, lo sforzo d’attrito e calcolatoin modo consistente solo nel caso di EF lineari a pezzi e con un errore pari all’erroredi quadratura con tecnica del mass lumping (metodo dei trapezi).

24

• void computeForceOnRigidBody(): calcola forze e momenti sul corpo rigido(principio di azione e reazione). Nel calcolo si sfrutta il fatto che gli EF utilizzatisiano una partizione dell’unita.

2.3.3 La classe SlipFriction

Tale classe descrive la curva coefficiente d’attrito in funzione della velocita di slip. Eimplementata come un funtore.

−15 −10 −5 0 5 10 15−0.5

−0.4

−0.3

−0.2

−0.1

0

0.1

0.2

0.3

0.4

0.5

vslip

µ

Figura 2.5: Curva coefficiente d’attrito versus velocita di slip

Attributi

• muDin: coefficiente d’attrito dinamico

• muMax : coefficiente d’attrito quasi statico, si ottiene per bassi valori di velocitadi sliding

• vSlipDin: valore della velocita di sliding oltre la quale si ha il coefficiente d’attritodinamico

• vSlipMax : valore della velocita di sliding alla quale si ottiene il massimo delcoefficiente d’attrito.

Metodi

• Costruttore: inizializza gli attributi da GetPot.

• Real operator(const Real & slipvel): valuta la curva coefficiente d’attrito infunzione della velocita sliding nel punto slipvel.

25

2.3.4 La classe GeneralContact

Tale classe e in grado di gestire il contatto tra un corpo deformabile e piu vincoli descritticome piani perpendicolari agli assi cartesiani. L’algoritmo di contatto e analogo a quelloimplementato dalla classe Contact, le strutture dati sono state invece modificate in mododa poter tener conto della presenza di un numero maggiore di vincoli.

Attributi

• constrain: puntatore ad un oggetto di classe VincoloHD, contiene le informazionirelative ai vincoli

• flexible: puntatore ad un oggetto di classe ElasticStructureSolver, contiene le infor-mazioni relative al corpo deformabile. La discretizzazione spaziale deve essere unapartizione dell’unita e l’avanzamento in tempo deve avvenire mediante lo schemadi Newmark

• lambda: vettore dei moltiplicatori di Lagrange associati al problema discreto,corrisponde alla valutazione dello sforzo normale di contatto σn contro le funzionitest vi

• A, APerturbed : matrice del sistema di Newmark libero e vincolato

• rhs, rhsPerturbed : termine noto del sistema di Newmark libero e vincolato

• a: vettore soluzione del sistema di Newmark vincolato

• boundaryDof : lista dei gradi di liberta di bordo, ossia i possibili punti di contatto.

Metodi

• Costruttore: inizializza gli attributi.

• void buildBoundaryDof(): costruisce la lista dei nodi di bordo. Tale funzionepresuppone che la numerazione dei nodi di griglia disponga per primi i nodi dibordo, come nel caso di mesh generate con NetGen.

• void timeAdavance(source type const & source, const Real & time):esegue uno step temporale risolvendo il problema di ottimizzazione vincolata. Cal-cola spostamenti, velocita, accelerazioni, sforzi di contatto nel corpo deformabile ele reazioni vincolari globali.

• void solvePerturbedProblem(): aggiunge i vincoli al sistema lineare e lo risolvetramite la libreria Atzec.

• void computeLambda(): calcola il valore dei moltiplicatori di Lagrange discreticome residuo della soluzione vincolata rispetto al problema non vincolato.

26

• void computeGlobalForce(): calcola le reazioni vincolari globali. Questa fun-zione presuppone che gli EF utilizzati per discretizzazione spaziale del corpodeformabile siano una partizione dell’unita.

2.3.5 La classe VincoloHD

Tale classe e un contenitore per la classe Vincolo. Il suo attributo principale e unostd::vector¡Vincolo¿, ed i metodi in essa implementati permettono di ciclare sui metodidella classe Vincolo.

Attributi

• vlist : standard vector contenente oggetti di classe Vincolo

• size: dimensione del vettore vlist

• nActiveConstrain: numero di gradi di liberta del corpo deformabile vincolati

• changedActiveSet : variabile booleana che indica se l’active set e stato modificatotra due iterazioni successive. Rappresenta il criterio d’arresto del ciclo iterativo.

Metodi

• Costruttore: inizializza una lista vuota.

• Setters&Getters per i principali attributi.

• void addVincolo(const Vincolo & a): aggiunge un oggetto di classe Vincoloa vlist.

• void computeVincolo(const BODY & flex, const vettore & xx): richiamal’omonima funzione di ciascun oggetto in vlist.

• UInt & refreshActiveSet(const vettore & lambda, const vettore & a,const std::list¡UInt¿ & boundaryDof): chiama l’omonima funzione di ciascunoggetto di vlist, conteggia i gradi di liberta vincolati, ricerca se sia stato modificatol’active set.

• void perturbedProblem(matrice & A, vettore & rhs): chiama l’omonimafunzione di ciascun oggetto di vlist.

2.3.6 La classe Vincolo

La classe Vincolo contiene la descrizione geometrica del vincolo (normale al piano,termine noto) ed i metodi per calcolare e imporre i vincoli cinematici.

27

Attributi

• cComp: direzione normale al piano del vincolo: 0 piano normale all’asse X, 1 asseY, 2 asse Z

• offsetVC : offset per raggiungere la componente cComp in un PhysVecUnknown;viene inizializzato in automatico dal costruttore

• dim: numero di elementi nella mesh del corpo deformabile; fornisce la dimensionedei vettori per il vincolo cinematico

• value: termine noto del piano di vincolo

• direction: indica il verso concorde all’asse (+1) o opposto (-1) della direzionenormale ammessa dal vincolo

• cvSmooth: parametro di rilassamento del vincolo di non compenetrazione

• penality : costante numerica di penalizzazione per la violazione del vincolo nelcalcolo dell’ActiveSet

• tgv : valore numerico da imporre agli elementi diagonali della matrice del sistemalineare nella procedura di esplicitazione dei vincoli

• g : vincolo cinematico sugli spostamenti; e un ScalUnknown di dimensione dim

• vincolo: vincolo cinematico sull’accelerazione; e un ScalUnknown di dimensionedim

• nActiveConstrain: numero di gradi di liberta del corpo deformabile vincolati

• changedActiveSet : variabile booleana che tiene traccia del fatto che l’ActiveSetsia stato modificato tra due chiamate della funzione UInt refreshActiveSet(...)

• name: nome del vincolo, e usato da GetPot nella lettura delle informazioni dafile.

Metodi

• Costruttori: e fornito un costruttore di default, un costruttore di copia, ed uncostruttore con lettura delle informazioni da GetPot.

• void computeVincolo(const BODY & flex, const vettore & xx): calcolai vincoli cinematici su spostamenti ed accelerazioni. La classe template BODYdeve ereditare da DataNewmark ed implementare il metodo vettore & disp()che ritorna il vettore spostamento del corpo deformabile. Il vettore xx contiene lecoordinate dei nodi della mesh.

28

• bool isViolated(const vettore & lambda, const vettore & a, const UInt& i): implementa il controllo sul grado di liberta i-esimo. In particolare tale gradodi liberta viene incluso nell’ActiveSet se

λi + c(vincoloi − ai) < 0

• UInt refreshActiveSet(const vettore & lambda, const vettore & a, conststd::list¡UInt¿ & boundarydof): aggiorna l’ActiveSet e ritorna il numero digradi di liberta vincolati.

• void perturbedProblem(matrice & A, vettore & rhs): modifica la matriceA tramite la funzione void MSRMatr::diagonalize(UInt comp, Real tgv,vettore rhs, Real value).

2.4 Classi di supporto

2.4.1 Class AnaliticalFunction

Tale classe ha lo scopo di fornire un metodo di interpolazione sulla mesh di calcolo perfunzioni analitiche. La classe AnaliticalFunction e una classe virtuale, priva di attributi,dotata dei metodi per l’interpolazione. Dichiara inoltre l’operatore virtuale operator()che verra implementato nelle classi figlie.

Parametri Template

La classe utilizza tre diversi parametri template:

• Mesh: rappresenta la classe a cui l’oggetto che contiene la mesh deve appartenere.

• vettore: solitamente le classi figlie assegnano a questo parametro la classeLifeV::ScalUnknow nel caso di grandezza scalare e LifeV::PhysUnknownVect nelcaso di grandezza vettoriale. Tale parametro template richiede che sia imple-mentato l’accesso ai suoi elementi tramite l’operatore Real & operator()(UInti).

• UInt NCOMP : numero di componenti della funzione: 1 scalare, 3 vettoriale;solitamente la classe figlia specifica il valore del parametro template.

Tali parametri template verranno prossimamente spostati dalla definizione della classea quella dei metodi che le utilizzano, per aumentare la flessibilita della classe. Il valoredi NCOMP sara passato ai vari metodi come parametro di input.

Metodi

• Costruttore di default: do-nothing constructor.

29

• virtual Real operator()(const Real &t, const Real &x, const Real &y,const Real &z, UInt component): valuta la funzione all’istante t, nel punto(x,y,z); nel caso di funzione vettoriale component specifica la componente che vieneritornata: 1 X, 2 Y, 3 Z.

• void interp(const Mesh & mesh, const RefFe & refFe, const GeoMap &geoMap, const QuadRule & Qr, const Dof & dof, vettore & vect, Realt) dove i primi cinque parametri servono a determinare le coordinate dei nodi doveoccorre interpolare, vect e il vettore che accoglie i valori interpolati e t e l’istantetemporale a cui valutare la funzione nel caso in cui essa presenti una dipendenzatemporale.

Classi figlie

Le classi figlie implementano alcune funzioni analitiche tra cui:

• identita vettoriale f(x) = x

• estrazione di una componente scalare: fi(x) = xi

• distanza dall’asse xi: f(x) =√

x2j + x2

k.

30

Capitolo 3

Simulazioni

3.1 Sfera su basamento fisso

Il primo caso test per il codice implementato consiste nella simulazione di una sferaviscoelastica che rimbalza su un basamento rigido. Per la sfera sono stati consideratidue diversi livelli di rigidezza:

• modulo di Young pari a .05 MPa (bassa rigidezza) per poter apprezzare le variazionidi forma del corpo elastico a causa dell’impatto con il basamento fisso;

• modulo di Young pari a 5 MPa (alta rigidezza) per valutare la stabilita del codicenel caso di elevate forze impulsive.

La seconda simulazione, in particolare, ha evidenziato l’importanza nell’adattivita intempo, soprattutto nella fase di rilascio del contatto, al fine di garantire il principio diconservazione dell’energia cinetica: senza un opportuno raffinamento del passo di discre-tizzazione temporale si assiste infatti ad un incremento dell’energia cinetica dopo l’urto,dovuto essenzialmente allo scarso dettaglio con cui si individua l’istante di separazionetra i due corpi in contatto.

I risultati riportati riguardano la seguente simulazione:

modulo Young .05÷ 5 [MPa]coefficiente Poisson 0.28 [·]densita 2.8·10−6 [ kg

mm3 ]raggio sfera 10 [mm]Coeff. Rayleigh a 0. [s−1]Coeff. Rayleigh b .01 [s]

In seguito sono riportati i grafici che mostrano la legge oraria della componente xdel baricentro, le forze scambiate tra i due corpi ed il passo di integrazione temporale.Si mostrano infine lo stato di deformazione e lo sforzo di contatto massimo.

31

0 2 4 6 8 10 12−25

−20

−15

−10

−5

0

sphereCM (5MPa)sphereCM (.05MPa)fondationunloaded radius

0 2 4 6 8 10 120

1

2

3

4

5

6

7

8

9

5MPa0.05MPa

Figura 3.1: A sinistra posizione del baricentro della sfera; a destra forze scambiate

Figura 3.2: Stato di sforzo massimo: E=5MPa a sinistra, E=.05MPa a destra

32

3.2 Sfera su basamento mobile

In questa simulazione si utilizza l’algoritmo presentato nella sezione 1.2.5. Si considera lasfera viscoelastica del test precedente (E=0.05MPa) poggiata su una molla precaricata.Si mostra sia la traiettoria del centro di massa del corpo rigido e della sfera deformabilesia le forze scambiate.

0 1 2 3 4 5 6 7 8 9 100

0.5

1

1.5Forze scambiate

0 1 2 3 4 5 6 7 8 9 10−35

−30

−25

−20

−15

−10

−5

0

5

10CM position

sphereCMrigidCMunloaded radius

Figura 3.3: A sinistra: forze scambiate; a destra: posizione del baricentro del corporigido (rosso), della sfera deformabile (blu). La distanza tra la linea gialla e quella blurappresenta la compressione del raggio della sfera

−40−20

0

−20

0

20−20

0

20

3d view; time =0.05 s

−40 −30 −20 −10 0 10−20

−10

0

10

20xy view; time =0.05 s

−40 −30 −20 −10 0 10−20

−10

0

10

20xz view; time =0.05 s

−20 −10 0 10 20−20

−10

0

10

20yz view; time =0.05 s

Figura 3.4: Animazione

33

3.3 Sfera in una scatola

In questa simulazione si utilizzano le classi GeneralContact, VincoloHD, Vincolo. Siconsidera una sfera di diametro unitario a bassissimo modulo di Young soggetta allaforza di gravita e vincolata a rimanere in una scatola di base quadrata (10 unita dilato) e aperta sulla faccia superiore. Si mostrano le coordinate del baricentro e le forzescambiate contro le pareti della scatola.

0 1 2 3 4 5 6 7 8 9 10−5

−4

−3

−2

−1

0

1

2

3

4

5

xyz

0 1 2 3 4 5 6 7 8 9 10−1.5

−1

−0.5

0

0.5

1

1.5x 10

−3

xyz

Figura 3.5: A sinistra: coordinate del baricentro della sfera; a destra: forze scambiatetra la sfera e le pareti della scatola

−5

0

5

−5

0

50

2

4

6

3d view

−5 0 5−5

0

5xy view

−5 0 5

0

1

2

3

4

5

6

xz view

−5 0 5

0

1

2

3

4

5

6

yz view

Figura 3.6: Animazione

34

3.4 Attrito disco-pastiglia

Descrizione di una prova al banco

Mostriamo infine i risultati di una simulazione per una prova al banco relativa al-l’impianto frenante anteriore della Ducati 999; l’impianto di questa moto e compostoda 2 dischi da 320 mm di diametro semiflottanti, 2 pinze, 2 pastiglie per ogni pinza,ognuna spinta da 2 pistoni idraulici.La prova consiste nel montare l’impianto su un volano che rappresenta l’inerzia equiva-lente della moto durante la frenata misurando poi la decelerazione imposta alla motoin funzione della pressione dell’impianto. Per simulare correttamente la decelerazione diun motoveicolo in una prova al banco, e necessario studiare i trasferimenti di carico infase di frenata e modificare l’inerzia del volano in funzione del maggior carico che gravasull’anteriore.

Considerando la Ducati 999 si ha

• Massa totale veicolo [kg]: Mtot=270

• Massa anteriore statica [kg]: Mant=125

• Altezza baricentro [m]: h=0.625

• Passo [m]: p=1.42

• Raggio rotolamento pneumatico[m]: Rrot=0.309

• Decelerazione lineare baricentro del motociclo[

ms2

]: a

Imponendo l’equilibrio a rotazione rispetto al punto di contatto strada-pneumaticoposteriore, otteniamo il trasferimento di carico dinamico in frenata

Ntraf = Mtoth

pa [N ], (3.1)

mediante il quale e possibile calcolare la massa dinamica ridotta all’anteriore:

Mrid−ant = Mant + Mtoth

p

a

g[kg] (3.2)

Da quest’ultima, infine, nell’ipotesi di contatto ideale tra pneumatico e strada (i.e.in assenza di scivolamento) si ottiene il momento d’inerzia equivalente per il volano:

J = Mrid−antR2rot [kgm2] (3.3)

35

Dati fisici e geometrici della pastiglia

La pastiglia e assimilata ad un corpo viscoelastico lineare. La mesh di calcolo e statagenerata con NetGet, conservando, nonostante alcune semplificazioni di forma, le carat-teristiche dell’oggetto reale. E inoltre stata implementata una routine di calcolo chepermette di apportare lievi modifiche alla geometria della mesh al fine di poter valutarel’effetto di caratteristiche come la rugosita o l’usura del materiale d’attrito.

Le caratteristiche fisico-geometriche della pastiglia sono le seguenti:

Modulo di Young 502 [MPa]Coeff di Poisson 0.20 [-]Densita 2.43e-6 [kg/m3]spessore 5 [mm]area di contatto 22 [cm2]

Risultati

Nel seguito riportiamo sia alcuni risultati di carattere globale - quali velocita e deceler-azione del motociclo, forza e di contatto e coppia d’attrito tra disco e pastiglia - sia alcunirisultati di carattere piu fine, come lo stato di sollecitazione all’interno della pastiglia.Dai grafici e possibile notare una certa difficolta del codice di descrivere la fase di primocontatto tra il disco e pastiglia (tempo: 0.5-0.7 sec). Si assiste infatti a grandi oscillazioninella forza di contatto tra pastiglia e disco. L’infittimento del passo temporale in talefase permette di ottenere una descrizione smooth della velocita del motociclo, nonostantela coppia frenante applicata presenti picchi elevati. Non e stata ancora perfettamentechiarita la natura delle oscillazioni del baricentro della pastiglia durante la fase a regimedella frenata. Vi sono in particolare alcune riserve sulla fisicita di tale comportamento:queste frequenze di vibrazione potrebbero essere dovute alla discretizzazione spazialeanziche alle proprieta fisiche del materiale. Si noti, inoltre, (tempo: 3.5 s) come al fer-marsi del motociclo anche la coppia frenante si annulli. Il picco della coppia frenantea basse velocita e facilmente interpretabile considerando la caratterizzazione del coeffi-ciente d’attrito in funzione della velocita relativa tra i due corpi in contatto. Infine, nelprofilo degli forzi di contatto scambiati, e possibile notare con chiarezza l’impronta deipistoni che premono la pastiglia contro il disco.

36

0 0.5 1 1.5 2 2.5 3 3.5 40

0.2

0.4

0.6

0.8

1

time [s]

pres

sure

[MP

a]

0 0.5 1 1.5 2 2.5 3 3.5 4−500

0

500

1000

1500

2000

time [s]

forc

e N

; tor

que

[Nm

]

normal forcefriction torque

0 0.5 1 1.5 2 2.5 3 3.5 4−0.05

0

0.05

0.1

0.15

0.2

0.25

0.3

time [s]

linin

g ce

nter

of m

ass

osci

llatio

n [m

m]

Figura 3.7: A sinistra: pressione in ingresso - condizione al bordo esogena al modello; alcentro:forze scambiate; a destra: posizione baricentro della pastiglia.

0 0.5 1 1.5 2 2.5 3 3.5 4−7

−6

−5

−4

−3

−2

−1

0

1

time [s]

acce

lera

tion

[m/s

2]

0 0.5 1 1.5 2 2.5 3 3.5 4

0

2

4

6

8

10

12

14

16

time [s]

velo

city

[m/s

]

Figura 3.8: Decelerazione (a sinistra) e velocita (a destra) del motociclo

Figura 3.9: Sforzo di contatto normale sulla pastiglia

37

Conclusioni

L’algoritmo di contatto, implementato dapprima facendo riferimento alla geometria delproblema specifico del sistema pastiglia-disco di un freno, e stato successivamente ge-neralizzato a geometrie piu complesse contenenti piu vincoli contemporaneamente. Irisultati delle simulazioni in linea con la fisica del problema ed i tempi di calcolo accetta-bili hanno dimostrato la bonta di un approccio variazionale basato sulla minimizzazionedel funzionale energia di deformazione nel trattare problemi di contatto in presenza omeno di attrito.

Si ottiene infatti una descrizione assai curata del moto vincolato di un corpo ela-stico, pur in presenza di oscillazioni numeriche per l’andamento della forza di contatto.L’adattivita in tempo ha lo scopo di attenuare tali oscillazioni. Essa e particolarmenteimportante nel caso di corpi elastici ad elevata rigidezza, soggetti per lo piu a forze im-pulsive nella ridotta frazione di tempo in cui sono a contatto con il vincolo. Utilizzandoun passo temporale superiore al tempo di contatto tra il corpo elastico ed il vincolosi introduce un grave errore sulla conservazione dell’energia cinetica, che tende ad au-mentare dopo l’urto. L’adattazione in tempo infittisce dunque il time-step nella fase diriallontanamento del corpo con il vincolo in modo da individuare con esattezza l’istantedi fine contatto. Va comunque notato che l’infittimento del passo di discretizzazionetemporale, pur migliorando la descrizione cinematica dello spostamento e velocita deicorpi ed il rispetto del principio di conservazione dell’energia cinetica, amplifica le oscil-lazioni della forza di contatto scambiata. Uno dei prossimi obiettivi sara dunque quellodi capire meglio questo trade-off e progettare metodi efficienti per lo smorzamento delleoscillazioni.

Sviluppi futuri del codice riguardano sia la parte modellistica che quella algoritmi-ca. Per quanto riguarda il primo aspetto, si vuol considerare l’effetto della rugositasuperficiale dei corpi a contatto (gia implementato un codice in grado di aggiungere unaperturbazione gaussiana nei nodi di bordo per la mesh) ed introdurre semplici modellidi usura. Per rimuovere la limitazione sulla direzione del vincolo, si vuole inoltre intro-durre una matrice di rotazione per esplicitare nel sistema lineare le componenti normalie tangenziali al vincolo. Per quanto riguarda l’aspetto algoritmico, sara implementato abreve un metodo di Newmark modificato a dissipazione controllata capace di tagliare leoscillazioni numeriche alle alte frequenze proprie del metodo average acceleration.

38

Descrizione eseguibili

Per compilare i sorgenti occorre disporre di un playground per LifeV correttamenteconfigurato. Con il comando make check vengono compilati tutti i makefile presentinella directory uvilla e nelle relative sottodirectory. Vengono inoltre creati dei linksimbolici ai file di input e di postprocess. Segue la descrizione degli eseguibili disponili.

• anis mesh: crea una nuova mesh anisotropa deformando una mesh isotropagenerata da NetGen.

• body interaction: risolve il caso test della sfera che rimbalza su un fondamentoelastico. Postprocess 3d (grafica.m) e 0d (elastfond.m) in Matlab, visualizzazionesforzi di contatto in Paraview.

• contact: risolve il caso test della sfera su fondamento rigido. Postprocess 0d(fixfond.m) in Matlab, visualizzazione sforzi di contatto in Paraview.

• gcontact: risolve il caso test della sfera in una scatola. Postprocess 3d (sfera3d.m)e 0d (baricentro.m) in Matlab, visualizzazione sforzi di contatto in Paraview.

• lap mesh: crea una nuova mesh introducendo una perturbazione aleatoria gaus-siana sulla superficie di contatto della pastiglia. Risolve un problema di laplacianoper determinare lo spostamento dei nodi interni.

• newmark: risolve un problema di elasticita lineare con il metodo di Newmark.

• pastigliadisco: risolve il problema del contatto con attrito tra pastiglia e dis-co. Postprocess 0d (postproc.m) in Matlab, visualizzazione sforzi di contatto inParaview.

39

Bibliografia

[1] L. FORMAGGIA, J. GERBEAU, and C. PRUD’HOMME. LifeV Developer Manual,Mag 2004. www.lifev.org.

[2] F. HECHT, O. PIRONNEAU, A. LE HYARIC, and K. OHTSUKA. FreeFem++,Nov 2006.

[3] T.J.R. HUGHES. The finire element method. Linear Static and Dynamic finiteelement analysis. Prentice-Hall, 1987.

[4] N. KIKUCHI and J. T. ODEN. Contact Problems in Elasticity: A Study of Vari-ational Inequalities and Finite Element Methods. Siam, Philadelphia, 1988. ISBN0-89871-202-5.

[5] A. QUARTERONI and A. VALLI. Numerical Approximation of Partial DifferentalEquations. Springer-Verlag, 1994.

[6] A. SCOTTI. Un approccio variazionale per il vincolo di galleggiamento. Master’sthesis, Politecnico di Milano, 2004.

[7] S. TUMINARO, M. HEROUX, A. HUTCHINSON, and N. SHADID. Official AztecUser’s Guide. v2.1. Sandia National Laboratories, Nov 1999.

[8] D. YANG. C++ and object-oriented numeric computing. Springer, New York, 2001.

40