Metodi numerici per lo studio di sistemi...

37
Universit ` a di Bergamo – Facolt ` a di Ingegneria – Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi , R. Riva , M. Camposaragna , R. Strada , B. Zappa Dipartimento di Progettazione e Tecnologie – Universit` a di Bergamo Dipartimento di Elettrotecnica – Politecnico di Milano

Transcript of Metodi numerici per lo studio di sistemi...

Universita di Bergamo

– Facolta di Ingegneria –

Metodi numerici per lo studio disistemi multicorpo

V. Lorenzi†, R. Riva†, M. Camposaragna‡, R. Strada†, B. Zappa†

† Dipartimento di Progettazione e Tecnologie – Universita di Bergamo‡ Dipartimento di Elettrotecnica – Politecnico di Milano

INDICE I

Indice

1 Introduzione 1

2 Concetti di base 22.1 Sistema multicorpo e vincoli . . . . . . . . . . . . . . . . . . . . . . . . . 22.2 Sistemi di coordinate dipendenti e indipendenti . . . . . . . . . . . . . . 32.3 Formulazione simbolica–numerica . . . . . . . . . . . . . . . . . . . . . . 42.4 Classi di problemi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

2.4.1 Analisi cinematica . . . . . . . . . . . . . . . . . . . . . . . . . . 42.4.2 Analisi dinamica . . . . . . . . . . . . . . . . . . . . . . . . . . . 52.4.3 Sintesi . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3 Coordinate ed equazioni di vincolo 63.1 Sistemi Multibody piani . . . . . . . . . . . . . . . . . . . . . . . . . . . 6

3.1.1 Coordinate relative . . . . . . . . . . . . . . . . . . . . . . . . . . 63.1.2 Reference point coordinate . . . . . . . . . . . . . . . . . . . . . 83.1.3 Coordinate naturali . . . . . . . . . . . . . . . . . . . . . . . . . 9

3.2 Sistemi multicorpo spaziali . . . . . . . . . . . . . . . . . . . . . . . . . 103.2.1 Coordinate relative . . . . . . . . . . . . . . . . . . . . . . . . . . 103.2.2 Reference point coordinate . . . . . . . . . . . . . . . . . . . . . 113.2.3 Coordinate naturali . . . . . . . . . . . . . . . . . . . . . . . . . 133.2.4 Vincoli anolonomi . . . . . . . . . . . . . . . . . . . . . . . . . . 13

4 Analisi cinematica 144.1 Assemblaggio . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144.2 Velocita . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 174.3 Accelerazione . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 184.4 Spostamenti finiti e simulazione cinematica . . . . . . . . . . . . . . . . 184.5 Sistemi con vincoli sovrabbondanti . . . . . . . . . . . . . . . . . . . . . 18

5 Analisi dinamica 195.1 Equazioni di moto e dinamica diretta . . . . . . . . . . . . . . . . . . . . 195.2 Posizione di equilibrio statico . . . . . . . . . . . . . . . . . . . . . . . . 225.3 Dinamica inversa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 235.4 Linearizzazione delle equazioni di moto: piccole oscillazioni attorno alla

posizione di equilibrio statico . . . . . . . . . . . . . . . . . . . . . . . . 24

6 Integrazione delle equazioni di moto 246.1 Concetti di base . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 256.2 Metodi tipo Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . 266.3 Metodi multistep . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 276.4 Metodo di Newmark . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 286.5 Integrazione diretta delle DAE . . . . . . . . . . . . . . . . . . . . . . . 28

7 Sintesi cinematica e dinamica 29

8 Elementi flessibili 31

9 Conclusioni 32

10 Bibliografia 33

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

1 INTRODUZIONE 1

1 Introduzione

L’analisi cinematica e dinamica di sistemi multicorpo costituisce una parte importantedel progetto delle macchine e rientra quindi tra gli strumenti che il CAD–CAE rendedisponibili al progettista.I sistemi meccanici che sono compresi nella definizione di sistemi multibody includono irobot, i veicoli spaziali, gli autoveicoli con i diversi sottosistemi che li compongono (so-spensioni, meccanismo di sterzo, trasmissione), le macchine tessili, le macchine utensili(figura 1).Si tratta in sostanza di un insieme di corpi, rigidi e flessibili, tra loro interconnessi inmodo da presentare possibilita di moti relativi.

Figura 1.1: Esempi di sistemi multibody

Normalmente gli elementi che compongono questi sistemi sono soggetti a grandi spo-stamenti e questo comporta che la configurazione e le condizioni al contorno dei sistemimulticorpo siano soggette a notevoli variazioni, contrariamente a quanto accade in ge-nere per problemi prettamente strutturali che vengono affrontati con il metodo deglielementi finiti. Inoltre le cadenze a cui operano le macchine moderne sono spesso elevatee comportano notevoli accelerazioni e quindi forze di inerzia: questo fa inevitabilmenteinsorgere problemi di tipo dinamico, di cui si deve tenere conto nella fase di progetto.La possibilita di simulare questi sistemi, offerta talvolta all’interno degli stessi pro-grammi CAD o da programmi dedicati che si interfacciano facilmente con questi ulti-mi, consente di predire il comportamento cinematico e dinamico di un generico sistema

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

2 2 CONCETTI DI BASE

multicorpo, ad ogni fase della progettazione, e di analizzare l’influenza dei differentiparametri di progetto in modo rapido ed economico.Mentre pero sono disponibili numerosi programmi per l’analisi di sistemi multicorpouna volta che siano state definite tutte le caratteristiche geometriche e dinamiche, almomento invece non si puo dire che esista uno strumento “general purpose” per la sintesi(o progettazione) “automatica” di questi sistemi, che consenta quindi l’ottimizzazionedei parametri di progetto rispettando i vincoli imposti.Sono pero disponibili strumenti che consentono di effettuare analisi di tipo parametrico(di sensitivita) che determinano la risposta di un sistema multicorpo al variare di unoo piu parametri di progetto. In ogni caso un programma di analisi costituisce la basedel processo di ottimizzazione o dell’analisi di sensitivita.I programmi multibody, proprio perche spesso il progetto si basa su analisi ripetute,sono utilizzati in modo interattivo anche nella fase di soluzione, a differenza di quan-to succede col FEM: l’operatore puo intervenire durante la simulazione modificandocarichi, vincoli, controlli.Inoltre le equazioni dinamiche possono essere sia utilizzate per progettare il control-lore del sistema che inserite (eventualmente in una forma semplificata) nell’anello dicontrollo real-time del sistema stesso. Per queste applicazioni la velocita di soluzione,importante comunque per ogni tipo di analisi, risulta cruciale.

2 Concetti di base

2.1 Sistema multicorpo e vincoli

Un sistema multicorpo e definito come un insieme di due o piu corpi collegati tra loroin modo che sia conservata la possibilita di moto relativo. Il collegamento, se dovutoall’accoppiamento cinematico di superfici, linee o punti, viene indicato col nome digiunto. Un giunto consente alcuni movimenti relativi e ne impedisce o limita altri.Alcuni vincoli lasciano un solo grado di liberta relativo, altri due e cosı via, fino ai seilasciati dal vincolo “nullo”. Un giunto rotoidale (cerniera), ad esempio, lascia un sologrado di liberta relativo (la rotazione attorno all’asse della cerniera), come il giuntoprismatico e quello elicoidale. Altri giunti spesso utilizzati sono quello cilindrico ecardanico (2 gdl, figura 2.1 e figura 2.2) e sferico (3 gdl, figura 2.3).

Figura 2.1: Giunto cardanico

I gradi di liberta del sistema corrispondono al numero di coordinate indipendenti chesono necessarie a definire la posizione del sistema stesso.I collegamenti tra due membri a volte avvengono non per un accoppiamento cinematico(di forma), ma attraverso elementi come molle, smorzatori etc. che non diminuiscono igradi di liberta del sistema, o mediante accoppiamenti di forza (cioe con contatti non

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

2 CONCETTI DI BASE 3

Figura 2.2: Giunto cilindrico

Figura 2.3: Giunto sferico

bilateri) che possono dare luogo a variazioni nel numero di gradi di liberta nel sistemaa seconda dei carichi agenti. Vincoli particolari, detti di mobilita o anolonomi, limitanosemplicemente le possibilita di movimento senza ridurre i gradi di liberta del sistema.Un sistema multicorpo puo essere a catena chiusa o aperta, a seconda che presenti omeno anelli chiusi: un bipendolo, o piu in generale sistemi ad albero sono esempi dicatene aperte, un quadrilatero articolato costituisce un esempio di sistema a catenachiusa con un solo anello.

2.2 Sistemi di coordinate dipendenti e indipendenti

Il primo passo per descrivere un sistema multicorpo consiste nel definire un set dicoordinate che consenta di individuare in modo univoco, e ad ogni istante, la posizionedi tutti i corpi che compongono il sistema. Diversi sono gli approcci possibili e ingenerale non si puo dire che esista un approccio in assoluto migliore di altri.Una prima scelta puo essere quella di un sistema di coordinate indipendenti, il cuinumero coincide con il numero di gradi di liberta del sistema, ed e percio minimo.Un inconveniente di queste coordinate e che definiscono direttamente la posizione dialcuni corpi, mentre per determinare la posizione di altri occorre procedere ad ulteriorianalisi (come ad esempio nel quadrilatero con la manovella considerata come movente:la coordinata indipendente definisce la posizione della manovella, mentre la posizionedi biella e bilanciere deve essere calcolata risolvendo un problema di posizione).

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

4 2 CONCETTI DI BASE

D’altro canto in alcuni casi l’efficacia computazionale e estremamente alta. Sonoutilizzate in applicazioni specifiche che riguardano principalmente sistemi a catenaaperta.La seconda possibilita e quella di usare un sistema “esteso” di coordinate dipendenti,collegate tra loro da equazioni in genere non lineari, note come equazioni di vincolo,il cui numero e pari alla differenza tra il numero di coordinate utilizzato e i gradi diliberta del sistema.Questo sistema di coordinate dipendenti puo essere costituito da coordinate relative(descrivono la posizione relativa di due corpi), reference point coordinate (descrivono laposizione di un singolo corpo rigido mediante la posizione di un suo punto e l’orienta-mento) o coordinate naturali (per descrivere la posizione di un corpo vengono utilizzatele coordinate di punti che appartengono al corpo stesso e/o direzioni di vettori ad essosolidali: la posizione di un corpo rigido nel piano e definibile ad esempio con due punti).

2.3 Formulazione simbolica–numerica

I programmi per l’analisi cinematica e dinamica si dividono in due grandi sottogruppi:i programmi “simbolici” e quelli “numerici”.I primi producono delle linee di codice (Fortran, Pascal, Matlab etc.) che contengo-no, in forma simbolica, le equazioni che descrivono la cinematica o la dinamica delsistema. La formulazione simbolica risulta vantaggiosa quando le equazioni di motosono valide per l’intero range di movimenti possibile. Uno dei problemi maggiori perquesti programmi deriva dal fatto che durante il movimento le equazioni di moto pos-sono cambiare radicalmente, non solo nei valori dei termini che le compongono, maanche qualitativamente, per la possibilita che alcune equazioni di vincolo scompaianoo appaiano (per la presenza di vincoli unilateri, giochi etc), o che si verifichino impattietc.I programmi numerici, che cioe formulano numericamente le equazioni di moto, con-sentono una gestione piu agevole di questi eventi ma non consentono di sviluppareapplicazioni “stand–alone”.

2.4 Classi di problemi

I problemi che possono essere affrontati e risolti tramite un programma multibody sonoi piu svariati; sono pero riconducibili ad alcune classi principali di cui sono riassunte leprincipali caratteristiche e peculiarita.

2.4.1 Analisi cinematica

L’analisi cinematica consiste nello studio della posizione o del movimento di sistemimulticorpo, a prescindere dalle forze agenti. E un problema puramente geometrico enon richiede percio le caratteristiche inerziali degli elementi che compongono il sistema(massa, momenti d’inerzia, posizione del baricentro).Si definiscono moventi quegli elementi la cui posizione o il cui movimento e noto (asse-gnato). Il numero dei moventi, nel caso di un problema cinematico, deve essere pari alnumero dei gradi di liberta del sistema multicorpo.I principali problemi che si affrontano sono:

Posizione iniziale o assemblaggio: consiste nel trovare la posizione di tutti i corpidel sistema nota quella dei moventi; comporta la soluzione di un sistema di equazioniin generale non–lineare e caratterizzato da piu soluzioni.

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

2 CONCETTI DI BASE 5

Spostamento finito: e una variante del precedente. Consiste nel trovare gli spo-stamenti dei corpi del sistema noti gli spostamenti degli elementi di input. E menodifficoltosa del precedente in quanto il calcolo degli spostamenti e fatto a partire dauna configurazione nota. Se il sistema non passa attraverso configurazioni singolarinon si verifica la possibilita di soluzioni multiple in quanto la modalita di assemblaggioe gia nota.

Analisi di velocita e accelerazione: consiste nel trovare velocita ed accelerazionedei corpi che compongono il sistema e dei punti rilevanti, note quelle dei moventi. Ilproblema e lineare ed ha in genere un’unica soluzione.

Simulazione cinematica: produce una vista d’insieme dell’intero range di movi-mento del sistema, comprende tutti i problemi sopra elencati. Permette di evidenziarecollisioni, studiare la traiettoria di punti etc.

2.4.2 Analisi dinamica

L’analisi dinamica consiste nello studio del movimento di sistemi in relazione alle forzeagenti; e in genere di soluzione piu impegnativa rispetto al problema cinematico. De-vono essere introdotti i carichi esterni agenti, le caratteristiche inerziali dei corpi, deglielementi elastici, viscosi, eventuali attriti etc.I principali problemi affrontati nella dinamica sono:

Posizione di equilibrio statica: consiste nel determinare la posizione del sistemain cui tutte le forze agenti si fanno equilibrio. La soluzione generale di questo problemaporta alla risoluzione di equazioni non lineari che devono essere risolte iterativamente.Possono essere presenti piu soluzioni (stabili o instabili)

Piccole oscillazioni: consiste nella determinazione dei modi e delle frequenze propriedel sistema nell’intorno della condizione di equilibrio statica (o dinamica). Questoproblema viene risolto linearizzando le equazioni di moto. Le frequenze proprie dannoun’idea della rigidezza del sistema; le equazioni linearizzate consentono di progettareun eventuale sistema di controllo.

Dinamica inversa: consiste nel calcolare le forze (coppie) richieste per effettuare undeterminato movimento. E’ necessario dapprima determinare le velocita e le accelera-zioni del sistema per calcolare le azioni di inerzia da cui dipendono, oltre che dai pesi,dalle forze negli elementi elastici etc. le forze (coppie) richieste agli attuatori.

Dinamica diretta (simulazione dinamica): fornisce il movimento di un sistemamulticorpo, in un dato intervallo di tempo, assegnate le forze applicate e le condizioniiniziali. Comporta l’integrazione di un sistema di equazioni differenziali nonlineari,a volte algebrico-differenziali, a partire da condizioni iniziali note. E’ un processo ingenere numericamente molto intensivo.Possono venire considerati problemi di impatto, cioe situazione in cui i corpi subisconovariazioni di velocita finite in un tempo infinitesimo.

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

6 3 COORDINATE ED EQUAZIONI DI VINCOLO

2.4.3 Sintesi

Nei casi precedenti si assume che le dimensioni del sistema siano note. Volendo proget-tare un sistema multicorpo, che risponda a certi requisiti, con soli strumenti di analisia disposizione, si deve procedere in modo iterativo, finche non si riesce a trovare unrisultato soddisfacente. L’esperienza e i risultati delle analisi ripetute possono guidareil progettista.Una analisi di sensitivita puo essere in questo caso di notevole aiuto per individuare iparametri che maggiormente hanno peso nella soluzione cercata.Una soluzione interamente automatizzata al progetto puo essere teoricamente ottenutamediante metodi numerici di ottimizzazione vincolata. A questo scopo deve esseredefinita una funzione (funzione obbiettivo) delle caratteristiche del sistema (variabili diprogetto) che il progetto finale deve rendere minima. La variazioni di queste variabilio di parametri ad esse legate in genere deve soddisfare delle condizioni che si possonotradurre con equazioni o disequazioni di “vincolo” (ad es. la lunghezza di un’asta deveessere compresa tra un minimo e un massimo, la distanza tra due punti deve rimanerecostante, etc.).Si tratta di un processo computazionalmente intensivo, tanto piu gravoso, quanto piuaumentano le variabili di progetto; non e poi garantita la convergenza del metodo versoun minimo assoluto.

3 Coordinate ed equazioni di vincolo

Esaminiamo ora in dettaglio i diversi tipi di coordinate dipendenti utilizzate – relative,reference point coordinate e naturali – con l’aiuto alcuni esempi che si riferiscono asistemi piani e 3D.

3.1 Sistemi Multibody piani

3.1.1 Coordinate relative

Definiscono la posizione di ciascun elemento relativamente a quello precedente nellacatena cinematica, utilizzando le coordinate corrispondenti ai gradi di liberta lasciatidal giunto che connette i due elementi (Paul (1970)). Nel caso di un sistema piano, sedue corpi sono uniti da una cerniera la posizione relativa e definita tramite un angolo,se da un giunto prismatico da una distanza (figura 3.1 e figura 3.2).Costituiscono un set con numero minimo di coordinate dipendenti; nel caso di unacatena aperta il numero di coordinate relative corrisponde ai gradi di liberta, non sonopresenti percio equazioni di vincolo.Nel caso di sistemi piani le equazioni di vincolo si derivano dalla condizione di chiusuravettoriale dei “loop” cinematici (indipendenti) contenuti e danno luogo a equazioni ingenere non lineari “piene”.Vettorialmente la condizione di chiusura per il quadrilatero di figura 3.1(b) risulta:

−→OA +

−−→AB +

−−→BD −−−→OD =

−→0

Questa equazione vettoriale e equivalente a due equazioni scalari che corrispondono allecomponenti lungo x e y:

{L1 cos Ψ1 + L2 cos(Ψ1 + Ψ2) + L3 cos(Ψ1 + Ψ2 + Ψ3)−OD = 0L1 sin Ψ1 + L2 sin(Ψ1 + Ψ2) + L3 sin(Ψ1 + Ψ2 + Ψ3) = 0

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

3 COORDINATE ED EQUAZIONI DI VINCOLO 7

Figura 3.1: Rappresentazione di un sistema articolato attraverso coordinate relative – (a)catena aperta, (b) catena chiusa

Figura 3.2: posizione relativa definita da una distanza

mentre per il meccanismo di figura 3.2 si ha:

−→OA +

−−→AB +

−−→BD −−−→OD =

−→0

Equivalente a:{

L1 cos Ψ1 + Ψ3 cos(Ψ1 + Ψ2) + L3 cos(Ψ1 + Ψ2 − π/2)−OD = 0L1 sin Ψ1 + Ψ3 sin(Ψ1 + Ψ2) + L3 sin(Ψ1 + Ψ2 − π/2) = 0

Gli aspetti positivi di questo approccio sono:

• un ridotto numero di coordinate

• una buona efficienza nei sistemi a catena aperta

• facilita nell’imporre moti relativi in corrispondenza dei giunti

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

8 3 COORDINATE ED EQUAZIONI DI VINCOLO

Svantaggi:

• la posizione di un elemento dipende da tutti quelli precedenti

• le equazioni di moto contengono matrici che, sebbene siano di ridotte dimensionisono piene e dispendiose da calcolare

• e richiesto un lavoro di pre-processing per individuare un set di equazioni divincolo indipendenti e di post-processing per determinare il moto assoluto diciascun elemento.

3.1.2 Reference point coordinate

Definiscono direttamente la posizione assoluta di un corpo con tre parametri (nel casopiano): le coordinate di un punto, spesso il baricentro, e l’orientamento del corporispetto ad un opportuno sistema di riferimento (figura 3.3).

Figura 3.3: Rappresentazione di un quadrilatero articolato con reference point coordinate –(a) cerniera, (b) giunto prismatico

In questa figura si nota come siano usate le stesse coordinate sebbene i meccanismisiano tra loro diversi.Le condizioni di vincolo vengono ottenute considerando i vincoli che i giunti introduconoal moto relativo tra elementi contigui. Una cerniera o un giunto prismatico da luogo adue equazioni di vincolo; un accoppiamento con due gdl relativi ad una sola equazione.Nel caso di un quadrilatero articolato piano (figura 3.3(a)), il sistema ha nove coordinatedipendenti e un solo grado di liberta, percio devono essere presenti otto equazionidi vincolo, due per ciascuna delle cerniere presenti. Immaginando che il punto diriferimento sia posto nella mezzeria delle aste le equazioni sono:

(x1 − x0)− L1/2 cos Ψ1 = 0(y1 − y0)− L1/2 sin Ψ1 = 0

(x2 − x1)− L1/2 cos Ψ1 − L2/2 cos Ψ2 = 0(y2 − y1)− L1/2 sin Ψ1 − L2/2 sin Ψ2 = 0(x3 − x2)− L2/2 cos Ψ2 + L3/2 cos Ψ3 = 0(y3 − y2)− L2/2 sin Ψ3 + L3/2 sin Ψ3 = 0

(x3 − xD)− L3/2 cos Ψ3 = 0(y3 − yD)− L3/2 sin Ψ3 = 0

In questo esempio si vede come le equazioni di vincolo siano sparse in quanto coinvolgonosolo le coordinate di due elementi adiacenti. Questa considerazione e comunque di

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

3 COORDINATE ED EQUAZIONI DI VINCOLO 9

carattere generale: le equazioni di vincolo di un giunto contengono solo le coordinatedi quei corpi che convergono nel giunto, sono definite cioe a livello locale; per ogni tipodi giunto le equazioni hanno sempre la stessa forma che non dipende dalla complessitadel sistema Si ottengono comunque equazioni non–lineari.In genere queste coordinate richiedono un numero piu elevato di variabili (9 inveceche 3, nel quadrilatero articolato). Catene chiuse e aperte vengono trattate allo stessomodo.

Vantaggi:

• la posizione di ciascun corpo e determinata direttamente

• le matrici che compaiono nelle equazioni di moto sono sparse, cioe hanno pochielementi non nulli e cio puo aumentare di molto l’efficienza numerica.

Svantaggi:

• numero elevato di coordinate, soprattutto nel caso di catene aperte

• “difficolta” nell’imporre un moto relativo ai giunti

3.1.3 Coordinate naturali

Possono essere viste come evoluzione delle reference point coordinate in cui i punti diriferimento vengono posizionati nei giunti o in altri punti importanti del corpo, in modoche ogni elemento sia caratterizzato da almeno due punti (Garcia (1981)). In questomodo la posizione e definita anche senza utilizzare variabili angolari, e le equazioni divincolo risultano spesso semplificate in quanto punti comuni hanno le stesse coordinate.Vediamo come puo’ essere definito un quadrilatero articolato tramite le coordinatenaturali. L’esempio si riferisce ad un quadrilatero articolato con un giunto prismatico(figura 3.4).

Figura 3.4: Quadrilatero articolato con giunto prismatico definito con coordinate naturali

Le coordinate sono 6: le coordinate dei punti 1,2 e 3; i gradi di liberta uno. Servonocinque equazioni di vincolo. Tre si ottengono imponendo che sia costante la lunghezzadei corpi 2 3 e 4:

(x1 − xA)2 + (y1 − yA)2 − L22 = 0

(x2 − x1)2 + (y2 − y1)2 − L23 = 0

(x3 − xB)2 + (y3 − yB)2 − L24 = 0

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

10 3 COORDINATE ED EQUAZIONI DI VINCOLO

Le altre equazioni vengono dal giunto prismatico: il punto 3 deve rimanere allineato aipunti 1 e 2 (il prodotto vettoriale tra

−→12 e

−→3B deve essere nullo):

(x3 − x1)(y2 − y1)− (x2 − x1)(y3 − y1) = 0

e inoltre l’angolo tra−→12 e

−→3B deve essere costante, cioe costante deve essere il prodotto

scalare tra questi vettori:

(x2 − x1)(x3 − xB) + (y2 − y1)(y3 − yB)− L3L4 cos Φ = 0

dove Φ e l’angolo formato dai due elementi.

3.2 Sistemi multicorpo spaziali

3.2.1 Coordinate relative

La maggiore differenza tra sistemi piani e 3D, riguardo all’uso di coordinate relative, stanella maggiore varieta di tipologie di giunti, che spesso permettono piu di un grado diliberta relativo. Per ognuno di essi e necessario introdurre una coordinata. Una cernierasferica introduce tre rotazioni, un giunto cilindrico una rotazione ed una traslazionelungo assi coincidenti etc.Alcuni autori “riducono” le tipologie possibili di vincoli, considerando alcuni di essicome combinazione di piu vincoli elementari a 1 gdl. Questo sostituzione e effettuataintroducendo elementi fittizi senza massa e di lunghezza nulla. Per esempio un giuntocilindrico puo essere sostituito da un giunto prismatico e da una cerniera introducendoun elemento intermedio.La robotica costituisce un campo di applicazione ove le coordinate relative sono moltoefficaci, in quanto i robot seriali sono costituiti da catene aperte con coppie rotoidali e/ogiunti prismatici, con un attuatore per ogni grado di liberta; in questo caso non sonopresenti equazioni di vincolo e le coordinate relative costituiscono un set di coordinateindipendenti.Per sistemi con catene chiuse le coordinate relative non sono indipendenti e le equazionidi vincolo vengono determinate dalla chiusura di anelli cinematici, formulata in generein termini matriciali invece che vettoriali come nel caso piano (ad esempio secondo lanotazione di Denavit e Hartenberg (1963)).

Figura 3.5: Rappresentazione di un elemento binario secondo Hartenberg e Denavit

Le equazioni di vincolo si ottengono applicando le matrici di trasformazione ad unanello chiuso (figura 3.6) e imponendo che il prodotto di queste matrici sia la matriceunita, come indicato in (3.1) (oppure interrompendo l’anello in corrispondenza di un

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

3 COORDINATE ED EQUAZIONI DI VINCOLO 11

elemento ed eguagliando la matrice di posizione dell’elemento rispetto alla base fissaottenuta con due diversi percorsi).

0T1(ψ1) 1T2(ψ2) 2T3(ψ3) 3T4(ψ4) 4T0(ψ0) = I (3.1)

Figura 3.6: Anello delle trasformazioni di Hartenberg e Denavit

Si genera in questo modo un set di equazioni sovrabbondanti, che si risolve con op-portuni metodi numerici (ad esempio nel senso dei minimi quadrati). In alternativa vaindividuato e risolto un sistema di equazioni tra loro indipendenti.

3.2.2 Reference point coordinate

Nel caso di sistemi 3D, per definire la posizione di un corpo si usano le coordinate carte-siane di un punto e l’orientamento di una terna solidale con esso, rispetto ad una ternafissa. Come noto il problema di definire l’orientamento relativo di due terne presentadiverse soluzioni (Argyris (1982)): per esempio con i nove elementi della matrice di ro-tazione, che rappresentano i coseni direttori di una terna rispetto all’altra (9 equazionilegate da 6 equazioni di vincolo, nessun subset di tre parametri sufficiente a definirela posizione in modo univoco), gli angoli di Eulero che portano il sistema mobile dal-la posizione iniziale a quella finale, mediante tre rotazioni successive attorno agli assiZ,x′,z′′ (figura 3.7), gli angoli Tait-Bryant (Cardano) di rollio, beccheggio e serpeggio(roll, pitch e yaw). Comunque tutti i sistemi di coordinate con tre parametri presentanouna singolarita, che viene evitata o con un posizionamento oculato della terna oppurepassando da una convenzione angolare ad un’altra in prossimita delle configurazionisingolari.Per questo motivo sono anche utilizzati sistemi basati su quattro parametri non indi-pendenti, che richiedono quindi l’introduzione di una equazione di vincolo.Il piu intuitivo fa uso dell’asse di rotazione finita: vengono utilizzati i versori dell’assedi rotazione finita e il valore dell’angolo di rotazione. L’equazione di vincolo tra iparametri impone che sia unitario il modulo del versore dell’asse di rotazione:

u2x + u2

y + u2z = 1

In pratica pero vengono utilizzate altre coordinate:

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

12 3 COORDINATE ED EQUAZIONI DI VINCOLO

Figura 3.7: Angoli di Eulero

• i parametri di Eulero

p1 = ux cos(Ψ/2) p2 = uy cos(Ψ/2) p3 = uz cos(Ψ/2) p4 = sin(Ψ/2)

legati dall’equazione di vincolo:

p21 + p2

2 + p23 + p2

4 = 1

• i parametri di Rodriguez -Hamilton;

• i quaternioni.

Alcuni programmi usano “coordinate” aggiuntive per le velocita (dette quasi-velocitao quasi-coordinate), soprattutto per quanto riguarda la velocita angolare, diverse dallederivate delle coordinate usate per descrivere la posizione. Si evitano in questo modoproblemi di singolarita sulle velocita e le accelerazioni. L’uso comunque di coordinateangolari e necessario in quanto la velocita angolare non e integrabile (o per lo meno ilsuo integrale non corrisponde all’orientamento del corpo).

Di seguito, come esempio, si riportano le equazioni di vincolo per un giunto sferico eper un giunto rotoidale.Due elementi uniti da un giunto sferico hanno sostanzialmente un punto in comune.Se P e questo punto, si e sj rappresentano la posizione del punto rispetto ai riferimentisolidali con i due corpi (il cui orientamento e espresso dalle matrici Ai e Aj), e ri e rj

rappresentano la posizione delle origini Oi e Oj dei due riferimenti, si ha:

ri + si − rj − sj = 0

dove tutti i vettori sono espressi nello stesso riferimento (in genere assoluto).Esprimendo i vettori s nel riferimento “locale” si ha:

ri + Aiisi − rj −Aj

jsj = 0

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

3 COORDINATE ED EQUAZIONI DI VINCOLO 13

che e equivalente a tre equazioni scalari, nelle coordinate dei punti O e nelle coordinateche esprimono l’orientamento dei corpi.Un giunto rotoidale impone che un punto sull’asse di rotazione relativo mantenga lastessa posizione rispetto ai due corpi connessi (come per il giunto sferico) e che i versoriiu e ju dell’asse di rotazione, riferiti ai singoli corpi, abbiano lo stesso orientamentorispetto ad un riferimento comune:

Aiiu−Aj

ju = 0

3.2.3 Coordinate naturali

Nel caso 3D le coordinate naturali descrivono la posizione di ciascun elemento permezzo delle coordinate cartesiane di punti “base” distribuiti per tutto l’elemento e conle componenti cartesiane di vettori unitari che definiscono l’orientamento del corpo.

3.2.4 Vincoli anolonomi

Negli esempi precedenti le equazioni di vincolo hanno tutte la forma generale

Φ(q, t) = 0

avendo indicato con q il vettore delle coordinate dipendenti.Nel caso dei vincoli anolonomi, le equazioni di vincolo legano tra loro semplicemente levelocita ed assumono la forma:

Φ(q,q, t) = 0

(in genere B(q)q+b = 0, con B matrice rettangolare) senza che esista una funzione Φtale che:

Bik =∂Φi

∂qkossia B = Φq

altrimenti significherebbe che le equazione di vincolo che legano le velocita sono ottenutesemplicemente dalla derivazione di Φ(q, t) = 0.Un esempio e dato da una ruota che rotola senza strisciare su un piano, mantenendosiad esso perpendicolare (figura 3.8).

Figura 3.8: Esempio di vincolo anolonomo

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

14 4 ANALISI CINEMATICA

La ruota ha 4 gdl, la posizione del centro, la rotazione attorno ad un asse perpendicolarealla ruota e la rotazione attorno ad un asse perpendicolare al piano. Tra queste 4coordinate esistono dei legami solo a livello di velocita:

{x = −R ϑ sin ϕ

y = R ϑ cosϕ

4 Analisi cinematica

Nella parte precedente si e visto, per alcuni casi, come vengono scritte le equazioni divincolo per un sistema multicorpo. Queste stesse equazioni di vincolo vengono poi uti-lizzate per risolvere problemi di carattere cinematico: posizione iniziale o assemblaggio,spostamento finito, analisi di velocita e accelerazione.I primi due problemi richiedono una soluzione iterativa, essendo le equazioni risolventinon lineari nelle coordinate. Particolari cautele dovranno essere utilizzate nei confrontidei sistemi sovra-vincolati, o in genere con equazioni di vincolo non indipendenti.

4.1 Assemblaggio

Consiste nel trovare la posizione di tutti i corpi nota la posizione degli elementi di input,tanti quanti sono i gradi di liberta del sistema.Il problema si riduce alla soluzione del sistema

Φ(q, t) = 0 (4.1)

dove q rappresenta il vettore delle coordinate del sistema.Immaginiamo che ci siano almeno tante equazioni quante sono le incognite, cioe lecoordinate.Per risolvere il sistema di equazioni non-lineari ci si basa sul metodo di Newton-Raphsonche ha una convergenza quadratica nell’intorno della soluzione e quindi, se si parte dauna buona approssimazione iniziale, sono sufficienti poche iterazioni per determinarela soluzione.Il metodo e basato sulla linearizzazione dell’equazione (4.1) attorno a una configurazioneapprossimata qi:

Φ(q, t) w Φ(qi,t) + Φq(qi,t)(q− qi) = 0 (4.2)

dove la matrice Φq e la matrice Jacobiana dell’equazione di vincolo:

Φq =

∂Φ1

∂q1

∂Φ1

∂q2........

∂Φ1

∂qn

∂Φ2

∂q1

∂Φ2

∂q2........

∂Φ2

∂qn

........ ........ ........ ........

∂Φm

∂q1

∂Φm

∂q2........

∂Φm

∂qn

(4.3)

in (4.3), m e il numero delle equazioni di vincolo e n il numero di coordinate. Se leequazioni di vincolo sono indipendenti, f = n −m corrisponde al numero di gradi diliberta del sistema.

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

4 ANALISI CINEMATICA 15

La (4.2) e un sistema di equazioni lineari e q, ottenuto risolvendo tale sistema, sara unaapprossimazione della soluzione della (4.1). Indicando questa soluzione approssimatacome qi+1 si ottiene una formula ricorsiva:

Φ(qi,t) + Φq(qi,t)(qi+1−qi) = 0 (4.4)

Il processo iterativo di soluzione verra arrestato quando l’errore in (4.1) diventa tra-scurabile e l’incremento delle q tra due iterazioni risulta inferiore ad una certa soglia.La figura 4.1 e una rappresentazione geometrica del metodo nel caso di un’equazionein una sola incognita.

Figura 4.1: Metodo di Newton–Raphson

Problemi di convergenza si possono avere se non si parte sufficientemente vicini allasoluzione o se le variabili di input non corrispondono ad una soluzione fisicamenterealizzabile.

Esempio: quadrilatero in coordinate naturali.Il sistema di figura 4.2 e caratterizzato da quattro coordinate naturali (x1, y1, x2, y2).Le equazioni di vincolo corrispondenti sono date dalla costanza della lunghezza delleaste:

Figura 4.2: Quadrilatero articolato in coordinate naturali

(x1 − xA)2 + (y1 − yA)2 − L22 = 0

(x2 − x1)2 + (y2 − y1)2 − L32 = 0

(x2 − xB)2 + (y2 − yB)2 − L42 = 0

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

16 4 ANALISI CINEMATICA

In questo caso l’equazione (4.4) assume l’espressione:

2

(x1 − xA) (y1 − yA) 0 0(x1 − x2) (y1 − y2) (x2 − x1) (y2 − y1)

0 0 (x2 − xB) (y2 − yB)

x1

y1

x2

y2

i+1

x1

y1

x2

y2

i

=

= −

(x1 − xA)2 + (y1 − yA)2 − L22

(x2 − x1)2 + (y2 − y1)2 − L32

(x2 − xB)2 + (y2 − yB)2 − L42

i

Per poter essere in grado di risolvere il problema, in questo sistema di equazioni almenouna delle 4 coordinate deve essere nota.Ad esempio se e nota x1 si puo scrivere:

(x1)i+1 − (x1)i = 0

e quindi la prima colonna della matrice jacobiana risulta moltiplicata per zero e puoessere percio eliminata.In figura 4.3 sono rappresentate la configurazione iniziale approssimata, le prime itera-zioni e la configurazione finale del quadrilatero.

Approssimazione iniziale

Configurazione finale

Prime iterazioni

Figura 4.3: Configurazioni relative alle diverse iterazioni

In questo tipo di problemi e fondamentale per una buona efficienza del programma lavalutazione e la triangolarizzazione dello jacobiano. Per migliorare i tempi di calcoloa volte si usano algoritmi modificati, che riutilizzano lo stesso jacobiano per un certonumero di iterazioni.In alcuni casi non e possibile individuare a priori quali siano i “moventi” del sistema, oaddirittura l’operatore puo “ignorare” il numero di gradi di liberta del sistema che staanalizzando: in questi casi si puo assegnare un valore iniziale q0 a tutte le coordinate delsistema, valori che in generale non saranno compatibili. La configurazione assemblatarisultera dalla soluzione di un problema di minimizzazione:

f(q) =12

(q− q0)T W(q− q0) (4.5)

soggetta ai vincoli espressi dallaΦ(q, t) = 0 (4.6)

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

4 ANALISI CINEMATICA 17

dove la matrice W in (4.5) e una matrice diagonale di pesi:

W = diag(w1, w2, . . . , wn)

Pesi elevati faranno sı che la coordinata corrispondente si discosti poco dal valoredesiderato, mentre pesi bassi permettono piu ampi scostamenti.La (4.6) viene linearizzata e definendo d = q− q0 il problema diventa minimizzare la

f(d) = dT Wd

soggetta a:Φ(q0,t) + Φq(q0,t)d = 0

Il problema si puo risolvere con il metodo dei moltiplicatori di Lagrange definendo:

F(d,λλλ) = f(d) + λλλT (Φ(q0,t) + Φq(q0,t)d)

e le condizioni di ottimo sono:(

∂F

∂d

)T

= 0(

∂F

∂λλλ

)T

= 0

che portano al sistema lineare:[

W ΦTq (q0)

Φq(q0) 0

] [dλλλ

]=

[0

Φ(q0)

]

La risoluzione di questo sistema porta ad un nuovo valore per q:

q = q0 + d

che costituisce il punto di partenza per un nuovo passo di iterazione, ponendo q0 = q.

4.2 Velocita

Per risolvere il problema di velocita (data la velocita degli elementi di input, trovarequella degli altri elementi), occorre innanzitutto definire il valore delle q, cioe la confi-gurazione del sistema. Le equazioni di vincolo possono poi essere differenziate rispettoal tempo, ottenendo l’equazione:

Φq(qi,t)q = −Φt (4.7)

dove Φq e lo Jacobiano definito dalla (4.3), le q sono le “velocita” (derivate dellecoordinate) e Φt e la derivata parziale rispetto al tempo delle equazioni di vincolo. Seil vincolo e “scleronomico”, cioe non dipende dal tempo, tale derivata e nulla.Definite le q dei moventi la (4.7) consente di trovare le ”velocita” del sistema.A differenza del problema di posizione la (4.7) risulta una equazione lineare nelle q, lasoluzione non deve essere iterativa e, in assenza di singolarita, e unica. Ad esempioper il quadrilatero articolato di figura 3.1(b), si ottiene:

[ −L1s1 − L2s12 − L3s123 −L2s12 − L3s123 −L3s123

L1c1 + L2c12 + L3c123 L2c12 + L3c123 L3c123

]

Ψ1

Ψ2

Ψ3

=

[00

]

dove: s1 = sin Ψ1, s12 = sin(Ψ1 + Ψ2), s123 = sin(Ψ1 + Ψ2 + Ψ3)c1 = cos Ψ1, c12 = cos(Ψ1 + Ψ2), c123 = cos(Ψ1 + Ψ2 + Ψ3)

Se una delle tre velocita e nota, le corrispondenti colonne dello Jacobiano possonoessere portate a destra dell’uguale. Si ottiene un sistema lineare in due equazioni condue incognite.

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

18 4 ANALISI CINEMATICA

4.3 Accelerazione

Il problema di accelerazione si risolve differenziando ulteriormente l’equazione sopra,ottenendo:

Φqq = −Φqq− Φt (4.8)

Se la posizione e la velocita sono note, risolvendo il sistema (4.8) si ottengono le ac-celerazioni q. Anche qui lo jacobiano Φq e lo stesso della (4.2) e (4.7), percio se egia stato “invertito” per la soluzione del problema di velocita, e sufficiente costruire iltermine noto. A differenza del problema di velocita, che e omogeneo se non ci sonovincoli reonomici (dipendenti dal tempo), il problema di accelerazione e sempre nonomogeneo, almeno finche le velocita non sono tutte nulle. Ulteriori derivazioni portanoalla determinazione del jerk e delle derivate successive delle q.

4.4 Spostamenti finiti e simulazione cinematica

Consiste nel trovare la configurazione del meccanismo, noti gli spostamenti (in generefiniti) dei moventi.E strettamente correlato al problema dell’assemblaggio ed e controllato dalle stesseequazioni. Per migliorare la stima della condizione iniziale da cui comincia il proce-dimento iterativo di soluzione delle equazioni di vincolo non lineari, vengono a volteeffettuate preliminarmente delle analisi di velocita ed accelerazione. Rispetto all’assem-blaggio il processo iterativo parte, di solito, da una buona stima della posizione inizialeche viene ottenuta a partire da una posizione che soddisfa esattamente le equazioni divincolo: sono cosı quasi esclusi i problemi derivanti soluzioni multiple. Ovviamente ilmetodo tendera a divergere se la posizione scelta non risulta raggiungibile (ad esempioper il posizionamento dell’end-effector di un robot al di fuori dello spazio di lavoro), ameno che non siano presi gli opportuni accorgimenti, come ad esempio, una soluzionenel senso dei minimi quadrati.

4.5 Sistemi con vincoli sovrabbondanti

In certe situazioni, soprattutto per motivi strutturali, i sistemi multicorpo vengonovincolati in modo sovrabbondante. A volte poi un meccanismo che risulta vincolatocorrettamente se considerato come piano risulta invece “iperstatico” se considerato in3D: ad esempio un quadrilatero articolato (costituito da 3 aste mobili tra loro collegatecon cerniere) possiede 1 gdl “nel piano” (3 gdl×3 aste−2 gdv×4 cerniere=1 gdl), men-tre nello spazio risulta doppiamente iperstatico (6 gdl×3 aste−5 gdv×4 cerniere= −2gdl). In effetti il quadrilatero si impunta a meno che gli assi delle cerniere non sianoconvergenti in un unico punto (al finito o all’infinito); nei sistemi reali la presenza dipiccoli giochi e della deformabilita fanno sı che questi sistemi comunque riescano afunzionare correttamente.Se un sistema ha n coordinate, e possiede f gradi di liberta effettivi, e vincolato dam = n− f equazioni indipendenti. Se invece m > n− f , il sistema presenta un numerodi vincoli sovrabbondante e, in genere, le equazioni di vincolo non ammettono soluzioni(l’assemblaggio del sistema risulta impossibile). Se pero le equazioni di vincolo sono“compatibili” la soluzione esiste e puo essere calcolata.Risolvendo il problema di posizione in modo iterativo, non e detto pero che il sistemalineare

Φ(qi,t) + Φq(qi,t)(qi+1−qi) = 0 (4.9)

comprenda equazioni tra loro compatibili (saranno compatibili solo in corrispondenzadella soluzione esatta). I metodi che si seguono sono sostanzialmente due: eliminazione

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

5 ANALISI DINAMICA 19

Figura 4.4: Traiettoria di un punto di biella di un quadrilatero articolato

delle equazioni linearmente dipendenti o risoluzione del sistema (4.9) nel senso deiminimi quadrati: (

ΦTqΦq

)i(qi+1−qi) = −(ΦT

q )i(Φ)i

5 Analisi dinamica

5.1 Equazioni di moto e dinamica diretta

La dinamica diretta consiste nella determinazione del movimento di un sistema multi-corpo risultante dall’applicazione di azioni esterne e/o dal movimento imposto da alcunigiunti azionati.A differenza della cinematica, dove ad ogni gdl era imposto un vincolo sul movimento,in dinamica il numero dei giunti controllati “cinematicamente” deve essere inferiore alnumero di gdl. Il movimento del sistema risulta dalle condizioni di equilibrio dinamicoche portano alla scrittura di un sistema di equazioni differenziali del secondo ordine,dette equazioni di moto.I procedimenti per generare le equazioni di moto sono molteplici: principio dei lavorivirtuali, principio di Hamilton, equazioni di Lagrange, potenze virtuali. In molti deisoftware commerciali viene utilizzata la formulazione di Lagrange che porta alla scrittu-ra del seguente sistema di equazioni, considerando un sistema descritto da n coordinate

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

20 5 ANALISI DINAMICA

dipendenti q e soggetto a m equazioni di vincolo:

d

dt

(∂L

∂q

)− ∂L

∂q+ ΦT

qλ = Qex (5.1)

Φ(q, t) = 0 (5.2)

che costituisce un sistema algebrico-differenziale (DAE).L e la Lagrangiana del sistema che e data dalla differenza tra l’energia cinetica T (q, q),che ha la forma generale:

T =12qT M(q)q

e potenziale V (q), dove Qex(q) e il vettore dei carichi generalizzati rispetto alle coordi-nate q. Il termine ΦT

qλλλ e introdotto dato che le coordinate non sono indipendenti e con-tiene le reazioni vincolari generalizzate. Il vettore λλλ (m×1) rappresenta i moltiplicatoridi Lagrange.Le equazioni (5.1) diventano cosı

Mq + ΦTqλ = Qex + Lq − Mq (con Lq = Tq −Vq) (5.3)

che rappresentano n equazioni in (n + m) incognite: gli n elementi del vettore q egli m elementi del vettore λλλ. Per avere un numero sufficiente di equazioni si possonoaccoppiare direttamente le m equazioni di vincolo, ottenendo un sistema DAE, oppurele loro derivate seconde, ottenendo l’equazione:

Φqq = −Φt − Φqq ≡ c

Operando in questo modo, e raggruppando alcuni termini si ottiene:[

M ΦTq

Φq 0

] [qλλλ

]=

[Qc

](5.4)

dove Q ≡ Qex + Lq − Mq(i termini Mq e Tq sono quadratici nelle velocita con coefficienti che possono dipenderedalle q. I termini contenenti q2

i sono detti centrifughi, quelli con (qiqj) di Coriolis. Iltermine Vq dipende da q ma non dalle sue derivate.)

Vediamo, a titolo di esempio, come puo essere calcolata l’energia cinetica di un corpo.Se r e la posizione di un generico punto e ρ la densita:

T =12

Vρ rT rdV

Per un corpo rigido la posizione r “assoluta”di un suo punto , la cui posizione rispettoad un riferimento solidale con esso sia u (figura 5.1), e data da:

r = R + A u

dove R rappresenta la posizione dell’origine della terna e la matrice A contiene i cosenidirettori della terna (espressi in funzione delle coordinate usate: angoli di Eulero, etc.)e la sovrasottolineatura denota matrici o vettori riferiti alla terna locale.La velocita del punto risulta essere:

r = R + A u

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

5 ANALISI DINAMICA 21

u

R

r

Figura 5.1: Riferimento assoluto e riferimento locale

ma A u, che rappresenta la velocita di trascinamento dovuta alla rotazione della terna,puo essere scritta come A(ω ∧ u), o anche −A ˜u ω, dove ˜u rappresenta la matriceemisimmetrica costruita col vettore u. In generale il legame tra la velocita angolare delcorpo e le derivate delle coordinate ”rotazionali” θ puo essere espresso come:

ω = G( θ) θ

Percio l’energia cinetica risulta essere:

T =12

[RT θT

] ∫

I −A ˜uG

−A ˜uG GT ˜uT ˜uG

dV

[Rθ

]

ossia:

T =12

qT

[mRR mRθ

mRθ mθθ

]q =

12

qT Mq

indicando con q l’insieme delle coordinate che definiscono la posizione del corpo.La matrice mRR risulta essere: diag(m,m, m), dove m e la massa del corpo, mRθ enulla se si usa un riferimento baricentrale (infatti questa matrice dipende dai momentistatici del corpo rispetto agli assi del riferimento) e mθθ = GT IθθG, dove Iθθ e il tensored’inerzia del corpo. Per un corpo e un movimento piano la matrice M ha dimensione(3×3) e il termine mθθ si riduce semplicemente al momento d’inerzia del corpo rispettoall’asse coordinato perpendicolare al piano del movimento.Effettuando le derivazioni previste dalle equazioni di Lagrange si ottengono i seguentitermini:

d

dt

(∂T

∂ q

)−

(∂T

∂ q

)= M q + M q−

(∂T

∂ q

)T

= M q + 2 ωT Iθθ˙G

Il sistema (5.4) puo essere risolto direttamente fornendo le accelerazioni e i moltiplicato-ri ad un certo istante e in una determinata configurazione q(t), q(t). Avendo a disposi-zione una routine di integrazione si puo determinare la configurazione q(t+∆t), q(t+∆t)del sistema al passo successivo.Altri metodi prevedono l’eliminazione dei moltiplicatori di Lagrange. Se con R(n, n−m)indichiamo una base per lo spazio nullo dello Jacobiano dei vincoli Φq, che consentedi definire tutti i movimenti possibili compatibili con i vincoli stessi, questa gode dellaproprieta che ΦqR = 0. Con la matrice R e possibile definire le velocita q del sistemain funzione di un set di velocita z tra loro indipendenti: q = Rz .

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

22 5 ANALISI DINAMICA

Premoltiplicando la (5.3) per RT si eliminano i moltiplicatori di Lagrange e si ottieneun sistema di (n−m) equazioni in n incognite:

RT M q + RT ΦTqλλλ = RT Q (5.5)

a cui vanno aggiunte le m equazioni di vincolo tra le accelerazioni. Si ottiene un sistemain cui le incognite sono solo le accelerazioni q.Entrambi questi metodi utilizzano pero le equazioni di vincolo differenziate due volte.Questo vuol dire che viene anche integrata l’equazione:

Φ(q, t) = Φqq + Φq q + Φt= 0 (5.6)

che ammette come soluzione generale

Φ(q,t) = a1t + a2 (5.7)

Se le condizioni iniziali soddisfano esattamente le equazioni di vincolo allora Φ(q,t) = 0e soddisfatta in ogni momento, in quanto a1 e a2 sono nulli.In realta la (5.7) e instabile perche per qualunque valore di a1 diverso da zero, lasoluzione cresce indefinitamente. Questo e quello che si verifica praticamente a causadegli inevitabili errori di arrotondamento.Per questo motivo non si integra mai semplicemente il sistema (5.5), con l’equazione divincolo (5.6), ma si ricorre all’integrazione del sistema DAE, alla chiusura ”periodica”dei vincoli, alla stabilizzazione dei vincoli (metodo di Baumgarte 1972), in cui la (5.6)e sostituita da

Φ+2αΦ + β2Φ = 0 (5.8)

con α e β scelti in modo che la (5.8) risulti un sistema smorzato criticamente (even-tuali errori di posizione si annullano nel tempo), o a metodi penalty (Bayo 1988) eLagrangiani “augmented”.In alternativa, utilizzando in modo appropriato le equazioni di vincolo, si possonoeliminare i moltiplicatori di Lagrange e scrivere le equazioni di moto in funzione di unset di coordinate indipendenti (Wehage e Haug (1982)). Si ottiene un sistema ODEche integrato da i valori delle coordinate indipendenti. Passo passo deve essere perorisolto un problema di posizione e velocita per determinare la posizione e la velocitadi tutti i corpi. Si tratta di un metodo particolarmente efficiente per sistemi in catenaaperta, per i quali , peraltro, le coordinate relative costituiscono direttamente un setdi coordinate indipendenti, dando luogo percio ad un sistema ODE.

5.2 Posizione di equilibrio statico

Consiste nel trovare la posizione di equilibrio statico di un sistema multibody che puocontenere elementi elastici ed e soggetto a un sistema di forze esterne come il peso, leforze centrifughe o in genere forze d’inerzia corrispondenti ad un campo noto di accele-razioni. (La posizione di equilibrio statico puo essere cercata ovviamente anche rispettoad un sistema di riferimento mobile). E’ un problema nonlineare. Nella posizione diequilibrio il sistema di forze deve essere equilibrato per ciascun elemento del sistema ele equazioni di vincolo devono pure essere soddisfatte.Per affrontare questo problema ci si basa sulla stazionarieta dell’energia potenzialetotale, a cui si aggiunge il contributo dei moltiplicatori di Lagrange (legato alle reazionivincolari), o su metodi tipo penalty (in cui le equazioni di vincolo sono introdottenell’energia del sistema come fattore di penalita) o potenze virtuali.

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

5 ANALISI DINAMICA 23

Una possibile traccia di soluzione impone che sia stazionaria l’energia potenziale Vdel sistema soggetto ai vincoli Φ(q) = 0. Questo problema, utilizzando il metodo deimoltiplicatori di Lagrange, coincide con la ricerca del minimo della funzione:

V ∗ = V + Φ(q)Tλλλ (5.9)

Questa e minima se e nulla la derivata di V ∗ rispetto a q e rispetto a λλλ:

∂V ∗

∂q=

∂V

∂q+ Φq(q)T λ = 0 ,

∂V ∗

∂λλλ= Φ(q) = 0 (5.10)

Si ottiene un sistema di equazioni nonlineari in q e λ:

Ψ(q, λ) = 0 (5.11)

Applicando il metodo di Newton-Raphson si ottiene un sistema di equazioni lineari:

[Ψq Ψλλλ

Φq 0

]

i

({qλλλ

}

i+1

−{

qλλλ

}

i

)= −

{ΨΦ

}

i

con Ψλλλ = ΦTq (5.12)

che consente di trovare le coordinate q e i moltiplicatori λλλ al passo i + 1. Si itera finoa convergenza.

5.3 Dinamica inversa

Consiste nel trovare le coppie e forze di azionamento e le reazioni ai giunti, dato ilmovimento del sistema e le forze esterne applicate. Diverse metodologie sono pos-sibili: metodo di Newton–Eulero (equilibrio di forze e coppie dei singoli membri odi sottoinsiemi di corpi opportunamente scelti), moltiplicatori di Lagrange e potenzevirtuali.Il metodo di Newton consiste essenzialmente nello scrivere le equazioni di equilibriodinamico del sistema, considerando le azioni motrici e le reazioni vincolari come in-cognite. Se non sono presenti vincoli sovrabbondanti c’e pareggio tra numero delleequazioni e delle incognite. E’ un metodo che puo essere applicato in modo moltoefficiente (ricorsivo) nel caso di catene aperte.Con il metodo dei moltiplicatori di Lagrange le equazioni di partenza sono:

Mq + ΦTqλ−Q = 0

Dove il primo termine contiene le azioni di inerzia dipendenti dalle q, il secondo rap-presenta le reazioni vincolari e le ”coppie” di azionamento, il terzo le forze esterne ele forze d’inerzia dipendenti dalla velocita. Ciascuna colonna di ΦT

qλλλ rappresenta ilvettore di forze (generalizzato) associato ad un vincolo.Saranno presenti anche le m equazioni di vincolo (sicuramente verificate) ed essendoil movimento del sistema noto (supponendo che non sia vincolato in modo sovrabbon-dante), ulteriori n−m equazioni di vincolo che impongono il movimento al sistema deltipo Φ(q, t) = 0. Si ottiene percio un sistema lineare in 2n equazioni nelle 2n incogniteq e λλλ.

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

24 6 INTEGRAZIONE DELLE EQUAZIONI DI MOTO

5.4 Linearizzazione delle equazioni di moto: piccole oscillazioni attor-no alla posizione di equilibrio statico

Le equazioni di moto di un sistema multibody sono in genere altamente non-lineari.La soluzione di queste equazioni e necessaria per la simulazione dinamica del sistemasoggetto a grandi spostamenti e rotazioni. Molti sistemi lavorano pero principalmentein prossimita di una configurazione di ”regime”. In questi casi risulta molto convenientelinearizzare le equazioni di moto attorno alla configurazione di equilibrio, cosı da poterutilizzare tutti gli strumenti che sono disponibili nel caso di sistemi lineari: sovrappo-nibilita degli effetti, calcolo di autovalori e autovettori, approccio modale, sintesi delcontrollore etc...Nel caso generale le equazioni di moto sono del tipo

H(q, q, q, Q) = 0

La configurazione di equilibrio dinamica, indicata dal pedice (0) soddisfera le equazionidi moto:

H0 ≡ H(q0, q0, q0, Q0) = 0

Considerando piccole perturbazioni l’equazione e linearizzata utilizzando i primi duetermini dello sviluppo in serie di Taylor:

H(q0 + ∆q0, q0 + ∆q0, q0 + ∆q0, Q0 + ∆Q0) ' H0 + ∆H0 =

= Hq∆q + Hq∆q + Hq∆q + HQ∆Q = 0(5.13)

I metodi per calcolare i termini della (5.13) sono differenti, possono essere convenientimetodi numerici, in quanto la linearizzazione e effettuata una sola volta e quindi puoessere accettata una bassa efficienza.

6 Integrazione delle equazioni di moto

Si e visto come l’applicazione delle leggi della dinamica ad un sistema multicorpo portialla scrittura di un sistema di equazioni algebrico-differenziale (DAE). Questo puo essereridotto a un sistema del secondo ordine di equazioni differenziali (ODE) o differenziandole equazioni di vincolo o utilizzando un set di coordinate indipendenti.Mentre per i sistemi ODE esistono da tempo algoritmi risolutivi, di cui sono ben notele caratteristiche di accuratezza, convergenza e stabilita, non e cosı per i sistemi DAE.Per questo molti programmi commerciali si basano sull’integrazione di sistemi ODEche possono essere messi nella forma generale:

q = F(q, q, t)

In genere gli algoritmi di integrazione si riferiscono a sistemi del primo ordine: a questoscopo le n equazioni del secondo ordine possono essere ridotte a 2n del primo ordinecon la posizione s = q, da cui {

s = F(q, s, t)q = s

(6.1)

cioe indicando con yT ={qT , sT

}

y = f(t, y)

Immaginando che la funzione f sia differenziabile rispetto a t e y esiste una sola soluzio-ne al problema (6.1), date le opportune condizioni iniziali y0. Differenti sono gli approc-ci per la soluzione numerica delle equazioni differenziali, con differenti caratteristichedi accuratezza e stabilita.

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

6 INTEGRAZIONE DELLE EQUAZIONI DI MOTO 25

6.1 Concetti di base

Un primo approccio per la ricerca della soluzione e quello di sviluppare in serie di Taylorla soluzione:

y(t) = y0 + ∆t y′(t0) +∆t2

2!y′′(t0) + · · ·

dove le derivate totali di y rispetto al tempo sono ottenute come:

y′ = f(t, y) y′′ = ft + fyf

La difficolta di valutare le derivate successive di f rende questo metodo troppo oneroso,a meno che lo sviluppo non venga troncato dopo pochi termini. Una di queste appros-simazioni e il ben noto metodo di Eulero che, nota la soluzione al passo n, approssimala soluzione al passo n + 1 con:

yn+1= yn+∆t f(tn, yn)

Per avere una buona accuratezza puo essere necessario usare passi molto piccoli, ren-dendo sensibile il problema dell’arrotondamento.Oltre all’accuratezza un’altra importante caratteristica del metodo di integrazione e lastabilita, cioe la capacita di mantenere limitati gli errori di integrazione. Un metodoinstabile fa crescere questi errori esponenzialmente e un overflow avviene dopo pochipassi.La stabilita dipende non solo dall’algoritmo, ma anche dal problema stesso. Per questo,parlando di stabilita, ci si riferisce all’equazione standard y′ = λy e la caratterizzazioneviene fatta definendo l’insieme di valori λ e ∆t per cui il metodo e stabile (figura 6.1).

Figura 6.1: Regioni di stabilita per metodi basati sulle formule di differenziazione all’indietro(BDF)

Gli algoritmi che sono stabili per un ristretto range di valori di (λ∆t) sono definiticondizionalmente stabili. Il passo di integrazione va scelto in base alle caratteristiche delproblema definite da λ. Ad esempio per il metodo di Eulero si deve avere ∆t <| λ | /2

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

26 6 INTEGRAZIONE DELLE EQUAZIONI DI MOTO

Algoritmi che sono stabili per Re(λ∆t) < −n, con n > 0 sono detti “stiffly stable” e sonoadatti per sistemi stiff, cioe sistemi con una dispersione di vari ordini di grandezza per λe in cui le equazioni con λ elevato non contribuiscono apprezzabilmente alla soluzione.L’integrazione con metodi condizionalmente stabili risulterebbe troppo onerosa.Se λ e immaginario, cioe, nel caso di un sistema multibody, se sono presenti motioscillatori non smorzati, i metodi “stiffly-stable” non sono adatti e si deve ricorrerea metodi A-stabili o a metodi condizionalmente stabili che includono parte dell’asseimmaginario.Un metodo e detto A-stabile se ha come regione di stabilita l’intero semipiano complessosinistro compreso l’asse immaginario. Questo vuol dire che non sono presenti limitazionisul passo per quanto riguarda la stabilita. Una sottoclasse e quella dei metodi cosiddettiL-stabili per i quali yn+1 = a(λ∆t)yn con a(λ∆t) → 0 se Re(λ) → ∞. Questi metodismorzano rapidamente le alte frequenze presenti in sistemi stiff. Comunque in genereintroducono smorzamento anche per frequenze piu basse.

6.2 Metodi tipo Runge-Kutta

Questi metodi ottengono lo stesso ordine di un metodo tipo Taylor sostituendo il calcolodelle derivate di f con la valutazione di f in vari punti del passo ∆t di integrazione. Ilpiu usato e l’algoritmo del quarto ordine:

yn+1 = yn +∆t

6(k1 + k2 + k3 + k4)

k1 = f(tn, yn)

k2 = f(tn +∆t

2, yn +

∆t

2k1)

k3 = f(tn +∆t

2, yn +

∆t

2k2)

k4 = f(tn + ∆t, yn + ∆tk3)

E un metodo cosiddetto esplicito perche tutti i ki dipendono da valori gia calcolati.Se cosı non e il metodo viene detto implicito e richiede la soluzione di un sistema diequazioni non lineare, come, ad esempio, il metodo dei trapezoidi (metodo del secondoordine):

yn+1 = yn +∆t

2[f(tn, yn) + f(tn+1, yn+1)]

L’errore di troncamento in questi metodi in genere viene valutato integrando lo stessointervallo con due passi diversi.

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

6 INTEGRAZIONE DELLE EQUAZIONI DI MOTO 27

6.3 Metodi multistep

I metodi di Taylor e Runge-Kutta sono detti a singolo passo, perche richiedono infor-mazioni solo sul passo n, n+1 per avanzare al passo successivo. Quando vengono usateinformazioni relative ai passi precedenti, il metodo e chiamato multistep.Questi metodi integrano l’equazione differenziale y = f(t, y) da tn−p a tn+1 dove p eun numero intero:

yn+1 = yn−p +∫ tn+1

tn−p

f(t, y)dt

Per valutare l’integrale si puo approssimare la funzione f(t, y) utilizzando i valori cal-colati ai passi precedenti. Si ottiene una famiglia di metodi che dipendono dll’ordinecon cui viene approssimata la f(t, y) e dal valore di p.Se l’approssimazione include il valore di f(t, y) al passo n + 1 il metodo e implicito.I metodi impliciti sono in genere piu accurati e maggiormente stabili rispetto a quelliespliciti, anche se sono piu difficile da usare in quanto per ricavare yn+1 e necessarioun processo iterativo.La formula generale di un metodo multistep e:

p+1∑

i=0

αiyn+1−i +k∑

i=0

∆t βi f(tn+1−i, yn+1−i) = 0 (6.2)

Se β0 = 0 il metodo e esplicito in quanto non compare il termine f(tn+1, yn+1).Un metodo esplicito molto utilizzato, condizionalmente stabile, e quello di Adams-Bashfort con p = 0. Il metodo del quarto ordine risulta essere:

yn+1 = yn +∆t

24(55fn − 59fn−1 + 37fn−2 − 9fn−3)

e l’errore di troncamento e esprimibile come E = (251/720)∆t5f IV (ξ).Richiede una sola valutazione di f per ogni passo di integrazione, al contrario dellequattro richieste dal metodo di Runge–Kutta dello stesso ordine. Il metodo, per poterpartire, richiede pero la valutazione di f in corrispondenza dei primi quattro passi, cosache puo essere fatta eventualmente ricorrendo a un metodo di RK.La famiglia di metodi impliciti con p = 0 sono detti di Adams-Moulton. Il metodo delquarto ordine e:

yn+1 = yn +∆t

24(9fn+1 + 19fn − 5fn−1 + fn−2) (6.3)

e l’errore di troncamento e esprimibile come E = −(19/720)∆t5f IV (ξ).L’uso della (6.3), o in genere di un metodo multistep implicito, richiede una soluzionedi tipo iterativo; in genere come punto di partenza si usa l’approssimazione data dalla(6.2), cioe da un metodo multistep esplicito, ottenendo un metodo di tipo predictor-corrector. Osservando l’errore di troncamento si puo dire che il metodo implicito e moltopiu accurato e quindi conviene generalmente ricorrere a questo metodo. Sfruttando ivalori di yn+1 ottenuti tramite il predictor e il corrector, e gli errori di troncamentodei due metodi, si ha una stima dell’errore dopo una iterazione, rendendo possibile unefficace controllo sul passo di integrazione.

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

28 6 INTEGRAZIONE DELLE EQUAZIONI DI MOTO

6.4 Metodo di Newmark

Nel caso di un sistema multibody con equazioni “stiff” puo essere conveniente ricorrerea metodi A-stabili, in quanto consentono di scegliere un passo adeguato per coglierelo svolgersi di un fenomeno (e quindi con buona accuratezza), senza dare problemi distabilita, come invece accade per i metodi solo condizionalmente stabili.Uno dei piu usati, soprattutto in campo strutturistico, e quello di Newmark che si basasulle seguenti interpolazioni che legano le posizioni q, le velocita q e le accelerazioni qdal passo n al passo n + 1:

qn+1 = qn + ∆t [(1− γ) qn + γqn+1] (6.4)

qn+1 = qn + ∆t qn +∆t2

2[(1− 2β) qn + 2β qn+1] (6.5)

i paramentri β e γ definiscono il metodo. Si tratta di un metodo implicito, A-stabilese 2β ≥ γ ≥ 1/2.La regola dei trapezoidi e un caso particolare con β = 1/4 e γ = 1/2, e implica chel’accelerazione rimanga costante nel passo e sia pari a: (qn+1 + qn)/2. Similmente aimetodi multistep, puo essere usato in modo predictor-corrector, oppure, introducendole interpolazioni direttamente nelle equazioni di moto, cosa che puo essere facilmenteillustrata applicando il metodo ad un problema strutturale lineare che assuma la formaseguente:

Mqn+1 + Cqn+1 + Kqn+1 = F(t)

dove M, C, K sono le matrici di massa, smorzamento e rigidezza e F il vettore deicarichi esterni. Sostituendo le (6.4) e (6.5) si ottiene:

[1

β∆t2 M+γ

β∆tC + K

]qn+1 =

= F(t) + M[

1β∆t2 qn +

1γ∆t

qn +

(1− 1

)qn

]+

+C[

γ

β∆tqn +

β− 1

)qn +

2β− 1

)∆tqn

](6.6)

Nel caso di un sistema multibody questa sostituzione porta ad un sistema nonlinearenelle qn+1 che va risolto in modo iterativo, o con il metodo di Newton–Raphson o conlo stile predictor–corrector.

6.5 Integrazione diretta delle DAE

Le equazioni di moto di un sistema multibody vincolato costituiscono un sistema DAEsemi-esplicito di indice 3 che puo essere posto nella forma:

F(t, y, y′) = 0

(Con indice di un sistema DAE si intende il numero di volte che il sistema deve esseredifferenziato per ottenere un sistema ODE).Due famiglie di integratori sono utilizzate: i BDF e i RK impliciti.In essenza si tratta di approssimare y′ tramite differenze finite o con metodi RK implicitie di risolvere il sistema di equazioni risultante con una procedura di tipo iterativo. Peresempio, usando la formula di Eulero implicito yn+1 = yn + ∆t yn+1, si ottiene:

F[tn+1, yn+1,

1∆t

(yn+1 − yn)

]= 0

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

7 SINTESI CINEMATICA E DINAMICA 29

In genere si usano approssimazioni del tipo

qn+1 =1

∆t β0

[qn+1 −

p∑

i=0

αi qn−i

]

sn+1 =1

∆t β0

[sn+1 −

p∑

i=0

αi sn−i

]

che vengono sostituite nelle equazioni di moto (5.3) [scritta come (6.1)] e (5.2) ottenendo

M [qn+1]1

∆t β0

[sn+1 −

p∑

i=0

αi sn−i

]= Qn+1 −ΦT

q [qn+1]λλλn+1

1∆t β0

[qn+1 −

p∑

i=0

αi qn−i

]= sn+1

Φ (t, qn+1) = 0

7 Sintesi cinematica e dinamica

Nei capitoli precedenti, in cui si e parlato dell’analisi cinematica e dinamica di un siste-ma multicorpo, si e supposto che le dimensioni e le caratteristiche del sistema fosseronote. Nel caso di progetto di un nuovo sistema, invece, queste ultime rappresentanole incognite del problema e devono essere determinate in funzione delle specifiche diprogetto.Nel caso in cui siano disponibili solo strumenti di analisi, si deve procedere per “tenta-tivi”, basandosi sui risultati di analisi ripetute. Il processo di sintesi puo essere moltolento e dipende fortemente dall’esperienza del progettista.Per questo motivo sono stati sviluppati metodi di sintesi o progetto consentono disuperare (parzialmente) questa difficolta in quanto portano direttamente alla soluzionedel problema.La sintesi cinematica di un meccanismo e un problema essenzialmente geometrico, sucui molto e stato scritto (ad esempio Angeles (1982), Erdman e Sandor (1978), Suhe Radcliffe (1978), Magnani Ruggieri (1986)), e puo essere ricondotta a tre classi diproblemi: generazione di funzioni, generazione di traiettorie e posizionamento di uncorpo rigido (figure 7.1, 7.2 e 7.3).

yj

Figura 7.1: Quadrilatero articolato generatore di funzione

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

30 7 SINTESI CINEMATICA E DINAMICA

Figura 7.2: Quadrilatero articolato generatore di traiettorie

Figura 7.3: Quadrilatero articolato generatore di moti piani

Sono disponibili sia metodi grafici, in genere molto specifici, che numerici di sintesiottima che presentano invece una buona generalita dal punto di vista teorico, anche sesono penalizzati in genere da lunghi tempi di calcolo e dalla difficolta di convergenza.Per applicare questi metodi “automatici” si deve definire la topologia del sistema,scegliere le variabili di progetto. identificare e minimizzare una funzione obbiettivo.La topologia del sistema e le specifiche di progetto determinano rispettivamente i vincolicosiddetti geometrici e funzionali:, in generale le equazioni di chiusura del meccanismoportano alla scrittura di equazioni del tipo:

Φ(q, b) = 0

dove q e b sono i vettori delle coordinate dipendenti e delle variabili di progetto.Le relazioni funzionali (se di eguaglianza), portano alla scrittura di equazioni del tipo

Φi(qi, b) = 0 i = 1, 2 · · ·ndove qi indica la posizione in cui il meccanismo deve soddisfare la i-esima condizione divincolo. In genere si cerca la soluzione a questo sistema nel senso dei minimi quadrati:

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

8 ELEMENTI FLESSIBILI 31

la funzione obbiettivo assume la forma:

Ψ(q1, q2, · · · , qN , b) =12

N∑

i=1

ΦiT (qi, b)Φi(qi, b)

o in forma compatta:

Ψ(q, b) =12ΦT (q, b)Φ(q, b)

Il problema di ottimo consiste nel trovare il minimo di Ψ rispetto a q e b.Il caso dinamico e ancora piu complesso. La funzione obbiettivo, in questo caso, edefinita spesso come l’integrale di una funzione o come una serie di condizioni da ri-spettare entro certi intervalli o in specifiche posizioni (o istanti) e quindi ogni “passo”di ottimizzazione richiede una simulazione completa del sistema.Anche qui si fa ricorso a metodi di ottimizzazione, che si basano (spesso) sulla cono-scenza delle derivate della funzione obbiettivo rispetto alle variabili di progetto. La de-terminazione di queste derivate e chiamata analisi di sensitivita. L’analisi di sensitivita,che determina la tendenza della funzione obbiettivo rispetto alle variabili di progettorisulta molto utile anche in un processo interattivo (non-automatico) di progetto.

8 Elementi flessibili

Quando le condizioni di funzionamento sono caratterizzate da prestazioni dinamichemolto elevate, il modello a corpi rigidi non e piu adeguato. Per studiare in manieraaccurata questi sistemi, si deve tener conto degli effetti della deformabilita dei varimembri e costrure quindi un modello ad elementi flessibili.Membri flessibili di forma complessa, nei programmi multibody, vengono introdotti ingenere interfacciando un programma di tipo FEM. Elementi snelli, tipo trave, possonoessere presenti nella “libreria” dei programma multibody.Sono seguiti principalmente due approcci di diversa complessita: un approccio cineto-statico ed uno di sintesi modale.Nel primo l’analisi dinamica del sistema multicorpo viene fatta considerando tutti imembri come rigidi: questo consente di ricavare le reazioni vincolare e le azioni d’i-nerzia. Questi carichi, insieme ai carichi esterni vengono “trasferiti” sul modello FEMche calcola sforzi, deformazioni e spostamenti dei corpi flessibili; il sistema di caricoviene applicato “staticamente” al corpo deformabile e quindi effetti dinamici come larisonanza non possono venire evidenziati con questo approccio.Il secondo approccio e piu complesso, ma consente di evidenziare anche questi fenomeni.Anch’esso poi consente di calcolare sforzi e deformazioni nei corpi flessibili.Brevemente il metodo “modale” (Craig-Bampton (1968)) prevede questi passi:la massa del corpo deformabile viene opportunamente concentrata nei nodi utilizzandoelementi con matrici di massa di tipo “lumped” (diagonali).Gli spostamenti dei nodi sono definiti mediante un set di coordinate generalizzate e dicorrispondenti “autovettori” che viene costruito dal programma FEM considerando imodi propri di vibrare del continuo vincolato in corrispondenza dei nodi di interfacciae le deformate statiche in corrispondenza di spostamenti unitari degli stessi nodi.Gli spostamenti dei nodi dovuti alla flessibilita sono percio esprimibili come:

u = N q

dove N e la matrice degli “autovettori” e q il vettore delle coordinate generalizzate.

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo

32 9 CONCLUSIONI

In genere i modi a frequenze piu alte, o comunque quelli non eccitati, portano uncontributo trascurabile alla deformazione del corpo: risulta percio possibile usare per lasimulazione dinamica un set ridotto di autovettori e quindi di coordinate generalizzate.La posizione r di un nodo, rispetto ad un riferimento fisso, risulta quindi (considerandoanche gli spostamenti rigidi del corpo):

ri= R + A(ui+Niqf )

dove ui e la posizione “indeformata” del nodo i e Ni e la parte della matrice N che siriferisce al nodo i.Vediamo, come esempio, il calcolo dell’energia cinetica T di un corpo flessibile. Lavelocita di un nodo risulta essere, utilizzando un procedimento analogo a quanto vistoper un punto generico di un corpo rigido:

ri = R−A(˜ui + Niqf )G(θ)θ + ANiqf

e quindi

T =12

[RT θT qT

f

] {N∑

i=1

mi

[I −A(˜ui + Niqf )G(θ) ANi

]T

[I −A(˜ui + Niqf )G(θ) ANi

] }[

R θ qf

]

Si ottiene:

T =12

qT

mRR mRθ mRf

mRθ mθθ mθf

mRf mθf mff

q =

12

qT Mq

dove gli elementi della matrice M dipendono dalla posizione q del corpo e da un set diinvarianti (9) che possono essere precalcolati all’inizio della simulazione.

9 Conclusioni

I programmi multibody e CAD, che contengono moduli per la simulazione cinematicae dinamica di sistemi meccanci, si vanno sempre piu diffondendo tra i progettisti; lacomplessita e il dettaglio dei sistemi che possono essere analizzati cresce di pari passo,sia per l’aumento della potenza di calcolo disponibile che per lo sviluppo di algoritmisempre piu efficienti.E cosı possibile affinare sempre piu i modelli tenendo conto anche dei giochi negli ac-coppiamenti, considerando la presenza di elementi flessibili di forma complessa e analiz-zando problemi di collisioni e contatto. Inoltre, per sistemi di particolare complessitae interesse come, per esempio, autoveicoli o parti di essi (sospensioni, trasmissione,pneumatici), veicoli ferroviari o aeromobili, i principali codici commerciali mettono adisposizione modelli “parametrici” che evitano al progettista di dover costruire “da ze-ro” ogni dettaglio del modello.Inoltre, sulla base della sempre piu diffusa integrazione tra meccanica, elettronica eazionamenti (meccatronica), i codici non si limitano piu ad analizzare sistemi pura-mente meccanici, ma offrono pacchetti per l’integrazione di azionamenti e di controllinel modello meccanico, oltre a permettere l’interfacciamento con software dedicati comeMatLab, Simulink, MatrixX etc.. Tutto cio fa parte del cosiddetto “virtual prototy-ping” che dovrebbe consentire di analizzare e predire il comportamento di un sistema

Metodi numerici per lo studio di sistemi multicorpo V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

10 BIBLIOGRAFIA 33

prima della sua realizzazione, evitando i costi per lo sviluppo e i test di modelli fi-sici. Bisogna pero ricordare che modelli complessi e dettagliati richiedono numerosiparametri di ingresso i cui valori numerici sono spesso di difficile valutazione. Inoltrel’interpretazione dei risultati puo risultare ambigua come puo esserlo per un sistemareale. Come in tutti i campi e percio necessario non perdere il contatto col mondo“reale”, sviluppando confidenza e sensibilita tramite modelli semplici, e mantenendosempre un atteggiamento critico nei confronti dei risultati.

10 Bibliografia

Argyris, J., “An Excursion into Large Rotations”, Computer Methods in Applied Me-chanics and Engineering, Vol. 32, pp. 85-155 (1982).

Garcia de Jalon, J., Unda, J., Avello, A., and Jimenez, J.M., “Dynamic Analysis ofThree-Dimensional Mechanisms in Natural Coordinates”, ASME J. of Mechanisms,Transmissions and Automation in Design, Vol. 109, pp. 460-465, (1987).

Hartenberg, R.S. and Denavit, I., “Kinematic Synthesis of Linkages”, McGraw-Hill,New York, (1964).

Paul, B. and Krajcinovic, K., “Computer Analysis of Machines with Planar Motion”,ASME Journal of Applied Mechanics, Vol. 37, pp. 697-712, (1970).

Nikravesh, P.E., “Computer-Aided Anslysis of Mechanical Systems”, Prentice-Hall(1988).

Shabana, A.A., Dynamics of Multibody Systems, Wiley (1989).

Schielen, W.O., Multibody System Handbook, Springer-Verlag (1990).

Wittenburg, J, Dynamics of Systems of Rigid Bodies, Teubner, Stuttgart, (1977).

R.R. Craig and M.C.C. Bampton, “Coupling of substructures for dynamics analyses,AIAA Journal, 6(7): 1313–1319 (1968).

G. Diana, F. Cheli, Cinematica e dinamica dei sistemi multi–corpo, ed. Spiegel (1998).

G.Legnani, F.Casolo, P.Righettini, B.Zappa, “A Homogeneous Matrix Approach to 3DKinematics and Dynamics - I . Theory”, Mechanisms and Machine Theory, Ed.Pergamon,Great Britain , Vol31, No.5, (1996).

G.Legnani, F.Casolo, P.Righettini, B.Zappa “A Homogeneous Matrix Approach to 3DKinematics and Dynamics - II . Applications to chains of rigid bodies and serial ma-nipulators”, Mechanisms and Machine Theory, Ed.Pergamon, Great Britain , Vol 31,No.5, (1996).

V. Lorenzi, R. Riva, M. Camposaragna, R. Strada, B. ZappaUniversita di Bergamo

Metodi numerici per lo studio di sistemi multicorpo