METODI DI DECOMPOSIZIONE DEI DOMINI APPLICATI ALLA ... · 3.8 In alto la curva statica di pull-in,...

POLITECNICO DI MILANO

Facoltà di Ingegneria Civile, Ambientale e Territoriale

Corso di Laurea Specialistica in Ingegneria Civile

Indirizzo Strutture

METODI DI DECOMPOSIZIONE DEI DOMINI

APPLICATI ALLA SOLUZIONE DI

PROBLEMI ELETTRO-MECCANICI

Relatori: Prof. Ing. Alberto CORIGLIANO

Dott. Ing. Federica CONFALONIERI

Tesi di Laurea di:

Martino DOSSI Matr. 735220

Matteo GORNATI Matr. 735370

Anno Accademico 2009/2010

Indice

1 Introduzione 1

2 La decomposizione dei domini: stato dell’arte 5

2.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5

2.2 Approccio per sottostrutture . . . . . . . . . . . . . . . . . . 7

2.3 Metodi di decomposizione dei domini . . . . . . . . . . . . . . 10

2.3.1 Formulazione del problema di interfaccia . . . . . . . . 11

2.3.2 Il metodo BDD . . . . . . . . . . . . . . . . . . . . . . 17

2.3.3 Il metodo FETI . . . . . . . . . . . . . . . . . . . . . . 19

2.4 L’algoritmo di Gravouil e Combescure . . . . . . . . . . . . . 27

2.4.1 L’algoritmo GC applicato a due sottodomini . . . . . 28

2.4.2 L’algoritmo GC applicato a s sottodomini con diversi

passi temporali . . . . . . . . . . . . . . . . . . . . . . 33

3 Problema elettro-meccanico accoppiato 35

3.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.2 Problema elettrostatico . . . . . . . . . . . . . . . . . . . . . 36

3.2.1 Equazioni di Maxwell ed equazione elettrostatica . . . 36

3.2.2 Formulazione ad elementi finiti del problema elettro-

statico . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

3.2.3 Esempio: problema elettrico 1D . . . . . . . . . . . . . 43

3.2.4 Le forze elettrostatiche . . . . . . . . . . . . . . . . . . 44

3.3 Problema elettro-meccanico . . . . . . . . . . . . . . . . . . . 47

3.3.1 Considerazioni iniziali . . . . . . . . . . . . . . . . . . 47

3.3.2 Formulazione del problema elettro-meccanico accoppiato 49

3.3.3 Dominio meccanico . . . . . . . . . . . . . . . . . . . . 51

3.3.4 Deformazione del dominio elettrico: problema pseudo-

meccanico . . . . . . . . . . . . . . . . . . . . . . . . . 51

I

INDICE II

3.3.5 Esempio: problema elettro-meccanico 1D . . . . . . . 55

3.4 Instabilita elettro-meccanica: fenomeno di

pull-in . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

3.4.1 Caso statico . . . . . . . . . . . . . . . . . . . . . . . . 57

3.4.2 Caso dinamico . . . . . . . . . . . . . . . . . . . . . . 60

4 Metodi di risoluzione 64

4.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

4.2 Metodo Monolitico . . . . . . . . . . . . . . . . . . . . . . . . 66

4.2.1 Esempio: metodo monolitico applicato ad un caso 1D 70

4.3 Metodi partizionati . . . . . . . . . . . . . . . . . . . . . . . . 74

4.3.1 Metodo staggered e metodo delle iterazioni simultanee 74

4.4 Voltaggio di pull-in . . . . . . . . . . . . . . . . . . . . . . . . 78

4.5 Esempio: problema 1D . . . . . . . . . . . . . . . . . . . . . . 79

5 La decomposizione dei domini applicata al problema elettro-

meccanico 88

5.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

5.2 Primo livello di decomposizione: contributi

free e link . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

5.2.1 Considerazioni iniziali . . . . . . . . . . . . . . . . . . 89

5.2.2 Formulazione del nuovo algoritmo . . . . . . . . . . . 91

5.3 Secondo livello di decomposizione: algoritmo GC applicato

alla risoluzione del problema meccanico . . . . . . . . . . . . 97

5.3.1 Decomposizione del problema meccanico free . . . . . 97

5.3.2 Decomposizione del problema meccanico link . . . . . 99

5.4 Terzo livello di decomposizione: metodo FETI applicato al

problema elettrico e pseudo-meccanico . . . . . . . . . . . . . 101

5.4.1 Metodo FETI applicato al problema pseudo-meccanico 103

5.4.2 Metodo FETI applicato al problema elettrico . . . . . 112

6 Applicazioni 120

6.1 Introduzione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

6.2 Casi prova . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

6.2.1 Trave a doppio incastro . . . . . . . . . . . . . . . . . 121

6.2.2 Trave a mensola . . . . . . . . . . . . . . . . . . . . . 145

6.2.3 Considerazioni finali . . . . . . . . . . . . . . . . . . . 156

6.3 Casi notevoli . . . . . . . . . . . . . . . . . . . . . . . . . . . 157

INDICE III

6.3.1 Esempio 1: trave a doppio incastro soggetta ad attua-

zione sinusoidale . . . . . . . . . . . . . . . . . . . . . 157

6.3.2 Esempio 2: portale semplice con attuazione laterale . 161

6.3.3 Esempio 3: trave a doppio incastro con smorzamento . 177

6.3.4 Esempio 4: trave a doppio incastro con doppia attua-

zione . . . . . . . . . . . . . . . . . . . . . . . . . . . . 180

7 Conclusioni 184

A Requisiti matematici degli algoritmi di decomposizione dei

domini 187

A.1 Metodo del gradiente coniugato . . . . . . . . . . . . . . . . . 187

A.2 Gradiente coniugato precondizionato (PGC) . . . . . . . . . . 189

B Integrazione delle equazioni del moto 192

B.1 Algoritmi impliciti . . . . . . . . . . . . . . . . . . . . . . . . 196

B.1.1 Metodo di Newmark dell’accelerazione costante . . . . 197

B.1.2 Metodo di Newmark dell’accelerazione lineare . . . . . 198

C Caratteristiche dei processori 199

D Descrizione dei codici 200

D.1 Decomposizione di 1 livello . . . . . . . . . . . . . . . . . . . 200



Elenco delle figure

1.1 Esempio di micro-sistema. . . . . . . . . . . . . . . . . . . . . 2

2.1 Problema di Schwarz ([10]). . . . . . . . . . . . . . . . . . . . 6

2.2 Esempio di suddivisione in sottostrutture. . . . . . . . . . . . 8

2.3 Problema di riferimento ([10]). . . . . . . . . . . . . . . . . . 12

2.4 Partizione in due sottodomini ([10]). . . . . . . . . . . . . . . 13

2.5 Interfaccia primaria e duale ([10]). . . . . . . . . . . . . . . . 14

2.6 Esempio di suddivisione in tre sottodomini ([10]). . . . . . . . 16

2.7 Definizione dei modelli di spostamento e delle forze di inter-

faccia. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.8 Esempio di partizione con sottodominio fluttuante. . . . . . . 23

2.9 Esempio di diverse discretizzazioni temporali ([2]). . . . . . . 32

3.1 Modello di accoppiamento elettro-meccanico ([22]). . . . . . . 36

3.2 Discretizzazione del dominio elettrico. . . . . . . . . . . . . . 40

3.3 Condensatore 1D con armature piane parallele ([22]). . . . . . 44

3.4 Problema elettro-meccanico di riferimento. . . . . . . . . . . . 49

3.5 Spostamenti meccanici, elettrici e di interfaccia. . . . . . . . . 55

3.6 Micro-risonatore monodimensionale ([22]). . . . . . . . . . . . 56

3.7 Andamento della curva spostamento-voltaggio dell’esempio

1D ([22]). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

3.8 In alto la curva statica di pull-in, in basso la risposta oscilla-

toria per un passo di voltaggio costante VD ([22]). . . . . . . . 60

3.9 Pull-in dinamico nel caso 1D: in alto la curva di pull-in statico,

in basso la risposta transitoria per V = VDPI ([22]). . . . . . 63

4.1 Approccio co-simulation (a sinistra) e fully coupled-simulation

(a destra). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

IV

ELENCO DELLE FIGURE V

4.2 Diagramma di flusso: metodo monolitico. . . . . . . . . . . . 69

4.3 Modello di riferimento monodimensionale. . . . . . . . . . . . 71

4.4 Procedura staggered ([22]). . . . . . . . . . . . . . . . . . . . 76

4.5 Diagramma di flusso: metodo staggered. . . . . . . . . . . . . 77

4.6 Ricerca del pull-in per potenziale imposto ([22]). . . . . . . . 78

4.7 Ricerca del pull-in con il metodo di Riks-Crisfield. . . . . . . 79

4.8 Esempio monodimensionale. . . . . . . . . . . . . . . . . . . . 81

4.9 Algoritmo di Newton-Raphson (istante n). . . . . . . . . . . . 82

4.10 Spostamento verticale nodo 4: soluzioni monolitica e staggered. 83

4.11 Potenziale nel nodo 2. . . . . . . . . . . . . . . . . . . . . . . 86

4.12 Pull-in dinamico nel caso monodimensionale. . . . . . . . . . 87

5.1 Decomposizione elettro-meccanica (a) e puramente meccanica

(b). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

5.2 Decomposizione di 1 livello: problemi free e link. . . . . . . . 95

5.3 Algoritmo di 1 livello: decomposizione secondo i due contri-

buti free e link. . . . . . . . . . . . . . . . . . . . . . . . . . . 96

5.4 Decomposizione della parte meccanica del problema di riferi-

mento. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

5.5 Decomposizione di 2 livello: schema di risoluzione del pro-

blema free meccanico. . . . . . . . . . . . . . . . . . . . . . . 100


blema link meccanico. . . . . . . . . . . . . . . . . . . . . . . 102

5.7 Decomposizione della parte elettrica del problema di riferimento.103

5.8 Teorema di composizione della velocita. . . . . . . . . . . . . 108


blema pseudo-meccanico. . . . . . . . . . . . . . . . . . . . . 111


blema elettrico free. . . . . . . . . . . . . . . . . . . . . . . . 116


blema elettrico link. . . . . . . . . . . . . . . . . . . . . . . . 117

5.12 Schema riassuntivo dei livelli di decomposizione introdotti. . . 119

6.1 Microbridge (a) e Comb Drive (b). . . . . . . . . . . . . . . . 121

6.2 Geometria del primo caso prova. . . . . . . . . . . . . . . . . 122

6.3 Mesh di riferimento. . . . . . . . . . . . . . . . . . . . . . . . 123

ELENCO DELLE FIGURE VI

6.4 Confronto della soluzione staggered con quella ottenuta con il

programma Comsol. . . . . . . . . . . . . . . . . . . . . . . . 124

6.5 1 modo di vibrare della trave doppiamente incastrata. . . . . 125

6.6 Confronto della freccia massima ottenuta con codici 0 e 1. . . 127

6.7 Confronto del potenziale di riferimento ottenuto con codici 0

e 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

6.8 Deformata qualitativa della trave a doppio incastro. . . . . . 128

6.9 Confronto degli spostamenti in direzione x (spostamenti in µm).129

6.10 Confronto degli spostamenti in direzione y (spostamenti in µm).130

6.11 Confronto degli sforzi σx (in N/m2). . . . . . . . . . . . . . . 131

6.12 Confronto degli sforzi σy (in N/m2). . . . . . . . . . . . . . . 132

6.13 Confronto degli sforzi τxy (in N/m2). . . . . . . . . . . . . . . 133

6.14 Confronto del potenziale (in V olt). . . . . . . . . . . . . . . . 134

6.15 Confronto del campo elettrico in direzione x (in V/m). . . . . 135

6.16 Confronto del campo elettrico in direzione y (in V/m). . . . . 136

6.17 Suddivisione del dominio meccanico in 2 sottodomini. . . . . 137



e 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

6.20 Spostamenti e velocita del nodo 3 di interfaccia fra i domini

1 e 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

6.21 Confronto dei tempi di analisi con diverso numero di sottodo-

mini meccanici. . . . . . . . . . . . . . . . . . . . . . . . . . . 141

6.22 Suddivisione della parte elettrica in 4 sottodomini e della

parte meccanica in 2 sottodomini. . . . . . . . . . . . . . . . 142

6.23 Confronto della freccia massima ottenuta con codici 0 e 3 (2

sottodomini meccanici, 4 sottodomini elettrici). . . . . . . . . 143


e 3 (2 sottodomini meccanici, 4 sottodomini elettrici). . . . . 144


mini elettrici e meccanici. . . . . . . . . . . . . . . . . . . . . 144

6.26 Geometria del secondo caso prova. . . . . . . . . . . . . . . . 145

6.27 Confronto della soluzione staggered con quella ottenuta con il

programma Comsol. . . . . . . . . . . . . . . . . . . . . . . . 146

6.28 1 modo di vibrare della trave a mensola. . . . . . . . . . . . 147


ELENCO DELLE FIGURE VII


e 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

6.31 Spostamenti in direzione x e y (in µm). . . . . . . . . . . . . 149

6.32 Sforzi normali in direzione x e y (in N/m2). . . . . . . . . . . 150

6.33 Sforzi di taglio (in N/m2) e potenziale (in V olt). . . . . . . . 151

6.34 Campi elettrici in direzione x e y (in V/m). . . . . . . . . . . 152


mini meccanici. . . . . . . . . . . . . . . . . . . . . . . . . . . 154


mini elettrici e meccanici. . . . . . . . . . . . . . . . . . . . . 155

6.37 Trave doppiamente incastrata soggetta a differenza di poten-

ziale variabile nel tempo. . . . . . . . . . . . . . . . . . . . . . 157

6.38 Attuazione sinusoidale. . . . . . . . . . . . . . . . . . . . . . . 158

6.39 Forzante elettrica nel nodo di mezzeria della trave. . . . . . . 160

6.40 Spostamento verticale del nodo di mezzeria della trave. . . . . 160

6.41 Geometria del portale. . . . . . . . . . . . . . . . . . . . . . . 161

6.42 Deformata del portale sollecitato a carico distribuito variabile

sul montante di sinistra. . . . . . . . . . . . . . . . . . . . . . 162

6.43 Mesh e suddivisione in sottodomini delle casistiche implemen-

tate. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163

6.44 Confronto degli spostamenti in direzione x del nodo A. . . . . 165

6.45 Confronto dei potenziali nel nodo B. . . . . . . . . . . . . . . 165

6.46 Confronto degli spostamenti in direzione x (in µm). . . . . . 167

6.47 Confronto degli spostamenti in direzione y (in µm). . . . . . . 168

6.48 Confronto degli sforzi σx (in N/m2). . . . . . . . . . . . . . . 169

6.49 Confronto degli sforzi σy (in N/m2). . . . . . . . . . . . . . . 170

6.50 Confronto degli sforzi τxy (in N/m2). . . . . . . . . . . . . . . 171

6.51 Confronto del campo elettrico in direzione x (in V/m). . . . . 172

6.52 Confronto del campo elettrico in direzione y (in V/m). . . . . 173

6.53 Confronto del potenziale (in V olt). . . . . . . . . . . . . . . . 174

6.54 Spostamenti in direzione x del nodo C (2 sottodomini). . . . 175

6.55 Spostamenti in direzione x del nodo E (3 sottodomini). . . . 175

6.56 Spostamenti in direzione x del nodo C (4 sottodomini). . . . 176

6.57 Spostamenti in direzione x del nodo E (6 sottodomini). . . . 176

6.58 Pulsazioni dei primi due modi della trave a doppio incastro. . 178

6.59 Andamento della freccia massima nel tempo. . . . . . . . . . 179

ELENCO DELLE FIGURE VIII

6.60 Andamento della forza elettrostatica nodale di mezzeria nel

tempo. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 179

6.61 Trave a doppio incastro soggetta a doppia attuazione. . . . . 181

6.62 Andamento nel tempo dei due potenziali di attuazione. . . . . 182

6.63 Andamento della freccia massima nel tempo. . . . . . . . . . 183

D.1 Schema del codice di decomposizione 1 livello. . . . . . . . . 204



Elenco delle tabelle

4.1 Geometria e proprieta dei domini elettrico e meccanico del

caso 1D. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80

4.2 Discretizzazione temporale . . . . . . . . . . . . . . . . . . . . 83

4.3 Periodo di oscillazione ottenuto dalla simulazione numerica. . 85

4.4 Confronto dei tempi analisi monolitico-staggered. . . . . . . . 85

6.1 Proprieta dei materiali delle parti meccanica ed elettrica. . . 122

6.2 Parametri dell’analisi dinamica. . . . . . . . . . . . . . . . . . 123

6.3 Confronto con la soluzione statica. . . . . . . . . . . . . . . . 124


6.5 Confronto dei tempi analisi staggered -decomposizione 1livello. 128

6.6 Confronto dei tempi analisi staggered -decomposizione 2livello

(2 sottodomini). . . . . . . . . . . . . . . . . . . . . . . . . . 139


(2 sd meccanici, 4 sd elettrici). . . . . . . . . . . . . . . . . . 142

6.8 Riassunto dei tempi analisi ottenuti con i vari algoritmi (trave

a doppio incastro). . . . . . . . . . . . . . . . . . . . . . . . . 143

6.9 Confronto con la soluzione statica. . . . . . . . . . . . . . . . 145


6.11 Confronto dei tempi di analisi staggered -decomposizione 1livello.153


(2 sottodomini). . . . . . . . . . . . . . . . . . . . . . . . . . 153


(2 sd meccanici, 4 sd elettrici). . . . . . . . . . . . . . . . . . 154

6.14 Riassunto dei tempi analisi ottenuti con i vari algoritmi (trave

a mensola). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155


6.16 Stima dei valori di picco della forza elettrica. . . . . . . . . . 159

IX

ELENCO DELLE TABELLE X

6.17 Stima dei valori di picco della freccia massima. . . . . . . . . 159

6.18 Confronto dei tempi di analisi staggered -decomposizione 1livello

(attuazione sinusoidale). . . . . . . . . . . . . . . . . . . . . . 160

6.19 Numero di nodi dei sottodomini meccanici. . . . . . . . . . . 164

6.20 Numero di elementi dei sottodomini meccanici. . . . . . . . . 164

6.21 Numero di nodi ed elementi del dominio elettrico. . . . . . . . 164


6.23 Valori massimi e minimi delle varie grandezze, valutate nel

nodo A, nell’istante t = 2 · 10−7 sec. . . . . . . . . . . . . . . 166

6.24 Confronto tra i tempi di analisi e relativo guadagno rispetto

all’algoritmo staggered. . . . . . . . . . . . . . . . . . . . . . . 177

6.25 Periodo di oscillazione ottenuto numericamente. . . . . . . . . 180


(trave smorzata). . . . . . . . . . . . . . . . . . . . . . . . . . 180


(doppia attuazione). . . . . . . . . . . . . . . . . . . . . . . . 182

C.1 Caratteristiche dei processori. . . . . . . . . . . . . . . . . . . 199

Capitolo 1

Introduzione

Il presente lavoro di tesi mette in relazione due grandi tematiche: da una

parte i metodi di decomposizione dei domini, intesi come particolari tecniche

di modellazione ad elementi finiti; dall’altra, il problema elettro-meccanico

accoppiato, appartenente alla categoria dei cosiddetti multi-field problems

(o problemi a piu campi), nei quali occorre tener conto dell’interazione

fra diversi fenomeni fisici; per questo motivo, si parla anche di problemi

accoppiati.

Il problema elettro-meccanico e tipico dei Micro-Sistemi Elettro-Mecca-

nici (MEMS), termine col quale vengono indicati dispositivi di varia natura

(meccanici, elettrici ed elettronici), integrati in forma altamente miniatu-

rizzata, la cui tecnologia e adottata negli ambiti applicativi piu vari; essi

trovano impiego come sensori, accelerometri, testine per stampanti, ecc.. A

causa delle loro piccole dimensioni, i carichi gravitazionali sono trascurabi-

li nei confronti delle forze elettriche, elastiche e adesive, preponderanti alla

microscala; queste di fatto determinano il comportamento meccanico di quei

dispositivi. La figura 1.1 mostra una particolare tipologia di micro-sistemi,

i comb-drives, costituiti da due serie di elettrodi, l’una fissa e l’altra mobile,

a costituire come due pettini, tra i quali e applicata una differenza di po-

tenziale; le conseguenti forze elettrostatiche che nascono in seguito a quella

attuazione sono proporzionali alla variazione di capacita fra le due serie di

elettrodi.

Nell’ambito della progettazione dei micro-sistemi e estremamente im-

portante riuscire a simulare l’accoppiamento elettro-meccanico in modo suf-

ficientemente preciso, per poter prevedere il comportamento del dispositivo

in esercizio e dimensionarlo di conseguenza; analogamente a quanto avviene

1

2

Figura 1.1: Esempio di micro-sistema.

alla macroscala, e necessario mantenersi lontani da una situazione di insta-

bilita che prende in questo caso il nome di pull-in, in corrispondenza della

quale gli elettrodi fra cui e imposta la differenza di potenziale vengono a

contatto fra loro, compromettendo il funzionamento del dispositivo.

La simulazione e quindi una fase estremamente importante nella pro-

gettazione dei micro-sistemi: essa avviene sfruttando dei metodi approssi-

mati, quale puo essere quello degli elementi finiti, e affidando il calcolo agli

elaboratori. Questi, avendo un limite fisico, necessitano di tempi di ela-

borazione che possono essere molto lunghi nel caso di analisi dinamiche,

problemi a geometria complessa, o di problemi fortemente accoppiati, come

quello elettro-meccanico.

Il punto di partenza del presente lavoro e proprio l’idea di applicare

i metodi di decomposizione dei domini, originariamente proposti per pro-

blemi meccanici sia statici che dinamici (per i quali e stata dimostrata la

maggiore efficienza nel calcolo rispetto ad un approccio di tipo monoliti-

co), al problema elettro-meccanico, con l’obiettivo di verificare sia l’entita

del vantaggio computazionale che essi possono comportare, sia l’accuratezza

della soluzione.

Si mostrera come la decomposizione possa essere applicata a quel tipo

di problema in corrispondenza di piu livelli, a partire dalla naturale suddi-

visione secondo le fisiche del problema fino a giungere alla partizione vera e

3

propria dei due domini meccanico ed elettrico. Occorrera tener conto delle

diverse caratteristiche dei due problemi (il meccanico risolto in dinamica,

l’elettrico in statica) e applicare di conseguenza gli algoritmi di decomposi-

zione piu adatti.

Nel capitolo 2 vengono passati in rassegna i principali metodi di decom-

posizione cosı come sono stati originariamente proposti, fornendo una pano-

ramica degli algoritmi di provata validita. Dopo un breve accenno all’approc-

cio per sottostrutture, precursore dei moderni metodi di decomposizione, si

procede ad introdurre i due grandi metodi proposti nei primi anni ’90, dai

quali sono stati dedotti tutti i successivi: il Balancing Domain Decomposi-

tion (BDD, [9]) e il Finite Element Tearing and Interconnecting (FETI, [1]).

Infine, viene affrontato l’algoritmo di Gravouil e Combescure per problemi

di dinamica strutturale, proposto recentemente e dedotto dal metodo FETI

([8]).

Nel capitolo 3 vengono dedotte le equazioni che governano il problema

elettrico, al quale viene associata una formulazione ad elementi finiti. Ven-

gono quindi messe in evidenza le caratteristiche dell’accoppiamento elettro-

meccanico; sulla base di queste, si procede alla sintesi dei due sistemi elettri-

co e meccanico, inizialmente disaccoppiati, a formare il sistema che governa

il problema accoppiato. Un esempio monodimensionale consente di com-

prendere meglio tale accoppiamento e di affrontare il tema dell’instabilita

elettro-meccanica (pull-in) sia nel caso statico che in quello dinamico.

Dedotto il sistema elettro-meccanico, il capitolo 4 presenta i meto-

di di risoluzione classici (monolitico e staggered) per la risoluzione di tale

problema. Dopo una trattazione principalmente teorica, viene presentato

un esempio monodimensionale risolto con le due procedure gia consolida-

te; per tale problema sono stati opportunamente implementati dei codici in

Matlabr, che hanno consentito un confronto fra i due diversi metodi.

Nel capitolo 5 vengono analizzate le modalita con cui la decomposizio-

ne dei domini puo essere applicata al problema elettro-meccanico. Partendo

dalla procedura di tipo staggered, si introduce innanzitutto una suddivisione

in base alle fisiche del problema, applicando l’algoritmo di Gravouil e Com-

bescure (GC) proposto al capitolo 2. Il nuovo algoritmo ottenuto e indicato

con il nome di decomposizione di 1 livello. Il passaggio successivo consiste

nell’introduzione della decomposizione del dominio meccanico, anche in que-

sto caso nello spirito di GC, dando vita al cosiddetto algoritmo di 2 livello.

4

Infine, l’ultimo passo consiste nell’introdurre la decomposizione della parte

elettrica, nella quale, come osservato nel capitolo 3, i problemi sono risolti in

statica. E quindi inevitabile l’adozione del metodo FETI, che viene in questa

sede esteso al problema elettrico secondo una nuova formulazione, dedotta

tenendo conto dell’analogia col problema meccanico. Il metodo FETI viene

anche applicato alla risoluzione del problema pseudo-meccanico, con il quale

viene aggiornata la configurazione del dominio elettrico sulla base degli spo-

stamenti indotti dalla parte meccanica. Questo terzo e ultimo algoritmo e

indicato con il nome di decomposizione di 3 livello.

Infine, nel capitolo 6 vengono analizzati alcuni casi semplici ai quali

sono applicati gli algoritmi proposti; questi sono confrontati con la procedura

di tipo staggered. Per poter effettuare queste analisi sono stati implementati

ex novo dei codici in Matlab per la procedura di tipo staggered e per i tre

nuovi algoritmi proposti (1, 2 e 3 livello); in realta, per la realizzazione

dell’ultimo codice, si e preso spunto dall’implementazione del metodo FETI

per lo studio di piastre, realizzato in un precedente lavoro di tesi ([17]). Quel-

le subroutines sono state rimaneggiate e adattate ai casi pseudo-meccanico

ed elettrico. Gli esempi proposti si distinguono in due categorie: innanzitut-

to si studiano dei casi prova, ai quali sono applicati tutti i tipi di algoritmi

affrontati nel corso dei precedenti capitoli; questo consentira la validazione

dei codici implementati. In secondo luogo, sono proposti alcuni casi parti-

colari nei quali sono messi a confronto gli algoritmi proposti e dove viene

dimostrato non solo il vantaggio computazionale derivante dall’introduzione

della decomposizione a piu livelli, ma anche l’equivalenza con gli algoritmi

classici di risoluzione del problema elettro-meccanico.

Capitolo 2

La decomposizione dei

domini: stato dell’arte

2.1 Introduzione

A partire dagli anni ’80 del secolo scorso grande importanza e stata attribui-

ta allo sviluppo di metodi efficienti per il calcolo della soluzione nei problemi

ad elementi finiti. Nelle simulazioni, infatti, i domini oggetto dell’analisi

possono risultare complessi per geometria o di notevoli dimensioni, oppure

possono presentare diversi comportamenti fisici (domini meccanici, elettrici,

. . . ); tutti questi aspetti, in genere, comportano un onere di calcolo eccessivo.

Nell’ambito delle tecniche che cercano di porre rimedio a questo problema

si situano i metodi di decomposizione dei domini ; ciascuna delle parti in cui

quelli sono suddivisi puo essere analizzata separatamente tenendo conto allo

stesso tempo dell’effetto-azione esercitata dalle parti adiacenti. In particola-

re, e necessario esprimere una condizione di interfaccia o di interazione fra le

parti, che varia a seconda del tipo di problema affrontato (es.: nei problemi

di dinamica strutturale, nei quali i domini sono solidi meccanici, tale condi-

zione corrisponde alla continuita del campo di spostamento e all’uguaglianza

delle trazioni all’interfaccia).

Il matematico tedesco Hermann Schwarz (1843-1921) e da molti con-

siderato il padre della decomposizione dei domini: risale infatti al 1869 un

articolo in cui egli propose un nuovo metodo per la risoluzione di equazioni

differenziali alle derivate parziali su un dominio complesso, costituito dalla

sovrapposizione di un disco e di un rettangolo (figura 2.1); con quell’articolo

furono poste le basi matematiche dei metodi di decomposizione. In realta,

5

2.1 Introduzione 6

Figura 2.1: Problema di Schwarz ([10]).

e solo un secolo piu tardi, a partire dagli anni ’80 del secolo scorso, con lo

sviluppo dei moderni processori, che ci si rese effettivamente conto dei van-

taggi che l’accoppiamento fra processamento parallelo su piu unita di calcolo

e metodi di decomposizione poteva garantire; con queste tecniche e possibile

risolvere problemi indipendenti su diversi processori, provvedendo contem-

poraneamente ad uno scambio di informazioni in modo sincronizzato. Un

approccio di questo tipo rappresenta oggi un aspetto essenziale, di fronte alla

necessita di effettuare simulazioni numeriche di fenomeni fisici sempre piu

complessi (non-linearita, accoppiamenti, . . . ), in grado di sfruttare al meglio

le potenzialita della moderna tecnologia.

Gia a partire dagli anni ’60 gli ingegneri aeronautici americani svilup-

parono un metodo per lo studio di strutture con elevato numero di gradi di

liberta, quali appunto gli aeromobili. Tale metodo comporta la suddivisione

del dominio di analisi in sottostrutture, per ciascuna delle quali si effettua

la condensazione dei gradi di liberta interni in modo tale da ottenere un si-

stema assemblato esclusivamente sui gradi di liberta di interfaccia (e quindi

ridotto).

E pero solo all’inizio degli anni ’90 che vengono sviluppati due me-

todi di decomposizione dei domini dai quali sono stati successivamente de-

dotti diversi algoritmi, generalmente classificati nelle seguenti due grandi

categorie:

• metodi primari (BDD: Balancing Domain Decomposition);

• metodi duali (FETI: Finite Element Tearing and Interconnecting).

Contrariamente al problema di Schwarz, i due tipi di metodi valgono per

domini non sovrapposti (non-overlapping) e si distinguono tra di loro per

2.2 Approccio per sottostrutture 7

le incognite del problema di interfaccia: spostamenti nei primi, trazioni nei

secondi.1

I metodi di decomposizione presentano vantaggi di diverso tipo: al di

la della naturale riduzione della dimensione dei sistemi risolti (legata al fatto

di considerare N +1 sistemi di dimensioni ridotte - con N numero dei sotto-

domini - in luogo di un unico grande sistema monolitico), metodi di questo

tipo consentono anche di considerare allo stesso tempo comportamenti fisici

differenti, potendo associare a ciascun sottodominio equazioni che governano

un determinato fenomeno fisico. Un esempio, che verra affrontato in segui-

to, e il cosiddetto problema elettro-meccanico accoppiato, che si presenta in

numerose applicazioni tra cui i Micro Sistemi Elettro-Meccanici (MEMS).

2.2 Approccio per sottostrutture

Come accennato in precedenza, la caratteristica principale dell’approccio per

sottostrutture e quella di ottenere una riduzione del numero dei gradi di li-

berta applicando il concetto di condensazione statica alle incognite nodali

interne di ciascuna sottostruttura in cui si immagina suddivisa la struttura

di partenza (figura 2.2); il risultato e la riduzione del problema monolitico

ad essa associato ad una serie di problemi di minore dimensione. Il sistema

algebrico discretizzato monolitico e il seguente:

KU = P (2.1)

in cui la matrice di rigidezza K e il vettore dei carichi nodali equivalenti

P hanno rispettivamente ordine e dimensione pari al numero totale n dei

gradi di liberta della struttura. Si immagini ora di operare una suddivisio-

ne in N parti. Per ciascuna di esse, il sistema di equilibrio viene scritto come:

KiUi = Pi +Pi0 i = 1, . . . , N (2.2)

in cui si mettono in evidenza i carichi nodali di interfaccia Pi0, cioe le forze

interne che la sottostruttura i-esima scambia con quelle ad essa adiacenti.

1Si fa riferimento al problema meccanico, per il quale questi metodi sono stati

originariamente proposti.


Figura 2.2: Esempio di suddivisione in sottostrutture.

Si riscrive il sistema (2.2) introducendo una partizione fra gradi di liberta

esterni (E) ed interni (I) (si omette l’indice i della sottostruttura per una

maggiore chiarezza della trattazione):

[

KEE KEI

KIE KII

][

UE

UI

]

=

[

PE

PI

]

+

[

P0E

0

]

(2.3)

Si osserva che P0I = 0 per definizione stessa di forze di interfaccia.

Dalla seconda delle (2.3) si esprime UI :

UI = K−1II (PI −KIEUE) (2.4)

che, sostituita nella prima, restituisce:

(KEE −KEIK−1II KIE)

︸︷︷︸

K

UE = (PE −KEIK−1II PI)

︸︷︷︸

P

+P0E

KiUiE = Pi +Pi

0Ei = 1, . . . , N (2.5)

Si osservi che nel calcolo di Ki e Pi l’operazione piu onerosa e l’inversione

della sottomatrice KiII , da effettuarsi per ogni sottostruttura: l’ordine di

ciascuna di quelle matrici e comunque ridotto e pari al numero di gradi di

liberta interni della corrispondente sottostruttura.


Si procede come nell’analisi ad elementi finiti usuale, considerando le

(2.5) come equazioni di equilibrio associate a singoli elementi finiti.2 L’ener-

gia potenziale totale del sistema puo essere espressa in modo immediato (si

veda [2]):

Π =N∑

i=1

(1

2UiT

E KiUiE −UiT

E Pi −

UiT

E Pi0E

)

(2.6)

in cui l’ultimo contributo, associato alle forze di interazione, e nullo per il

principio di azione-reazione: in particolare, il lavoro di quelle forze lungo

ogni linea di interfaccia da un duplice contributo uguale ed opposto. Si in-

troduce il vettore globale UE dei gradi di liberta del sistema, comprensivo

degli spostamenti sul contorno esterno della struttura e sulle interfacce; le

matrici di connettivita consentono di esprimere gli spostamenti UiE di cia-

scuna sottostruttura in funzione di UE :

UiE = LiUE (2.7)

Sostituendo la (2.7) nella (2.6) si ottiene:

Π =1

2UT

E

( N∑

i=1

LiT KiLi)

︸︷︷︸

K

UE −UTE

( N∑

i=1

LiTPi)

︸︷︷︸

P

(2.8)

Imponendo la condizione di stazionarieta dell’energia potenziale totale, si

ricava in definitiva:

δΠ = 0, ∀ δUE = 0 ⇒ KUE = P (2.9)

Il sistema (2.9) ha dimensione pari al numero di incognite nodali di contorno

e interfaccia. Il sistema monolitico di partenza (2.1) e quindi ricondotto a

N sottoproblemi di dimensione ridotta, a cui si aggiunge la risoluzione del

problema (2.9).

2L’analogia infatti e evidente: la suddivisione della struttura globale in sottostrutture

e del tutto analoga alla discretizzazione ad elementi finiti di un generico dominio.

2.3 Metodi di decomposizione dei domini 10

Una volta determinato il vettore UE , e possibile risalire a UiE con la

(2.7) e agli spostamenti nodali interni attraverso la (2.4). 3

Nonostante sia di qualche decennio antecedente ai metodi di decom-

posizione veri e propri, l’approccio per sottostrutture possiede gia alcune

delle caratteristiche degli algoritmi proposti nei primi anni ’90: suddivisione

della struttura di partenza in piu parti, sostituzione del sistema monolitico,

computazionalmente troppo oneroso, con piu sottoproblemi di dimensione

ridotta e importanza attribuita all’interfaccia intesa come quella zona, fisi-

camente individuabile, attraverso cui avviene lo scambio di informazioni fra

sottodomini adiacenti.

2.3 Metodi di decomposizione dei domini

Come anticipato i metodi di decomposizione dei domini si suddividono nelle

due grandi categorie di metodi primari e metodi duali a seconda delle inco-

gnite assunte nel problema di interfaccia, spostamenti nei primi e trazioni nei

secondi. I metodi BDD ([9]) e FETI ([1]) sono i primi, in ordine cronologi-

co, ad essere stati proposti e appartengono rispettivamente alla prima e alla

seconda categoria. Essi sono caratterizzati da un’interpretazione meccanica

molto semplice e da un problema di interfaccia che, in vista della paralle-

lizzazione su piu processori, viene, in genere, affrontato con un risolutore di

tipo iterativo; ogni iterazione, d’altra parte, richiede la risoluzione di un pro-

blema ad elementi finiti che puo avvenire tramite un risolutore diretto. In

questo modo, i metodi di decomposizione combinano l’utilizzo di risolutori

diretto e iterativo, cercando di ottenere dal primo robustezza e dal secondo

risparmio computazionale.

I metodi primari e duali sono stati estesi ai problemi eterogenei e

di elasticita del quarto ordine (piastre e gusci); si parla in questo caso di

algoritmi BDDC e FETIDP ([10]), che oggi risultano essere efficienti tanto

quanto i due metodi originari ([9],[1]).

Ulteriori algoritmi sono stati proposti per ambiti diversi da quello pret-

tamente statico-meccanico; tra questi il transitorio dinamico ed i problemi

multi-fisici (mezzo poroso, fluido incomprimibile, contatto,. . . ; si veda ad

esempio [15]).

3Si osservi che la matrice K−1II e gia stata calcolata per l’assemblaggio della matrice K.


Dato lo stretto legame fra metodi primari e duali, sono stati proposti

anche algoritmi che cercano di fondere le due categorie in una sola. E il caso

dell’approccio ibrido ([10]), che consente di specificare per ciascun grado di

liberta di interfaccia il tipo di incognita (spostamento o trazione); nel caso

dei problemi multi-fisici un approccio di questo tipo consente di definire degli

opportuni risolutori che tengano conto del diverso comportamento fisico che

si presenta. Gli approcci di tipo misto invece ricercano una combinazione

lineare di spostamenti di interfaccia e di sforzi a seconda di una rigidezza

fittizia/artificiale; in particolare gli approcci misti ricadono nella prima o

nella seconda categoria nel caso in cui la rigidezza sia infinita (primari) o

nulla (duali). I metodi misti presentano vantaggi soprattutto per la maggio-

re semplicita con cui vengono affrontate interfacce complesse (per esempio

in presenza di fenomeni di contatto o di attrito) rispetto ai metodi originari.

2.3.1 Formulazione del problema di interfaccia

Problema di riferimento

Si consideri un generico dominio meccanico Ω in Rn, con n = 1, 2, 3 (figura

2.3). Le equazioni che governano il problema elasto-meccanico, usando la

notazione tensoriale compatta proposta in [10], sono:

div(σ) + f = 0 in Ω

σ = Dε in Ω

ε = 12

(

∇u+ (∇u)T)

in Ω

σn = g su ∂gΩ

u = u0 su ∂uΩ

(2.10)

Le equazioni del sistema (2.10), come e noto, esprimono rispettivamente

l’equilibrio, il legame costitutivo (assunto elastico lineare) e la congruenza,

oltre alle condizioni al contorno naturali (equilibrio al contorno) ed essenziali

(di congruenza cinematica fra spostamento e cedimento sulla parte vincola-

ta) rappresentate dalle ultime due equazioni.


Figura 2.3: Problema di riferimento ([10]).

Decomposizione in due sottodomini

Si consideri, in primo luogo, una suddivisione del dominio Ω di partenza in

due parti Ω1 e Ω2 come in figura 2.4. L’interfaccia Υ tra le due parti e

definita nel modo seguente:

Υ = ∂Ω1

⋂

∂Ω2 (2.11)

Le equazioni del problema elastico (2.10) possono essere scritte per ciascuna

delle due sottostrutture:

s = 1, 2

div(σs) + fs = 0 in Ωs

σs = Dsεs in Ωs

εs =12

(

∇us + (∇us)T)

in Ωs

σsns = gs su ∂gΩ⋂∂Ωs

us = u0s su ∂uΩ⋂∂Ωs

(2.12)

I due sistemi (2.12), scritti rispettivamente per Ω1 e Ω2, non sono sufficienti

a sostituire il sistema (2.10); essi devono infatti essere corredati da condizio-

ni all’interfaccia, che esprimano la connessione fra i due sottodomini:

u1 = u2 su Υ (2.13)

σ1n1 = σ2n2 su Υ (2.14)

Le (2.13) e (2.14) esprimono rispettivamente la continuita del campo di spo-

stamento e l’equilibrio degli sforzi (principio di azione-reazione) all’interfac-


Figura 2.4: Partizione in due sottodomini ([10]).

cia.

Le equazioni (2.12), (2.13) e (2.14) sono equivalenti al sistema monolitico

(2.10).

Decomposizione in N sottodomini

Si immagini ora di suddividere il dominio Ω in N parti indicate con Ωs.

Data la maggiore complessita della decomposizione geometrica, occorre in-

trodurre i concetti di interfaccia fra due sottodomini, interfaccia totale di un

sottodominio e interfaccia della struttura globale, definiti come:

Υi,j = Υj,i = ∂Ωi⋂∂Ωj

Υs =⋃

j Υs,j

Υ =⋃

sΥs

(2.15)

In genere ci si riferisce alle Υs e Υ come interfacce locali e interfaccia globale

rispettivamente; Υi,j rappresenta invece l’interfaccia di comunicazione tra i

due sottodomini i e j.

Nel caso di suddivisione in N parti possono formarsi i cosiddetti cros-

spoints, cioe nodi condivisi da piu di due sottodomini. La loro esistenza

rende necessaria la definizione di due nuove entita: la cosiddetta interfac-

cia geometrica, coincidente con Υ, e l’interfaccia di connettivita, costituita

dall’insieme di tutte le interfacce associate alle varie coppie di sottodomini

adiacenti (Υi,j , con 1 ≤ i < j ≤ N). Ciascuna delle due categorie di metodi

di decomposizione fa riferimento ad una di queste due tipologie di intefaccia,

per questo in genere ci si riferisce all’interfaccia geometrica come interfaccia

primaria e a quella di connettivita come interfaccia duale (figura 2.5).


Figura 2.5: Interfaccia primaria e duale ([10]).

Definizione degli operatori booleani

Per poter scrivere in modo compatto le equazioni associate ai vari sottodo-

mini e per seguire allo stesso tempo una trattazione rigorosa che sia la piu

generale possibile (come quella proposta in [10]), e necessario definire diversi

operatori.

Innanzitutto si definisce l’operatore ts che consente di estrarre dall’in-

sieme dei gradi di liberta del sottodominio Ωs quelli della sua interfaccia.

Il suo trasposto consente l’operazione inversa, cioe il passaggio dai gradi di

liberta di interfaccia a quelli dell’intero sottodominio.

In secondo luogo, si definisce un operatore che consente il passaggio

di dati da un sottodominio a quelli adiacenti; a seconda che il metodo di

decomposizione (e quindi l’interfaccia) sia primario o duale si definiscono gli

operatori di assemblaggio corrispondenti, indicati rispettivamente con As e

As. Il primo e di tipo booleano, il secondo e, invece, dotato di segno: in

quest’ultimo caso, infatti, se ad un grado di liberta corrisponde 1 da un

lato dell’interfaccia, al corrispondente grado di liberta pensato come appar-

tenente all’altro sottodominio e associato invece −1. As e As godono delle

seguenti proprieta:


∑

sAsATs = 0

ATs As = IΥs

ATs As = diag(molteplicita−1)Υs

AsATs =

∣∣∣∣∣

I su Υs

0 altrove

(2.16)

Nella terza delle (2.16) si e indicato con molteplicita il numero di sottodo-

mini a cui appartiene il generico punto di interfaccia del dominio s. Dalle

precedenti definizioni deriva il seguente operatore booleano composto:

Bs = Asts (2.17)

L’esempio che segue ([10]) consente di comprendere meglio il significato di

cio che e stato definito e che puo sembrare, di primo acchito, eccessivamente

astratto.

Esempio

Si consideri la partizione di figura 2.6. Considerando la numerazione ripor-

tata, si scrivono gli operatori booleani in questo caso particolare:

t1 =

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

t2 =

0 0 1 0 0

0 0 0 1 0

0 0 0 0 1

t3 =

1 0 0 0

0 1 0 0

0 0 1 0

A1 =

0 0 0

0 1 0

1 0 0

0 0 1

A2 =

1 0 0

0 0 1

0 0 0

0 1 0

A3 =

0 0 1

0 0 0

1 0 0

0 1 0


Figura 2.6: Esempio di suddivisione in tre sottodomini ([10]).

A1 =

0 0 0

0 1 0

1 0 0

0 0 0

0 0 1

0 0 1

A2 =

1 0 0

0 0 −1

0 0 0

0 1 0

0 −1 1

0 0 0

A3 =

0 0 −1

0 0 0

−1 0 0

0 −1 0

0 0 0

0 −1 0

Si indichino con ns, nbs e ni

s rispettivamente i gradi di liberta totali, di

interfaccia ed interni del sottodominio s-esimo.4 Si definiscano poi nΥ e

nΥi,ji gradi di liberta associati rispettivamente all’interfaccia primaria e

duale (quattro e sei per l’esempio che si sta considerando; la differenza,

4In questo esempio si assume, per semplicita, che a ciascun nodo corrisponda un unico

grado di liberta.


come spiegato in precedenza, consiste nel fatto che la seconda tipologia di

interfaccia tiene conto degli accoppiamenti fra i vari sottodomini, ne segue

che nΥ e nΥi,jrisultano in genere diversi). Si deduce che:

- gli operatori ts hanno dimensioni nbs × ns;

- gli operatori As hanno dimensioni nΥ × nbs;

- gli operatori As hanno dimensioni nΥi,j× nb

s.

Gli operatori booleani composti Bs:

B1 =

0 0 0 0 0

0 0 0 1 0

0 0 1 0 0

0 0 0 0 0

0 0 0 0 1

0 0 0 0 1

B2 =

0 0 1 0 0

0 0 0 0 −1

0 0 0 0 0

0 0 0 1 0

0 0 0 −1 0

0 0 0 0 0

B3 =

0 0 −1 0

0 0 0 0

−1 0 0 0

0 −1 0 0

0 0 0 0

0 −1 0 0

Gli operatori booleani composti (di seguito booleani) hanno quindi dimen-

sioni nΥi,j× ns, consentono cioe il passaggio dalla numerazione interna del

sottodominio s-esimo a quella dei gradi di liberta (nodi in questo caso) del-

l’interfaccia duale.

2.3.2 Il metodo BDD

Il Balancing Domain Decomposition method e stato proposto da J.Mandel

in [9]; esso appartiene alla categoria dei metodi primari, in quanto adotta gli

spostamenti come incognite del problema di interfaccia.

Il punto di partenza per la formulazione del metodo e dato dalla scrit-

tura dei sistemi algebrici discretizzati associati a ciascuno degli N sottodo-

mini Ωs in cui si immagina suddiviso il dominio Ω di partenza:


KsUs = fs + λs s = 1, . . . , N (2.18)

in cui sono messe in evidenza le forze di interazione λs di ciascun sottodomi-

nio. Operando una partizione fra gradi di liberta interni (i) e di interfaccia

(b), il sistema (2.18) viene riscritto come:

[

Kiis Kib

s

Kbis Kbb

s

][

Uis

Ubs

]

=

[

f isf bs + λ

bs

]

(2.19)

dove le forze di interazione λ agiscono ovviamente nei soli nodi di interfaccia.

Con un procedimento del tutto analogo a quello utilizzato nell’approccio per

sottostrutture (si veda il paragrafo 2.2), il sistema (2.19) puo essere riscritto

in termini dei soli spostamenti di interfaccia. Dalla prima delle (2.19) si

ricava:

Uis = Kii−1

s (f is −Kibs U

bs) (2.20)

Sostituendo nella seconda:

(Kbbs −Kbi

s Kii−1

s Kibs )

︸︷︷︸

Sps

Ubs = (f bs −Kbi

s Kii−1

s f is)︸︷︷︸

bps

+λbs

SpsUbs = bps + λ

bs s = 1, . . . , N (2.21)

Procedendo ancora una volta come gia fatto nell’approccio per sottostruttu-

re, si esegue l’assemblaggio delle N equazioni (2.21), eliminando in tal modo

il contributo delle forze di interfaccia; cio che si ottiene e la formulazione

primaria del problema di interfaccia:

SpUb = bp (2.22)

in cui Sp rappresenta il complemento primario globale di Schur della strut-

tura decomposta, avente espressione:

Sp =N∑

s=1

LsTSpsLs (2.23)


con Ls le consuete matrici di connettivita.

Il sistema (2.22) puo essere risolto direttamente attraverso l’inversione

del complemento primario di Schur, le cui dimensioni sono pari al numero

totale di incognite nodali di interfaccia e che risulta essere in generale poco

sparso, essendo ciascun grado di liberta connesso a gradi di liberta apparte-

nenti agli stessi sottodomini. Per questo motivo Mandel propose di risolvere

quel sistema in modo iterativo, in modo tale da non dover calcolare l’inversa

diretta di Sp. In particolare, il risolutore utilizzato e il metodo del gradiente

coniugato precondizionato (si veda l’appendice ), che riduce la risoluzione

del sistema (2.22) ad una serie di prodotti matrice-vettore.

2.3.3 Il metodo FETI

Formulazione del metodo

Il Finite Element Tearing and Interconnecting method fu originariamente

proposto da Farhat e Roux in [1] per il problema elasto-statico, ma risulta

comunque applicabile a tutti i problemi governati da equazioni differenziali

alle derivate parziali, trattati mediante il metodo degli elementi finiti.

Come e noto, il problema elastico puo essere riformulato con riferimento alla

proprieta di stazionarieta del funzionale energia potenziale totale in corri-

spondenza della soluzione. La EPT per il generico solido elastico Ω soggetto

a forze F sul proprio volume e a trazioni f sul contorno non vincolato ΓF e

definita come segue ([2]):

Π(u) =1

2

∫

ΩεT d ε dV −

∫

ΩFT u dV −

∫

ΓF

fT u dΓF (2.24)

in cui il simbolo 2 indica il rispetto delle equazioni di congruenza.

Si immagini ora di suddividere il dominio Ω in due parti, ciascuna

delle quali puo essere considerata come indipendente dall’altra, mettendo

in evidenza le forze di interazione che esse si scambiano in corrispondenza

dell’interfaccia Γ1−2; tenendo conto della (2.24) e considerando il contribu-

to aggiuntivo delle forze di interazione, indicate con λ, l’energia potenziale

totale per i due sottodomini Ω1 e Ω2 risulta:


Π(u1) =1

2

∫

Ω1

εT1 d1 ε1 dV −

∫

Ω1

FT1 u1 dV −

−∫

ΓF1

fT1 u1 dΓ−∫

Γ1−2

λT1 u1 dΓ (2.25)

Π(u2) =1

2

∫

Ω2

εT2 d2 ε2 dV −

∫

Ω2

FT2 u2 dV −

−∫

ΓF2

fT2 u2 dΓ−∫

Γ1−2

λT2 u2 dΓ (2.26)

Per ripristinare la continuita dell’intero dominio Ω e necessario imporre le

seguenti condizioni:

u1 = u2 su Γ1−2 (2.27)

λ1 + λ2 = 0 su Γ1−2 (2.28)

Le (2.27) e (2.28) rappresentano le condizioni di continuita del campo di

spostamento e l’equilibrio delle trazioni sull’interfaccia Γ1−2.

Il problema elastico associato al solido Ω e quindi ricondotto alla determina-

zione del punto di stazionarieta dei funzionali (2.25) e (2.26), tenendo conto

del vincolo costituito dalle (2.27) e (2.28). A tal fine si definisce un nuovo

funzionale lagrangiano nel quale compare anche la variabile vettoriale λ che

assume il ruolo di moltiplicatore di Lagrange:

L(u1, u2,λ) = Π(u1) + Π(u2) +

∫

Γ1−2

λT (u1 − u2) dΓ (2.29)

La soluzione e ottenuta cercando il punto di stazionarieta del funzionale Lal variare di u1, u2,λ. Si puo scrivere:

δL,δu1= δΠ(u1) +

∫

Γ1−2

λT δu1 dΓ = 0 ∀ δu1 (2.30)

δL,δu2= δΠ(u2)−

∫

Γ1−2

λT δu2 dΓ = 0 ∀ δu2 (2.31)

δL,δλ =

∫

Γ1−2

(u1 − u2)T δλ dΓ = 0 ∀ δλ (2.32)


Figura 2.7: Definizione dei modelli di spostamento e delle forze di interfaccia.

Soluzione ad elementi finiti

Per ciascuno dei due sottodomini Ω1 e Ω2 si introduce una discretizzazione

spaziale definendo un modello per il campo di spostamento sul singolo ele-

mento finito e attraverso l’utilizzo delle funzioni di forma (figura 2.7):

ue1 = N1U

e1 ∀ e = 1, . . . , ne

1 (2.33)

ue2 = N2U

e2 ∀ e = 1, . . . , ne

2 (2.34)

avendo indicato con Ue1 e Ue

2 gli spostamenti nodali del generico elemento e

appartenente alla prima o alla seconda parte; ne1 e ne

2 rappresentano invece

il numero di elementi in cui si immaginano suddivisi i due sottodomini.

Accanto alle (2.33) e (2.34) si definisce anche un modello per il vettore dei

moltiplicatori di Lagrange λ su ciascun elemento in cui e suddivisa l’inter-

faccia Γ1−2 (si assume che le mesh siano compatibili):

λe = NΓΛ

e ∀ e = 1, . . . , neΓ (2.35)

avendo indicato con Λe le forze di interfaccia nodali del generico elemento e

con neΓ il numero di elementi in cui e suddivisa l’interfaccia stessa.

Introducendo i modelli (2.33), (2.34) e (2.35) nelle condizioni di staziona-

rieta (2.30), (2.31) e (2.32), e tenendo conto della definizione degli operatori


booleani (2.17), si ottiene il seguente sistema:

K1U1 = P1 +BT1 Λ

K2U2 = P2 +BT2 Λ

B1U1 +B2U2 = 0

(2.36)

in cui Ks e Ps sono rispettivamente la matrice di rigidezza ed il vettore dei

carichi nodali equivalenti del sottodominio s-esimo; Λ rappresenta il vettore

delle forze nodali di interfaccia, di dimensione pari al numero di gradi di

liberta lagrangiani, cioe associati all’interfaccia duale. Le matrici di connet-

tivita booleane Bs consentono il passaggio dai gradi di liberta di interfaccia

a quelli del sottodominio s-esimo; per definizione esse sono dotate di segno

(sono matrici contenenti 0, 1,−1) che dipende dalla convenzione assunta. Si

osservi che, nel caso particolare di due sottodomini, il numero totale di gra-

di di liberta di interfaccia corrisponde a quello dei due sottodomini che vi

insistono, si ha cioe nΥ1,2 = nb1 = nb

2.

Se i sottodomini sono opportunamente vincolati, le corrispondenti ma-

trici di rigidezza sono non singolari; di conseguenza la soluzione del sistema

(2.36) puo essere ottenuta attraverso la loro diretta inversione: in parti-

colare, esprimendo U1 e U2 in funzione di Λ dalle prime due equazioni e

sostituendo nella terza, risulta:

U1 = K−11 (P1 +BT

1 Λ)

U2 = K−12 (P2 +BT

2 Λ)

(B1K−11 BT

1 +B2K−12 BT

2 )Λ = −(B1K−11 P1 +B2K

−12 P2)

(2.37)

Dalla terza equazione del sistema (2.37) e possibile ricavare direttamente i

moltiplicatori di Lagrange Λ e quindi, sostituendo questo vettore nelle prime

due, risalire ai campi di spostamento dei due sottodomini.

Estensione al caso di N sottodomini

Quanto ottenuto finora per la suddivisione del dominio Ω in due parti puo

essere facilmente generalizzato al caso di N sottodomini (per ipotesi non so-

vrapposti). Il sistema ottenuto a valle della discretizzazione spaziale risulta:

KsUs = Ps +BTs Λ s = 1, . . . , N

∑Ns=1BsUs = 0

(2.38)


Figura 2.8: Esempio di partizione con sottodominio fluttuante.

Analogamente al caso di due sottodomini (sistema (2.37)), se le matrici Ks

sono non singolari, la soluzione puo essere ottenuta per inversione diretta

delle stesse, previa determinazione dei moltiplicatori di Lagrange. Qualora

le matrici presentino singolarita e necessario adottare la particolare strategia

computazionale esposta nel seguito.

Singolarita locali

Dalla suddivisione del dominio Ω di partenza possono ottenersi dei sottodo-

mini non sufficientemente vincolati o completamente non vincolati, indicati

in genere come sottodomini fluttuanti ; l’esempio tipico e quello di una menso-

la in cui i sottodomini sono ottenuti da una partizione verticale della stessa:

il dominio Ω2, pensato come a se stante, non e vincolato. Si faccia riferimen-

to per semplicita al caso prima citato della mensola suddivisa in due parti

(figura 2.8): il sistema risolvente e costituito dal sistema (2.36). L’inversione

nel sistema (2.36) della matrice di rigidezza del dominio fluttuante Ω2, sin-

golare per le labilita presenti, non e piu possibile direttamente. E necessario

ricorrere al concetto di matrice pseudo-inversa, definita come quella matrice

K+2 tale da soddisfare alle seguenti condizioni ([1]):

K2K+2 K2 = K2

K+2 K2K

+2 = K+

2

(K2K+2 )

T = K2K+2

(K+2 K2)

T = K+2 K2

(2.39)

La seconda delle equazioni del sistema (2.37) puo essere quindi sostituita


dalla seguente:

U2 = K+2 (P2 +BT

2 Λ) +R2α (2.40)

dove R2 e una matrice di dimensioni (ni2+nb

2)×nr2 ed ha il significato fisico

di modi rigidi del dominio fluttuante Ω2; il vettore α di dimensione nr2 rap-

presenta invece una combinazione lineare dei precedenti. Ne segue che nr2

coincide con il numero di moti rigidi consentiti al sottodominio, al massimo

pari a 3 e 6 rispettivamente per problemi bidimensionali e tridimensionali.

Sostituendo quindi la (2.40) nel sistema (2.37) si ricava:

U1 = K−11 (P1 +BT

1 Λ)

U2 = K+2 (P2 +BT

2 Λ) +R2α

(B1K−11 BT

1 +B2K+2 B

T2 )Λ = −

(

B1K−11 P1 +B2(K

+2 P2 +R2α)

)

(2.41)

La scrittura del sistema (2.41) comporta l’introduzione di ulteriori incognite

contenute nel vettore α; le equazioni di (2.41) sono pertanto insufficienti per

la risoluzione del problema. Un’ulteriore condizione puo pero essere intro-

dotta considerando la seconda equazione delle (2.36):

K2U2 = P2 +BT2 Λ (2.42)

Essendo la matrice K2 simmetrica, l’equazione (2.42) ammette almeno una

soluzione se e solo se il suo secondo membro non ha componenti nello spazio

nullo di K2; esso deve cioe risultare ortogonale al nucleo della stessa matri-

ce.5 Quanto affermato si traduce nella seguente condizione (di ortogonalita):

RT2 (P2 +BT

2 Λ) = 0 (2.43)

Combinando tra di loro le equazioni (2.41) e (2.43), e queste con le prime

due delle (2.41), il sistema risolvente viene riscritto nel modo seguente:

5Si ricorda che lo spazio nullo, o nucleo, di una matrice A e definito come l’insieme di

quei vettori v tali che Av = 0.


U1 = K−11 (P1 +BT

1 Λ)

U2 = K+2 (P2 +BT

2 Λ) +R2α

[

(B1K−11 BT

1 +B2K+2 B

T2 ) B2R2

BT2 R

T2 0

][

Λ

α

]

=

=

[

−(B1K−11 P1 +B2K

+2 P2)

−RT2 P2

]

(2.44)

In presenza di N sottodomini, di cui Nf fluttuanti, il sistema risolvente

(2.38) viene riformulato in modo piu generale per tener conto delle eventuali

singolarita; con riferimento a quanto ottenuto per il caso particolare di due

sottodomini, risulta:

Us = K+s (Ps +BT

s Λ) +Rsαs s = 1, . . . , Nf

Us = K−1s (Ps +BT

s Λ) s = 1, . . . , (N −Nf )

[

Sd −G

−GT 0

][

Λ

α3

]

=

[

d

−e3

](2.45)

dove

Sd =N∑

s=1

(BsK+s B

Ts ) (2.46)

G = [−B1R1 · · · −BNfRNf

] (2.47)

d =N∑

s=1

(−BsK+s Ps) (2.48)

e3 = [PT1 R1 · · · PT

NfRNf

]T (2.49)

α3 = [αT

1 · · · αTNf

]T (2.50)


in cui la sottomatrice Sd, che prende il nome di complemento duale di Schur,

e simmetrica e definita positiva mentre i vettori Rs (con s = 1, . . . , Nf ) han-

no rango pari alla loro dimensione. Ne segue che il sottosistema di (2.45)

risulta simmetrico e non singolare e ammette un’unica soluzione [Λ,α3]T ;

questa, sostituita nelle prime di (2.45), consente di risolvere il problema. Si

osserva infine che, nel caso in cui non siano presenti domini fluttuanti, il

sistema (2.45) e equivalente al (2.38) poiche questo corrisponde al caso par-

ticolare di α3 = 0.

Risoluzione del problema di interfaccia

Il sistema (2.45) puo essere, di fatto, risolto direttamente, essendo il proble-

ma ben posto (numero di equazioni che uguaglia il numero delle incognite)

e date le caratteristiche del sistema di cui si e discusso in precedenza. I van-

taggi computazionali associati ai metodi di decomposizione emergono pero

in modo piu significativo associando ad essi dei risolutori di tipo iterativo;

in particolare, nel caso del FETI, si fa riferimento al metodo del gradiente

coniugato precondizionato (PCG, [10]) applicato al problema di interfac-

cia costituito dal sottosistema di (2.45). L’utilizzo di tale metodo iterativo

consente di non calcolare l’inversa della sottomatrice Sd.

Il PCG consiste in un algoritmo per la risoluzione di sistemi lineari

del tipo Ax = b sotto l’ipotesi che la matrice A risulti simmetrica e definita

positiva (ipotesi soddisfatta dal problema di interfaccia). La soluzione del

sistema puo essere determinata attraverso la minimizzazione di una funzione:

f(x) =1

2xTAx− xTb+ c

Per trovare il minimo si procede in modo iterativo partendo da un punto xi

e trovando il punto successivo attraverso la minimizzazione della funzione f

lungo una direzione scelta di uscente dallo stesso punto xi. La differenza del

PCG dal metodo del gradiente consiste proprio nella direzione di che viene

scelta A-ortogonale alle superfici di livello (cio costituisce di conseguenza un

vincolo). Sulla base di queste considerazioni, il problema di interfaccia e

equivalente al seguente problema di minimizzazione vincolata:

2.4 L’algoritmo di Gravouil e Combescure 27

min[

Φ(Λ) = 12Λ

TSdΛ− dTΛ]

GTΛ = e3

(2.51)

Le due equazioni di (2.51) costituiscono rispettivamente la funzione da mini-

mizzare ed il vincolo; quest’ultimo prende il nome di condizione di ammissi-

bilita. La variabile Λ (vettoriale in questo caso) viene espressa come somma

di un valore di inizializzazione e di un contributo che viene aggiornato ad

ogni iterazione:

Λ = Λ0 +PΛ∗ (2.52)

P e l’operatore di proiezione che garantisce il soddisfacimento della condizio-

ne di ammissibilita e che la ricerca di Λ∗ avvenga in un sottospazio prefissato

(corrispondente alla direzione di nell’esempio di partenza) ed e tale per cui:

GTP = 0 (2.53)

Si rimanda comunque all’appendice A per un maggiore approfondimento.

2.4 L’algoritmo di Gravouil e Combescure

Nel 2003 Gravouil e Combescure ([4]) formularono un nuovo algoritmo per

la risoluzione di problemi di dinamica strutturale; tale algoritmo e stato

recentemente riproposto in un’altra forma ([3]), senza pero cambiarne le

caratteristiche fondamentali.

Queste ultime possono essere riassunte nei seguenti punti:

- e ottenuto a partire dal metodo FETI e come tale appartiene alla

categoria dei metodi duali. In particolar modo, le incognite adottate

nel problema di interfaccia sono trazioni, che assumono allo stesso

tempo il ruolo di moltiplicatori di Lagrange;

- e di tipo multi-time-step, cioe consente di associare a ciascuno dei sot-

todomini in cui si immagina suddivisa la struttura un diverso schema di


integrazione, nonche un diverso passo temporale. Questo si traduce nel-

la possibilita di infittire la discretizzazione temporale dove necessario

e di adottare passi temporali piu grandi in altre zone della struttura;

- essendo stato ricavato direttamente dal metodo FETI, l’algoritmo GC

ne eredita tutti quegli aspetti positivi derivanti dalla decomposizione

dei domini, legati per lo piu all’abbattimento dei tempi di analisi e

quindi ad una riduzione dell’onere computazionale. Quanto affermato

e sicuramente valido per problemi meccanici; cio che ci si propone e

di dimostrare i vantaggi del metodo anche nell’ambito dei problemi

accoppiati, quale puo essere quello elettro-meccanico.

- un’ipotesi alla base dell’algoritmo originario e quella di interfaccia

perfetta.

Di seguito si riporta, a titolo di esempio, il caso particolare di suddi-

visione della struttura in due parti con diversi passi temporali, applicando

l’algoritmo al generico problema meccanico (che rappresenta il tipo di pro-

blema per il quale l’algoritmo e stato originariamente proposto). Il caso

generale di N sottodomini potra essere dedotto successivamente.

2.4.1 L’algoritmo GC applicato a due sottodomini

Si consideri il sistema di equazioni di equilibrio dinamico per un sistema a

k gradi di liberta scritto nell’istante tn+1:

Mun+1 +Kun+1 = f extn+1 (2.54)

dove M, K e f extn+1 rappresentano rispettivamente la matrice delle masse, la

matrice di rigidezza (entrambe simmetriche e definite positive) ed il vettore

dei carichi esterni all’istante tn+1. Tale sistema e ottenuto a valle della di-

scretizzazione spaziale del dominio di riferimento ed e costituito da equazioni

differenziali ordinarie nella variabile tempo: si tratta quindi di un sistema

semi-discretizzato.

Si suddivide l’intervallo di tempo in cui si vuole effettuare l’analisi in sot-

tointervalli ∆t; adottando lo schema di integrazione temporale di Newmark

(appendice B) e possibile esprimere i vettori accelerazione e spostamento

che compaiono nella (2.54) in funzione della velocita, secondo la formulazio-

ne proposta in [2] (si veda anche l’appendice B):


un+1 =pun +

β

γ∆tun+1 (2.55)

un+1 =pun +

1

γ∆tun+1 (2.56)

dove i predictor pun e pun sono definiti nel modo seguente:

pun = un +∆t(γ − β

γ

)

un +(γ − 2β

2γ

)

∆t2un (2.57)

pun =(γ − 1

γ

)

un − 1

γ∆tun (2.58)

Sostituendo le (2.55) e (2.56) nella (2.54) si ottiene il sistema algebrico in

cui l’unico vettore incognito e quello delle velocita:

Mun+1 = Fn+1 (2.59)

dove

M =1

γ∆tM+

β∆t

γK (2.60)

Fn+1 = f extn+1 −Mpun −Kpun (2.61)

Stesso passo temporale: ∆t1 = ∆t2

Suddividendo il dominio Ω in due sottodomini A e B non sovrapposti ai quali

e associato, in una prima trattazione, uno stesso passo temporale, il sistema

(2.59) puo essere scritto per ciascuno di essi mettendo contemporaneamente

in evidenza le forze di interfaccia, rappresentate dal vettore dei moltiplicatori

di Lagrange:

MAuAn+1 = FA

n+1 + LTAλn+1 (2.62)

MBuBn+1 = FB

n+1 + LTBλn+1 (2.63)

dove L rappresenta l’operatore di interfaccia che consente il passaggio dai

gradi di liberta del singolo sottodominio a quelli di interfaccia e viceversa.


Accanto ai precedenti sistemi e necessario considerare la condizione di con-

tinuita delle velocita all’interfaccia fra i due sottodomini:

LAuAn+1 + LBu

Bn+1 = 0 (2.64)

La strategia proposta in [2] per risolvere il sistema costituito dalle (2.62),

(2.63), (2.64) consiste nel decomporre la soluzione in una parte free (otte-

nuta risolvendo un problema non vincolato in cui intervengono i soli carichi

esterni) e in una parte link (soluzione di un problema vincolato in cui sono

applicate le forze di interfaccia):

uin+1 = ui

n+1/free + uin+1/link ∀i ∈ A,B (2.65)

In termini matriciali, il problema non vincolato puo essere espresso nel

seguente modo:

MA 0 0

0 MB 0

0 0 0

uAn+1/free

uBn+1/free

λn+1

=

FAn+1

FBn+1

0

(2.66)

Dalla struttura del sistema precedente si osserva come sia possibile risolvere

il problema non vincolato separatamente per ciascun sottodominio, poiche i

sottosistemi che compaiono non risultano accoppiati fra loro. Note le velocita

free e possibile risalire alle corrispondenti quantita cinematiche (spostamenti

e accelerazioni) attraverso le relazioni (2.55) e (2.56). D’altra parte il pro-

blema vincolato puo essere scritto come segue:

MA 0 −LTA

0 MB −LTB

−LA −LB 0

uAn+1/link

uBn+1/link

λn+1

=

0

0∑

i Liuin+1/free

(2.67)

in cui la terza del sistema (2.67) corrisponde alla (2.64), avendo suddiviso

fra contributo free e link.

Il sistema (2.67) consente, una volta nota la soluzione del problema non vin-

colato, di determinare i termini correttivi da aggiungere alle quantita free

ottenute dal sistema (2.66). E pero necessario determinare in via preliminare

le forze di interfaccia risolvendo il cosiddetto problema condensato, che


puo essere espresso nel modo seguente:

Hλn+1 = bn+1 (2.68)

dove H e bn+1 rappresentano rispettivamente l’operatore di condensazione

ed il termine noto non nullo nel sistema (2.67):

H =∑

i

Li(aiMi + biKi)−1LT

i (2.69)

bn+1 = −∑

i

Liuin+1/free (2.70)

dove ai e bi sono due coefficienti dipendenti dallo schema di integrazione

adottato. Si osserva che l’operatore di condensazione H e ottenuto dai con-

tributi forniti da ciascuno dei sottodomini che concorrono all’interfaccia;

l’aspetto forse piu interessante e pero il fatto che il calcolo dell’operatore

puo essere effettuato una sola volta all’interno di un codice di calcolo, non

essendo esso dipendente dall’istante temporale considerato.6

Diverso passo temporale: ∆t1 6= ∆t2

L’algoritmo GC ha fra i suoi punti di forza la possibilita di associare un

diverso passo temporale a ciascun sottodominio, sotto l’ipotesi che il passo

temporale maggiore sia proporzionale a quello minore attraverso un intero

m > 1:

∆T = m∆t

Inevitabilmente l’attribuzione di diversi passi temporali ai due sottodomini

comporta un’incompatibilita fra gli istanti temporali individuati dalle due

diverse discretizzazioni; e quindi necessario definire un opportuno operatore

di interpolazione che consenta di passare dalla macroscala alla microscala.

Le equazioni di equilibrio dei due sottosistemi A e B scritte rispettivamente

negli istanti tm (della macroscala) e tj (della microscala) sono espresse come:

6Quanto affermato vale nel caso di comportamento elastico dei sottodomini o, nel caso

in cui essi siano non-lineari, assumendo uno schema di integrazione di tipo esplicito.


Figura 2.9: Esempio di diverse discretizzazioni temporali ([2]).

MAuAm = FA

m + LTAλm

MBuBj = FB

j + LTBλj ∀j ∈ 1 . . .m

(2.71)

L’operatore di interpolazione e definito nel modo seguente:

uj =

(

1− j

m

)

u0 +j

mum (2.72)

Le operazioni che devono essere eseguite possono essere schematizzate nel

modo seguente:

a. Risoluzione del problema non vincolato sul sottodominio a cui e asso-

ciata la macroscala, cioe A:

MAuAm/free = FA

m (2.73)

Le altre quantita cinematiche free sono ricavate attraverso le relazioni

(2.55) e (2.56).

b. Ciclo j = 1 . . .m.

– Risoluzione del problema non vincolato sul sottodominio B (scala

piu fitta) all’istante tj :

MBuBj/free = FB

j (2.74)

– Interpolazione delle velocita free di A all’istante tj (passaggio

dalla macroscala alla microscala temporale):

uAj/free =

(

1− j

m

)

uA0/free +

j

muAm/free (2.75)


– Risoluzione del problema condensato all’interfaccia per la deter-

minazione dei moltiplicatori di Lagrange:

Hλj = −∑

i∈A,B

Liuij/free (2.76)

– Risoluzione del problema vincolato sul sottodominio B all’istante

tj :

MBuBj/link = LT

Bλj (2.77)

Si ricavano anche le altre quantita cinematiche link di B all’istante

tj con le consuete relazioni (2.55) e (2.56).

– Deduzione della soluzione totale di B all’istante tj :

uBj = uB

j/free + uBj/link

uBj = uB

j/free + uBj/link

uBj = uB

j/free + uBj/link

(2.78)

– Fine del ciclo quando j = m.

c. Risoluzione del problema vincolato per il sottodominio a all’istante tm:

MAuAm/link = LT

Aλm (2.79)

E ora possibile calcolare anche le altre quantita cinematiche link di A

all’istante tm con le (2.55) e (2.56).

d. Calcolo della soluzione totale di A all’istante tm:

uAm = uA

m/free + uAm/link

uAm = uA

m/free + uAm/link

uAm = uA

m/free + uAm/link

(2.80)

2.4.2 L’algoritmo GC applicato a s sottodomini con diversi

passi temporali

La sequenza di operazioni per il caso di due sottodomini puo essere estesa

anche al caso in cui gli stessi sottodomini siano in numero pari a s, mante-

nendo l’ipotesi di passi temporali multipli di quello della scala piu fitta. Si

riporta di seguito l’algoritmo cosı come proposto in [3]:


1. Inizializzazione dei tempi per ciascun sottodominio: tk = 0, k = 1, . . . , s, t =

0.

2. Ciclo sui passi temporali della scala piu fitta ∆t : t = t+∆t.

• Ciclo sugli s sottodomini.

- if tk < t ⇒ tk = tk +∆tk; risoluzione del problema free.

- if tk 6= t ⇒ interpolazione delle velocita free.

• Fine ciclo sui sottodomini.

• Risoluzione del problema condensato all’interfaccia.

• Ciclo sugli s sottodomini.

- if tk = t ⇒ risoluzione del problema link.

• Fine ciclo sui sottodomini.

3. Fine ciclo sui passi temporali della scala piu fitta.

Capitolo 3

Problema elettro-meccanico

accoppiato

3.1 Introduzione

Le moderne e complesse applicazioni scientifiche possono comportare l’inte-

razione fra diversi rami della fisica; si parla in questo caso di problemi a piu

campi o semplicemente multifisici. Un tipico esempio riguarda i cosiddetti

Micro-Sistemi Elettro-Meccanici (MEMS): si tratta di piccoli dispositivi per

i quali le forze gravitazionali risultano trascurabili nei confronti di quelle

elastiche, adesive e soprattutto elettrostatiche. A queste piccole scale, di

conseguenza, si presenta inevitabilmente un forte accoppiamento fra diversi

fenomeni fisici, in particolar modo fra gli spostamenti meccanici subiti dal

dispositivo e le forze elettrostatiche agenti sullo stesso: si tratta del cosid-

detto accoppiamento elettro-meccanico. Per comprendere le caratteristiche

fondamentali di questo tipo di problema si consideri il dispositivo di figura

3.1. Il sistema e costituito da un elettrodo deformabile di tipo mensola al

quale viene applicata una differenza di potenziale rispetto ad un secondo

elettrodo fisso. Tale attuazione provoca un immediato accumulo di carica

di segno opposto rispettivamente sulla superficie sottostante della mensola,

considerata come conduttore, e sull’altro elettrodo; le forze di Coulomb che

nascono per l’attrazione tra le cariche determinano la deformazione del di-

spositivo.

L’accoppiamento elettro-meccanico e caratterizzato da non-linearita: le ca-

riche elettriche si concentrano nelle zone del contorno dove la distanza tra

gli elettrodi e minore; di conseguenza in questi punti l’intensita delle forze e

35

3.2 Problema elettrostatico 36

Figura 3.1: Modello di accoppiamento elettro-meccanico ([22]).

maggiore e si genera una deformazione aggiuntiva dell’elettrodo. Il processo

deformativo termina al raggiungimento della situazione di equilibrio, quando

le forze elettrostatiche sono bilanciate dalle forze di richiamo elastico della

struttura.

Nei paragrafi che seguono verranno dedotte le equazioni che governano

i due problemi elettrico e meccanico, procedendo poi ad una loro integrazione

per ottenere le equazioni del problema elettro-meccanico accoppiato.

3.2 Problema elettrostatico

3.2.1 Equazioni di Maxwell ed equazione elettrostatica

Le relazioni fondamentali che governano il problema elettro-magnetico sono

costituite dalle leggi di Maxwell ; si tratta di equazioni differenziali alle deri-

vate parziali che descrivono i campi elettrico e magnetico e la loro interazione

([31]):

∇ ·D = ρ

∇ ·B = 0

∇×E = −∂B∂t

∇×H = i+ J+ ∂D∂t

(3.1)

dove E rappresenta il campo elettrico, B l’induzione magnetica, H il campo

magnetico, D il campo di spostamenti elettrici, ρ la densita di carica vo-

lumetrica, J la densita di carica superficiale dovuta alle cariche elettriche

libere, i e la densita di corrente imposta, ∇× e l’operatore rotore e ∇· l’o-peratore divergenza. Le precedenti si immaginano riferite ad un generico

dominio elettrico Ω. Per un materiale lineare, omogeneo e isotropo, le leggi

costitutive sono date dalle seguenti:


J = σE

B = µH

D = εE

(3.2)

dove la conduttivita σ, la permittivita ε e la permeabilita magnetica µ sono

parametri del materiale.

Se le variazioni che le quantita che compaiono nelle (3.1) subiscono

nel tempo sono trascurabili, il problema elettro-magnetico puo considerarsi

statico. Aggiungendo l’ipotesi di stazionarieta, le equazioni di Maxwell ri-

sultano disaccoppiate; in particolare la prima e la terza delle (3.1) governano

il problema elettrostatico:

∇ ·D = ρ

∇×E = 0(3.3)

La seconda delle (3.3) esprime la natura irrotazionale del campo elettrico; E

e conservativo e quindi puo essere espresso come gradiente di un potenziale

scalare φ:

E = −∇φ (3.4)

Tenendo conto invece della prima delle equazioni (3.3) e della (3.4), nonche

del legame costitutivo (3.2), si ottiene la seguente equazione differenziale nel

potenziale φ, altresı nota come equazione elettrostatica:1

∇2φ = −ρ

εin Ω (3.5)

dove ∇2(2) indica l’operatore di Laplace. Se non sono presenti cariche con-

centrate nel dominio elettrico, come nel caso in cui si tratti di vuoto o aria,

la densita di carica ρ e nulla per cui la (3.5) si riduce ad un’equazione di

Laplace:

∇2φ = 0 in Ω (3.6)

1La divergenza di un gradiente corrisponde al laplaciano.


Occorre osservare che le equazioni (3.1)-(3.6) sono valide quando Ω e compo-

sto da materiale dielettrico, caratterizzato da un valore finito di permittivita;

alle cariche e quindi impedito il movimento all’interno del dominio elettrico

e di conseguenza si accumulano sul suo bordo.

L’equazione differenziale (3.5) (o analogamente la (3.6)) deve essere corre-

data da condizioni al contorno. Esse sono essenzialmente di due tipi:

- potenziale imposto sulla porzione di contorno ΓV (Dirichlet):

φ = V su ΓV (3.7)

- densita superficiale di carica imposta sulla porzione di contorno ΓE

(Neumann):2

−ε∇φ · n = ρ su ΓE (3.8)

Come gia emerso nell’esempio di introduzione, nell’ambito dei MEMS

si ha a che fare con due tipi di materiali che differiscono per la mobilita delle

loro cariche elettriche: conduttori e isolanti (o dielettrici ideali). La mensola

di figura (3.1) si configura come conduttore: l’applicazione del potenziale tra

i due elettrodi di fatto carica la mensola ma le cariche non rimangono loca-

lizzate in una particolare regione della struttura; esse infatti sono libere di

muoversi da una parte all’altra del corpo, tanto che questo tipo di materiali

possono a tutti gli effetti essere caratterizzati da permittivita infinita. Il

movimento delle cariche avviene a velocita pressoche infinita, tale da annul-

lare una qualsiasi differenza di potenziale che si puo creare tra zone diverse

della mensola: ne segue quindi che al conduttore e associato un potenziale

costante. Il campo elettrico al suo interno e quindi nullo cosı come la sua

componente tangenziale sul contorno: le cariche infatti, potendo muover-

si liberamente, annullerebbero una qualsiasi differenza di potenziale anche

in direzione tangenziale. L’unica componente non nulla del campo elettri-

co e quindi quella perpendicolare allo stesso contorno. Dalle considerazioni

precedenti si deduce che il contorno rivolto verso la parte interna gap e da

considerarsi equipotenziale (condizione di tipo Dirichlet).

Si immagini poi che lo spazio compreso fra i due elettrodi di figura

(3.1) sia riempito con dell’aria, che rappresenta il tipico esempio di mate-

riale dielettrico. Questo tipo di materiali e caratterizzato da permittivita

2La condizione di Neumann corrisponde di fatto al teorema di Coulomb.


elettrica estremamente bassa se non addirittura nulla, e quindi da assenza

di cariche al loro interno. Sulla superficie del materiale isolante, dove non

sono presenti concentrazioni di carica, si impone, in accordo con la legge di

Gauss, che il campo elettrico in direzione normale sia nullo (condizione di

tipo Neumann).

Nel caso particolare in cui ci sia una porzione di frontiera che divide mate-

riali differenti, le due tipologie di condizioni al contorno vanno intese come

imposizione della continuita del potenziale elettrico e dell’uguaglianza della

densita superficiale di carica.

3.2.2 Formulazione ad elementi finiti del problema elettro-

statico

Definizione del modello per il potenziale

L’equazione (3.5) puo essere risolta in forma chiusa solo in alcuni semplici

casi; per tutti gli altri e necessario adottare un metodo approssimato quale

puo essere quello degli elementi finiti. Si introduca quindi una discretiz-

zazione spaziale del generico dominio elettrico Ω (figura 3.2) e si definisca

contemporaneamente sul singolo elemento finito un modello per il potenzia-

le, espresso in funzione dei potenziali nodali attraverso le funzioni di forma

Nφ:

φe =

nen∑

j=1

Nφjφej = NφΦ

e e = 1, . . . , nel (3.9)

con nel numero degli elementi della discretizzazione e nen numero dei nodi

dell’elemento e. In particolare si osserva che il potenziale φe sul singolo

elemento e una funzione scalare della posizione ed e espresso come prodot-

to fra i due vettori Nφ delle funzioni di forma e Φe dei potenziali nodali

dell’elemento. Si osserva che, contrariamente al caso meccanico in cui, a di-

scretizzazione introdotta, a ciascun nodo corrisponde un numero di incognite

(spostamenti nodali) pari alla dimensione del problema, nel caso elettrico ad

ogni nodo e associata un’unica incognita potenziale.

Analogamente al problema elasto-meccanico, si definisce il vettore glo-

bale delle incognite nodali Φ tale che:


Figura 3.2: Discretizzazione del dominio elettrico.

Φe = LeφΦ (3.10)

in cui le consuete matrici di connettivita Leφ consentono il passaggio dalla

numerazione locale a quella globale dei gradi di liberta e viceversa.

Si consideri ora il legame fra potenziale e campo elettrostatico espres-

so dalla (3.4). La componente i-esima del vettore E puo essere espressa nel

modo seguente:

Eei = − ∂φ

∂xi= −

nen∑

j=1

∂Nφj

∂xiφej (3.11)

In forma compatta:

Ee = −BφΦe (3.12)

in cui nella matrice Bφ sono contenute le derivate delle funzioni di forma.

Formulazione debole del problema elettrostatico

Si consideri l’equazione elettrostatica (3.5). Si moltiplichino entrambi i mem-

bri per una funzione test δφ che soddisfi alle condizioni al contorno essenziali:

in particolare tale funzione deve essere nulla in corrispondenza della porzio-

ne di contorno ΓV su cui e imposto il potenziale. Si integrino entrambi i

membri dell’equazione ottenuta sul dominio elettrico Ω:

∫

Ωε∇2φ · δφ dΩ = −

∫

Ωρ · δφ dΩ ∀ δφ = 0 su ΓV (3.13)


Applicando il teorema di Gauss (della divergenza) all’integrale del primo

membro, dopo alcuni passaggi algebrici, l’equazione diventa:

−∫

Ω(∇φ)T ε∇(δφ) dΩ = −

∫

Ωρ δφ dΩ+

∫

Γ(−ε∇φ · n) · δφ dΓ (3.14)

∀ δφ = 0 su ΓV

dove n e la normale alla frontiera di Ω, indicata con la lettera Γ. Introdu-

cendo le condizioni al contorno non omogenee (3.7) e (3.8) ed eliminando gli

integrali nulli sulla frontiera Γ, si ottiene la formulazione debole del proble-

ma:

∫

Ω(∇φ)T ε∇(δφ) dΩ =

∫

Ωρ δφ dΩ−

∫

ΓE

ρ δφ dΓ (3.15)

∀ δφ = 0 su ΓV

Sostituendo il modello per il potenziale (3.9) e per il campo elettrico (3.12)

nella (3.15) e considerando il contributo di tutti gli elementi della discretizza-

zione, si ottiene un sistema di equazioni algebriche che esprimono l’equilibrio

elettrico. In forma matriciale risulta:

KφΦ = Q (3.16)

dove:

- Kφ e la matrice di rigidezza elettrica associata al dominio Ω, simme-

trica e definita positiva:

Kφ =nel∑

e=1

LeT

φ

(∫

Ωe

BTφ εBφ dΩ

)

Leφ (3.17)

- Q e un vettore dipendente dalla densita di carica nel volume e sulla

superficie libera:

Q =nel∑

e=1

LeT

φ

(∫

Ωe

ρNT dΩ−∫

ΓEe

ρNT dΓ

)

(3.18)


in cui Ωe e ΓEe rappresentano rispettivamente il dominio e la porzione

di contorno ΓEe associati al singolo elemento finito e.

Risolto il sistema (3.16) e quindi possibile risalire al campo elettrico

con la seguente:

Ee = −Bφ LeφΦ (3.19)

Nel caso in cui il dominio elettrico sia costituito da vuoto o aria, non sono

presenti cariche concentrate al suo interno e le condizioni al contorno pos-

sono essere essenzialmente di due tipi: carica nulla o potenziale assegnato.

Analogamente al problema meccanico in cui si rende necessaria una distin-

zione fra gradi di liberta liberi e vincolati, si procede con una partizione del

vettore delle incognite Φ che distingue i nodi a voltaggio libero e a voltaggio

imposto. La partizione riguarda anche le altre quantita che compaiono nel

sistema risolvente (3.16); in particolare si indicano con Φ i gradi di liberta

liberi e con Φ quelli vincolati. La partizione consente di riscrivere il sistema

nel modo seguente:

[

Kφφ Kφφ

Kφφ Kφφ

][

Φ

Φ

]

=

[

0

q

]

(3.20)

Il sottovettore q, associato ai nodi con potenziale imposto, corrisponde alla

carica presente sulla superficie in comune col conduttore; quelle cariche sono

del tutto analoghe alle reazioni vincolari del problema meccanico. E questo

il caso di quella parte di contorno della mensola dell’esempio introduttivo

a contatto con l’aria (e rivolta all’elettrodo fisso), sulla quale e imposto il

potenziale V . Per calcolare i valori nodali di potenziale Φ, e sufficiente espli-

citare la prima riga del sistema matriciale (3.20), dalla quale si ricavano i

potenziali nei nodi liberi:

KφφΦ = −KφφΦ (3.21)

Sostituendo quindi i valori trovati nel secondo blocco di equazioni del sistema

(3.20) si ottengono i valori delle cariche-reazioni q:

q = KφφΦ+KφφΦ (3.22)


3.2.3 Esempio: problema elettrico 1D

In questo paragrafo viene presentato un modello schematizzato monodimen-

sionale, formato da due piastre parallele, alle quali viene applicata una dif-

ferenza di potenziale V (figura 3.3). Il gap tra le due piastre e riempito da

aria per cui il problema elettrostatico e governato dall’equazione di Laplace

(3.6), che nel caso monodimensionale assume la forma:

ε∂2φ

∂ξ2= 0 (3.23)

dove φ e il potenziale elettrico, ξ la coordinata spaziale e ε la permitti-

vita elettrica. Il dominio nel quale l’equazione (3.23) deve essere risolta e

ξ = [0, x] mentre le condizioni al contorno sono le seguenti:

φ = V per ξ = x

φ = 0 per ξ = 0(3.24)

La soluzione del problema puo essere ottenuta per integrazione diretta della

(3.23); ricavando i valori delle costanti di integrazione attraverso le (3.24),

si ricava:

φ =V

xξ (3.25)

L’espressione del campo elettrico risulta:

E = −∂φ

∂ξ= −V

x(3.26)

Infine le forze elettrostatiche, distribuite uniformemente sui piatti ed in di-

rezione perpendicolare rispetto alla superficie del conduttore, sono calcolate

come segue ([31]):

Fes =1

2εE2 =

1

2εV 2

x2(3.27)


Figura 3.3: Condensatore 1D con armature piane parallele ([22]).

3.2.4 Le forze elettrostatiche

Le forze elettrostatiche rivestono un ruolo fondamentale nell’analisi dei mi-

crosistemi proprio per le dimensioni che li caratterizzano; ne segue quindi

l’esigenza di una loro corretta valutazione. Sono essenzialmente due gli ap-

procci seguiti per il calcolo delle forze elettrostatiche: il metodo energetico e

il calcolo tramite il tensore degli sforzi di Maxwell.

Approccio energetico

Un primo modo per calcolare le forze elettrostatiche fa riferimento all’energia

elettrica totale del sistema (questo e costituito da un generico dominio elet-

trico Ω):

We =1

2

∫

Ωε∇φ∇φ dV (3.28)

Per valutare le forze elettrostatiche si utilizza il principio di conservazione

dell’energia: si scrive il lavoro virtuale compiuto dal sistema in seguito ad

una variazione δu del campo di spostamenti della struttura sul quale le forze

agiscono:

δL = fe · δu =∂We

∂uδu (3.29)

L’equazione (3.29) mette in evidenza come la forza fe sia legata alla variazio-

ne di energia elettrica dovuta ad un cambiamento virtuale di configurazione

(Ω ⇒ Ω′

). Dalla (3.29) si ricava quindi:


fe =∂We

∂u(3.30)

Di fatto il metodo piu immediato per il calcolo delle forze elettrostatiche e

il seguente.

Tensore degli sforzi di Maxwell

Un modo alternativo per determinare le forze elettrostatiche comporta il

calcolo del tensore degli sforzi di Maxwell e ne calcola il flusso rispetto alla

superficie considerata. Quel tensore e definito a partire dalla forza di Loren-

tz; essa, nel caso piu generale di problema elettro-magnetico, per una carica

(puntiforme) che si muove con velocita v, e data dalla seguente: 3

F = q(E+ v ×H) (3.31)

che nel caso di una distribuzione di carica nota ρ, puo essere scritta come

integrale di una forza per l’unita di volume:

F =

∫

f dV (3.32)

Alla base di questo approccio vi e l’idea di esprimere la forza f tramite le

derivate di un tensore, in modo tale che, conoscendo questo tensore sulla

superficie che racchiude un certo volume, sia possibile calcolare le forze che

agiscono sul volume stesso attraverso il teorema della divergenza o di Gauss.

In particolare se valesse la seguente:

fα = divTα =∑

β

∂

∂xβTαβ (3.33)

si avrebbe:

Fα =

∫

Vfα dV =

∫

VdivTα dV =

∫

Γ(V )Tα dΓ =

∫

Γ(V )

∑

β

Tαβ dΓβ (3.34)

dove Γ(V ) e la superficie che racchiude il volume V in cui il tensore Tαβ e

3Si suppone che la carica di prova q non alteri in modo apprezzabile i campi elettrico

E e magnetico H.


noto. Si osservi che Tαβ e rappresentato da una matrice 3 × 3 simmetrica

(caso tridimensionale, piu generale) ed e indicato come tensore degli sforzi.

Nel caso di dominio elettrico immerso nel vuoto, la forza di Lorentz

agente sulle cariche e dovuta solamente al campo elettrostatico:

Fα =

∫

Vfα dV =

∫

VρEα dV = ε

∫

V(∇ ·E)Eα dV (3.35)

Nell’ultimo passaggio della (3.35) si e utilizzata la prima delle leggi di Max-

well (3.3). Dalla (3.35) si ricava:

fα = ε(∇ ·E)Eα

che puo essere rimaneggiata per ottenere il tensore T; infatti, sfruttando le

proprieta delle derivate di un prodotto di funzioni, si ottiene:

fα = εEα(∇ ·E) = εEα∂

∂xβEβ = ε

[∂

∂xβ(EαEβ)− Eβ

∂

∂xβEα

]

(3.36)

Un ulteriore passo avanti nell’equazione (3.36) puo essere compiuto se si con-

sidera la relazione seguente:

∂Eα

∂xβ=

∂Eβ

∂xα(3.37)

La (3.37), che risulta un’identita se α = β, e valida anche nel caso generale

di α 6= β, come conseguenza delle proprieta irrotazionali del campo elettro-

statico (3.3).

Il secondo termine tra parentesi della (3.36) viene quindi modificato e con

esso l’intera equazione:

fα = ε

[∂

∂xβ(EαEβ)− Eβ

∂

∂xβEβ

]

=

=∂

∂xβε

[

EαEβ − 1

2δαβ|E| · |E|

]

=∂

∂xβTαβ (3.38)

Si giunge cosı alla definizione del tensore degli sforzi nel caso tridimensionale:

3.3 Problema elettro-meccanico 47

T = ε

E2x − 1

2 |E|2 ExEy ExEz

EyEx E2y − 1

2 |E|2 EyEz

EzEx EzEy E2z − 1

2 |E|2

(3.39)

dove Ei = − ∂φ∂xi

e il campo elettrico in direzione xi. Analogamente allo sfor-

zo meccanico, la forza elettrostatica nel dominio Ω e esprimibile come ∇·T;

quindi la forza esercitata da Ω sui domini circostanti e:

fe = −∫

Ω∇ ·T dV (3.40)

L’integrale (3.40) risulta comunque scomodo per la determinazione delle for-

ze elettrostatiche in quanto comporta il calcolo delle derivate spaziali seconde

del potenziale elettrico φ. Per superare tale difficolta, come accennato all’i-

nizio del paragrafo, viene utilizzato il teorema di Gauss, in modo da ricavare

le forze elettriche come integrale sul contorno Γ.

fe = −∫

Ω∇ ·T dV =

∫

ΓT · n dS (3.41)

dove n e la normale interna al bordo del dominio elettrico Ω e allo stesso

tempo normale esterna per i domini adiacenti.

Infatti, con riferimento all’esempio di inizio capitolo, sul contorno della

mensola il potenziale assume valore costante e l’unica componente del campo

elettrico diversa da zero e quella perpendicolare al bordo stesso; l’equazione

(3.41) si riduce alla scrittura seguente:

fe =1

2

∫

Γε |∇φ|2 n dS =

1

2

∫

Γε |E|2 n dS (3.42)

3.3 Problema elettro-meccanico

3.3.1 Considerazioni iniziali

Nel paragrafo precedente sono state dedotte le equazioni che governano il

problema elettrico. Si e anche mostrato che il dominio Ω in cui quel problema

viene risolto e costituito in genere da materiale dielettrico, caratterizzato


dall’assenza di cariche interne, e puo confinare con dei conduttori, nei quali

il potenziale risulta costante. Per questo motivo, come gia messo in evidenza,

su quella parte di contorno del dominio elettrico in comune con il conduttore

viene assegnata una condizione di tipo Dirichlet (potenziale imposto).

Il principio che sta alla base dell’accoppiamento elettro-meccanico e

gia stato esposto nell’esempio di introduzione al capitolo: una differenza

di potenziale imposta dall’esterno provoca un accumulo di carica di segno

opposto sui bordi dei due elettrodi e inevitabilmente nascono delle forze

attrattive che agiscono su di essi; deformandosi gli elettrodi provocano un

cambiamento nella configurazione del dominio elettrico (costituito dal gap

compreso fra di essi) e quindi nella distribuzione di potenziale dal quale,

come e noto, dipendono a loro volta le stesse forze elettrostatiche. Questo

fenomeno e caratterizzato da non-linearita, poiche l’accumulo di cariche in

quelle zone degli elettrodi per i quali il gap e minore provoca a sua volta un

aumento dell’intensita delle forze di Coulomb in quei punti.

In definitiva, l’accoppiamento fra i due problemi meccanico ed elettrico

si esplica nei due seguenti modi:

- dal punto di vista del problema meccanico, la parte elettrica influisce in

termini di forze elettriche agenti su quella parte di contorno in comune

con il dominio elettrico; quelle forze dipendono dalla distribuzione di

potenziale;

- d’altra parte, gli spostamenti del dominio meccanico provocano una

deformazione non solo dell’interfaccia ma anche dell’intera parte elet-

trica; questo cambiamento di configurazione provoca a sua volta una

modifica nell’andamento del potenziale e quindi influisce sulla soluzione

del problema elettrico.

Dalle precedenti osservazioni emerge una dipendenza reciproca dei due pro-

blemi meccanico ed elettrico, tanto che non e possibile ricavare la soluzione

di un problema senza aver ricavato anche quella dell’altro. Data poi la non-

linearita, risulta inevitabile l’utilizzo di una procedura di tipo iterativo per

risolvere il problema elettro-meccanico.


Figura 3.4: Problema elettro-meccanico di riferimento.

3.3.2 Formulazione del problema elettro-meccanico accop-

piato

In questo paragrafo si vuole procedere con una sintesi delle equazioni che

governano i due problemi, sulla base delle precedenti osservazioni. Si consi-

derino due domini meccanico ed elettrico Ωm e Ωe (figura 3.4). Assumendo

che non ci siano cariche libere in Ωe, l’equazione del problema elettrostatico

e la seguente (di seguito si adotta la notazione di Einstein ad indici ripetuti):

∂

∂xi

(

ε∂φ

∂xi

)

= 0 (3.43)

D’altra parte, l’equazione di equilibrio dinamico del problema meccanico

(costituito dal generico solido Ωm), nell’ipotesi di piccoli spostamenti e de-

formazioni, e data da:

−ρmuj +∂

∂xiσji + Fj = 0 (3.44)

dove ρm e la densita del materiale solido, σ e il tensore degli sforzi meccanici e

Fj rappresenta le forze esterne. Per semplificare la trattazione del problema,

si assume che le uniche forze esterne agenti siano quelle di Coulomb.


Per calcolare il campo di spostamenti dovuto ai carichi applicati e ne-

cessario introdurre una legge costitutiva del materiale, che leghi gli sforzi σij

alle deformazioni εij . Si esprime innanzitutto εij come funzione del campo

di spostamenti (congruenza):

εij =1

2

(∂uj∂xi

+∂ui∂xj

)

(3.45)

σij = Cijhkεhk (3.46)

dove Cijhk e il tensore di Hooke, composto da termini costanti per l’ipotesi

di isotropia del materiale. Si puo osservare che il problema meccanico, con-

siderato indipendente da quello elettrico, ha carattere lineare.

La distribuzione delle forze elettrostatiche segue la legge di Coulomb:

Fj = qEj (3.47)

dove q e la densita di carica sulla superficie ed Ej e il campo elettrico gene-

rato dalle cariche poste sull’altro elettrodo. In prossimita delle cariche sulla

superficie del conduttore, meta del campo elettrico e generato dalle stesse

cariche, che non danno contributo alla forza elettrostatica sulla superficie

dove sono concentrate; quindi, se Ej e il campo elettrico totale nel dominio,

Fj = q(Ej/2).

In accordo con la legge di Gauss, la carica q sulla superficie di contatto

fra le due parti puo essere espressa in funzione del campo elettrico al bordo:

q = εEini (3.48)

dove ni e il versore normale all’interfaccia rivolta verso il dominio elettrico.

Ricordando l’equazione (3.4), la forza elettrostatica puo essere espressa come

funzione del potenziale elettrico:

Fj =1

2ε

(∂φ

∂xi

∂φ

∂xi

)

nj (3.49)

Anche utilizzando il tensore di Maxwell, emerge la dipendenza delle forze

elettrostatiche dal quadrato del campo elettrico (equazione (3.42)). Quindi,

se da una parte i due problemi meccanico ed elettrico, analizzati in modo


indipendente, presentano un comportamento lineare, dall’altra la risposta

del problema elettro-meccanico e complessivamente non-lineare. Come di-

mostra l’equazione (3.49), questo e dovuto alla dipendenza quadratica della

forza elettrostatica dal potenziale; a questo si aggiunge la dipendenza della

forma del dominio elettrico dagli spostamenti di quello meccanico.

3.3.3 Dominio meccanico

Si consideri il consueto solido deformabile Ωm adiacente ad un dominio elet-

trico Ωe, con il quale confina attraverso l’interfaccia Γ. Si introduca una

discretizzazione spaziale analogamente a quanto fatto nella deduzione del

sistema delle equazioni di equilibrio elettrico. Definendo un opportuno mo-

dello per lo spostamento sul singolo elemento della mesh ed introducendo

tale modello nell’espressione della variazione prima del funzionale energia

potenziale totale ([2]), si giunge alla scrittura del sistema di equilibrio dina-

mico semi-discretizzato, costituito da equazioni differenziali ordinarie nella

variabile tempo:

Muuu+Kuuu = f (3.50)

dove Muu e la matrice delle masse e Kuu e la matrice di rigidezza meccanica,

entrambe simmetriche e definite positive nel caso di solido ben vincolato. 4

3.3.4 Deformazione del dominio elettrico: problema pseudo-

meccanico

Le (3.43) e (3.44) rappresentano le equazioni differenziali alle derivate par-

ziali che governano rispettivamente il problema elettrico e quello meccanico.

Una volta discretizzate, esse danno luogo ai due sistemi elettrico (3.16) e

meccanico (3.50) di cui gia si e discusso. Questi sistemi devono ora esse-

re accoppiati, sulla base delle considerazioni del paragrafo 3.3.1, per dar

luogo alle equazioni del problema elettro-meccanico; infatti i due sistemi

4I pedici uu sono una prerogativa di questo capitolo e consentono di identificare in

maniera piu immediata, al prezzo di una maggior pesantezza della trattazione, le quantita

in gioco ed in particolare di comprendere se esse si riferiscono alla parte meccanica oppure

a quella elettrica.


(3.16) e (3.50) considerati singolarmente governano i due problemi elettrico

e meccanico disaccoppiati.

Sulla base delle considerazioni iniziali del paragrafo 3.3.1 il sistema

governante il problema elettro-meccanico e espresso dalla seguente:

Muuum = f elec(ue, φ, V ) su Ωm

Kφ(ue)Φ = Q(ue, V ) su Ωe

(3.51)

Il sistema (3.51) richiede alcune precisazioni:

- il sottosistema meccanico e stato discretizzato anche temporalmente

adottando uno schema di integrazione, quale puo essere quello di New-

mark (appendice B); in questo modo si e ottenuto un sistema algebrico

nel quale la matrice M e espressa in funzione delle matrici delle masse

Muu e di rigidezza Kuu e dei parametri di integrazione;

- si e ipotizzato che il sistema meccanico sia soggetto esclusivamente a

forzanti elettriche tralasciando quelle gravitazionali, ipotesi realistica

nel campo dei micro-sistemi. Quelle forze dipendono dalla configu-

razione del dominio elettrico (attraverso gli spostamenti nodali ue),

dalla distribuzione di potenziale sullo stesso dominio e dal potenziale

imposto V ;

- per quanto riguarda il sistema elettrico, e stata evidenziata la dipen-

denza della corrispondente matrice di rigidezza dalla configurazione del

dominio elettrico stesso attraverso gli spostamenti nodali ue. E stata

inoltre messa in evidenza la dipendenza del vettore dei termini noti

dagli spostamenti ue e dal potenziale imposto.

Il calcolo delle forze elettrostatiche che compaiono nel primo sottosi-

stema di (3.51) e gia stato discusso nel paragrafo 3.2.4; si vuole mostrare ora

un metodo per tener conto degli spostamenti nodali del dominio elettrico.

Si consideri ancora il generico problema elettro-meccanico di figura 3.4: gli

spostamenti um della parte meccanica inducono, attraverso l’interfaccia Γ,

una deformazione della parte elettrica; sono essenzialmente due i metodi che

consentono di tener conto di tale cambiamento di configurazione: re-meshing

e mesh deformation. Il primo comporta la generazione di una nuova mesh

elettrica ogni volta che si modifica la configurazione di Ωe per effetto di um;

il secondo aggiorna la mesh esistente attraverso il calcolo degli spostamenti

dei nodi della mesh elettrica stessa.


Nel seguito si fara riferimento a questo secondo metodo, che prevede

quindi il calcolo degli spostamenti nodali ue della parte elettrica in funzione

di quelli um della parte meccanica. In realta, tale dipendenza si riduce ai

soli spostamenti dell’interfaccia Γ, indicati con ui, che possono essere estratti

dal vettore um degli spostamenti dell’intero dominio meccanico Ωm attra-

verso una matrice di connettivita opportunamente definita (infatti ui ⊆ um):

ui = Lium (3.52)

In generale, la relazione che lega gli spostamenti di interfaccia con quelli

elettrici puo essere espressa come segue:

g(ue,ui) = 0 (3.53)

La scelta del tipo di legame (3.53) e varia; si preferisce, per semplicita, risol-

vere un problema elasto-meccanico sul dominio elettrico (indicato appunto

con problema pseudo-meccanico) nel quale siano assenti carichi esterni e sia-

no solamente imposti i cedimenti ui in corrispondenza del bordo coincidente

con Γ; in termini matematici si tratta di risolvere il sistema discretizzato:

Kelecue = Fps su Ωe (3.54)

in cui ue sono gli spostamenti nodali incogniti, Kelec e la matrice di rigidez-

za del dominio pseudo-meccanico (geometricamente coincidente con quello

elettrico); Fps, infine, rappresenta il vettore dei carichi nodali equivalenti. Il

sistema (3.54) e corredato dalla seguente condizione al contorno:

ue = ue = ui su Γ

Il vettore Fps e un vettore di termini nulli poiche non agiscono carichi sul

dominio pseudo-meccanico. Si esegue una partizione delle grandezze presen-

ti in (3.54) secondo i gradi di liberta liberi (ue) e vincolati (ue), mettendo

in evidenza le reazioni vincolari incognite R in corrispondenza dei gradi di

liberta vincolati:


[

Kelec(ueue)Kelec(ueue)

Kelec(ueue)Kelec(ueue)

][

ue

ue

]

=

[

0

R

]

(3.55)

Dalla prima delle (3.55) e possibile esplicitare la relazione (3.53) nel caso di

problema pseudo-meccanico:

Kelec(ueue)ue = −Kelec(ueue)

ui (3.56)

In realta gli spostamenti ui dell’interfaccia Γ che compaiono nella (3.54)

possono essere ricondotti al vettore globale degli spostamenti meccanici um

attraverso la (3.52); per questo la (3.56) viene aggiornata come segue:


Lium (3.57)

Risolvendo si ottiene:

ue = −K−1elec(ueue)

Kelec(ueue)Lium (3.58)

La (3.58) esplicita la dipendenza della configurazione del dominio elettri-

co (rappresentata dal vettore di spostamenti nodali ue) dagli spostamenti

meccanici um.

La (3.57) deve essere aggiunta alle equazioni iniziali (3.51) a formare

il seguente sistema:

Muuum = f elec(ue, φ, V ) su Ωm


Lium su Ωe

Kφ(ue)Φ = Q(ue, V ) su Ωe

(3.59)

Combinando le equazioni in (3.59), condensando in particolar modo la se-

conda, si ottiene il sistema accoppiato:

Muuum = f elec(ue(um), φ, V ) su Ωm

Kφ(ue(um))Φ = Q(ue(um), V ) su Ωe

(3.60)

Il (3.60) puo essere infine scritto in maniera compatta:


Figura 3.5: Spostamenti meccanici, elettrici e di interfaccia.

Muuum = f elec(um, φ, V ) su Ωm

Kφ(um)Φ = Q(um, V ) su Ωe

(3.61)

Si e in questo modo ottenuta la diretta dipendenza delle quantita elettriche

che compaiono nel (3.61) dagli spostamenti meccanici um, esprimendo ma-

tematicamente l’effetto della parte meccanica su quella elettrica.

3.3.5 Esempio: problema elettro-meccanico 1D

Si riporta ora un esempio monodimensionale che vuole mettere in evidenza

la caratteristica di non-linearita del problema elettro-meccanico accoppia-

to, in particolare la relazione tra le forze elettrostatiche e gli spostamenti

meccanici.

La figura 3.6 mostra un condensatore a facce piane parallele immerso

nel vuoto. Le due piastre sono corpi rigidi che si trovano inizialmente ad una

distanza x; in seguito all’applicazione di una differenza di potenziale V , le

cariche elettriche si distribuiscono in modo uniforme su di esse. La risposta

elastica e affidata ad una molla di costante k, collegata alla massa meccanica

m.

Il problema elettrostatico e governato dall’equazione di Laplace (3.6),

che nel caso monodimensionale assume la forma:


Figura 3.6: Micro-risonatore monodimensionale ([22]).

ε∂2φ

∂ξ2= 0 (3.62)

dove φ e il potenziale elettrico, ξ e la coordinata spaziale e ε0 la permit-

tivita elettrica. Il dominio nel quale l’equazione (3.62) deve essere risolta

e ξ = [0, x], dove x descrive la posizione dell’elettrodo mobile che nel caso

dinamico puo essere espresso come x(t). Le condizioni al contorno per il caso

in esame sono:

φ = V per ξ = x

φ = 0 per ξ = 0(3.63)

Per quanto riguarda il problema elettrico, i risultati sono del tutto analoghi

a quelli ricavati nell’esempio 1D del paragrafo (3.2.3). L’andamento del po-

tenziale:

φ =V

xξ (3.64)

Il campo elettrostatico:

E = −∂φ

∂ξ= −V

x(3.65)

3.4 Instabilita elettro-meccanica: fenomeno dipull-in 57

Infine, le forze elettrostatiche agenti sulla massa m:

Fes =1

2εE2 =

1

2εV 2

x2(3.66)

La forza meccanica applicata all’elettrodo a causa della presenza della molla

vale:

Fmech = k(x− x0) (3.67)

dove k e la costante di rigidezza della molla e x0 e la posizione iniziale di

riferimento della stessa.

Si scrive quindi l’equazione di equilibrio dinamico del problema mono-

dimensionale:

mx+ k(x− x0) +1

2εV 2

x2= 0 (3.68)

Nonostante sia valida per il semplice caso monodimensionale, l’equazione di

equilibrio (3.68) esprime la non-linearita del problema elettro-meccanico, che

risiede, in questo caso, nella presenza di x2 a denominatore del termine della

forza elettrica.

3.4 Instabilita elettro-meccanica: fenomeno di

pull-in

3.4.1 Caso statico

Per comprendere il significato fisico del pull-in, si consideri il problema mo-

nodimensionale descritto nel paragrafo 3.3.5 e schematizzato in figura 3.6.

L’equilibrio statico e governato dalla relazione:

k(x− x0)︸︷︷︸

forza elastica

+1

2εV 2

x2︸︷︷︸

forza elettrostatica

= 0 (3.69)

che fornisce l’espressione V = V (x) del voltaggio imposto in funzione della


distanza tra le piastre; la (3.69) mostra come la forza elastica sia propor-

zionale al valore di spostamento (u = x − x0), mentre la forza elettrostati-

ca inversamente proporzionale al quadrato della distanza tra le piastre del

condensatore (x). Come dimostrato sperimentalmente per dispositivi reali,

anche nel modello qui adottato gli elettrodi si avvicinano progressivamente

al crescere del valore di voltaggio imposto fino a quando la forza di richiamo

elastico della molla non e piu in grado di contrastare la forza elettrostati-

ca; l’equilibrio non e piu soddisfatto e l’armatura mobile del condensatore

collassa sull’elettrodo fisso, col verificarsi di una situazione di cortocircuito.

Questo fenomeno prende il nome di pull-in.

Risulta fondamentale, soprattutto nella progettazione dei micro-sistemi,

conoscere la coppia di valori (VPI , xPI) corrispondenti alla situazione limite,

che conduce il sistema all’instabilita. Analiticamente la ricerca di quei valo-

ri puo avvenire imponendo la stazionarieta della curva V (x) ottenuta dalla

(3.69):

dV (x)

dx= k − ε0

V 2PI

x3PI

= 0 ⇒ VPI =

√

x3PI · kε0

(3.70)

Il valore ottenuto VPI rappresenta il massimo voltaggio imponibile. La di-

stanza fra le piastre del condensatore (gap) corrispondente al pull-in si ot-

tiene sostituendo VPI nella (3.69):

k(xPI − x0) +1

2ε0

V 2PI

x2PI

= 0 ⇒ xPI =2

3x0 (3.71)

Quindi l’espressione esatta di VPI risulta:

VPI =

√

8x30 · k27ε0

(3.72)

Si osservi che il valore di pull-in e una caratteristica intrinseca del problema,

in quanto dipende da grandezze invarianti: rigidezza e gap iniziale. Il grafico

di figura 3.7 mostra l’andamento della funzione (3.69) normalizzata rispetto

al valore VPI :

(V

VPI

)2

=27

4

(x0 − x

x0

)[(x0 − x

x0

)

− 1

]2

(3.73)


Figura 3.7: Andamento della curva spostamento-voltaggio dell’esempio 1D

([22]).

La (3.73) e un’equazione algebrica di secondo grado, che rispecchia la non-

linearita del problema elettro-meccanico; ad un qualsiasi valore di V corri-

spondono in genere due soluzioni in termini di spostamento. Dal grafico 3.7

si osserva che tra le due situazioni di equilibrio corrispondenti al generico V

quella con spostamento maggiore (e quindi con elettrodo deformabile pros-

simo a quello fisso superiore) e instabile; il ramo instabile e rappresentato

dalla linea tratteggiata. Infatti, imprimendo una piccola perturbazione al-

l’elettrodo mobile, questo collassa sulla piastra superiore. Viceversa, i punti

corrispondenti a spostamenti minori rappresentano situazioni di equilibrio

stabile (linea continua). La coppia di valori di pull-in (VPI , xPI) costituisce

lo spartiacque tra i due rami della curva.

Nonostante i risultati ottenuti siano validi per il caso monodimensio-

nale analizzato, e comunque interessante notare che anche nei microsistemi

reali il gap di pull-in e prossimo ai due terzi della distanza iniziale tra le

armature ([22]).

Si sottolinea ancora una volta l’importanza della corretta valutazio-

ne del voltaggio e spostamento di pull-in nell’ambito della progettazione

dei microsistemi elettro-meccanici (MEMS); l’instabilita elettro-meccanica,

che comporta il collasso dell’elettrodo mobile su quello fisso, e infatti una

situazione che deve essere scongiurata in quanto, una volta verificatasi, il

dispositivo risulta irrimediabilmente compromesso.


Figura 3.8: In alto la curva statica di pull-in, in basso la risposta oscillatoria

per un passo di voltaggio costante VD ([22]).

3.4.2 Caso dinamico

Si consideri il problema monodimensionale di figura 3.6 e se ne scriva l’equa-

zione di equilibrio dinamico:

mx+ k(x− x0) +1

2ε0

V 2

x2= 0 (3.74)

La (3.74) e un’equazione differenziale ordinaria del secondo ordine nella va-

riabile tempo; note le condizioni iniziali e possibile, attraverso l’uso di me-

todi di integrazione numerica, calcolarne una soluzione approssimata. Un

secondo modo di procedere ([22]) mantiene fissa l’attuazione V e lineariz-

za la soluzione intorno ad una configuarazione di equilibrio di riferimento:

x = xeq + ∆x, dove xeq rappresenta la soluzione statica. Sostituendo la

linearizzazione nella (3.74), si ottiene:

mxeq +m∆x+ k(xeq +∆x− x0) +

1

2ε0

V 2

(xeq +∆x)2= 0 (3.75)

dove il termine mxeq e nullo, essendo xeq soluzione statica e quindi costan-


te. Esprimendo il termine elettrostatico attraverso uno svilupppo in serie di

Taylor, l’equazione diventa:

m∆x+ k∆x+ k(xeq − x0) +1

2ε0

V 2

(xeq)2− ε0

V 2

(xeq)3∆x+ t.o.s. = 0 (3.76)

Considerando che xeq e soluzione della (3.69) e trascurando i termini di or-

dine superiore, l’equazione di equilibrio dinamico si riduce a:

m∆x+

(

k − ε0V 2

x3eq

)

∆x = 0 (3.77)

La (3.77) e un’equazione differenziale ordinaria del second’ordine. Per il pro-

blema considerato, la frequenza propria e:

ω =

√

k

m− ε0

V 2

mx3eq(3.78)

La (3.78) mostra la dipendenza della frequenza propria del sistema ω dal

voltaggio imposto V. Ponendo (V, xeq) = (VPI , xPI) la frequenza propria si

annulla. Se V < VPI e xeq e nella zona stabile della curva di figura (3.7),

ω e reale; viceversa, se xeq si trova nella parte instabile della curva, ω e un

numero complesso.

Si analizza ora la risposta transitoria del modello sottoposto a vol-

taggio crescente. Per valori di potenziale compresi tra 0 e un generico VD

l’equazione del problema in esame e la seguente (insieme alle condizioni ini-

ziali):

mx = −k(x− x0)− 12ε0

V 2D

x2 = 0 t > 0

x = 0 t = 0

x = 0 t = 0

x = 0 t = 0

(3.79)

Non e possibile calcolare la soluzione del problema (3.79) in forma chiusa;

attraverso l’utilizzo di un algoritmo di integrazione numerica nel tempo,

come il metodo di Runge-Kutta del quarto ordine ([22]), e possibile ottenere

la risposta del sistema.


Se VD e relativamente basso, la risposta oscillatoria e periodica; l’am-

piezza della soluzione dipende dal valore del voltaggio applicato, dalla ri-

gidezza meccanica, dalla permittivita elettrica e dalla distanza iniziale tra

le armature, ma non dipende dalla massa dell’elettrodo mobile. Al contra-

rio il periodo T e funzione della massa m. Nella parte inferiore di figura

3.8 si osserva l’andamento della risposta transitoria per un voltaggio pari

a VD = 0.9VPI , mentre la parte superiore mostra la soluzione di equilibrio

statico (cerchio rosso) per lo stesso valore di potenziale. Le frecce poste al

centro indicano la direzione delle forze elettro-meccaniche; essa cambia in

corrispondenza della soluzione di equilibrio statico.

La risposta dinamica oscilla asimmetricamente intorno alla soluzione di equi-

librio stabile. In dinamica e possibile ottenere uno spostamento maggiore

rispetto a quello di pull-in statico (triangolo verde in figura 3.9).

Durante il transitorio del sistema, lo spostamento puo raggiungere valori

maggiori di quello che si avrebbe nella situazione di equilibrio instabile; a

questo punto le forza elettrostatica cambia direzione e l’elettrodo mobile col-

lassa sull’armatura fissa. Per valori di voltaggio inferiori a 0.9VPI questo non

puo accadere; per l’esempio analizzato lo spostamento di equilibrio istabile

viene raggiunto per valori VD = 0.92VPI . La Figura 3.8 mostra il transi-

torio per questo valore di voltaggio, chiamato pull-in dinamico (VDPI). In

precedenza si e osservato come sia complessa la ricerca della soluzione della

(3.79) in corrispondenza di un passo a voltaggio imposto; risulta invece piu

semplice il calcolo del valore esatto di VDPI . In corrispondenza del pull-in

dinamico, l’accelerazione e nulla, trattandosi di un punto di equilibrio del

sistema:

k(x− x0)2 +

1

2ε0

V 2DPI

x2DPI

= 0 (3.80)

In accordo con il principio di conservazione dell’energia, durante il transito-

rio l’energia stessa deve mantenersi costante; in corrispondenza di (VDPI , x0)

e (VDPI , xDPI) l’energia e la stessa:

1

2k(xDPI − x0)

2 − 1

2ε0

V 2DPI

xDPI= −1

2ε0

V 2DPI

x0. (3.81)

Nelle equazioni (3.80) e (3.81) le uniche due quantita incognite sono VDPI e

xDPI . Risolvendo si ottiene:


Figura 3.9: Pull-in dinamico nel caso 1D: in alto la curva di pull-in statico,

in basso la risposta transitoria per V = VDPI ([22]).

VDPI =

√

kx304ε0

=3

8

√6VPI ≈ 0.92VPI (3.82)

xDPI =1

2x0. (3.83)

Si osserva come il valore di spostamento di pull-in dinamico, analogamente

al caso statico, dipenda unicamente dalla distanza iniziale delle armature del

condensatore.

Capitolo 4

Metodi di risoluzione

4.1 Introduzione

Le relazioni che governano i problemi elettrico e meccanico (equazioni 3.43

e 3.44) descritte nel paragrafo 3.3 sono equazioni differenziali alle derivate

parziali. La soluzione puo essere ottenuta in forma chiusa solo per alcuni

casi particolari; in generale e necessario introdurre una discretizzazione dei

domini su cui quelle stesse equazioni sono scritte, consentendo il passaggio

da un numero teoricamente infinito di gradi di liberta (continuo) ad un

numero finito (discreto). Esistono diversi metodi di discretizzazione spaziale:

il piu diffuso per domini meccanici e senza dubbio quello degli elementi finiti

(FEM); per la parte elettrica sono invece egualmente utilizzati sia il FEM

sia il metodo degli elementi di contorno (BEM), con i quali viene modellata

l’interfaccia fra le due parti meccanica ed elettrica. Nelle analisi che verranno

svolte si considerera il primo metodo anche per il dominio elettrico, come

peraltro gia introdotto nel paragrafo 3.2.2.

Il problema elettro-meccanico accoppiato e governato dalle seguenti

equazioni, dedotte al capitolo 2. In particolare, per seguire la trattazione di

[22], si escludono effetti dinamici sulla parte meccanica:

Kmechuu u = f elec(u,Φ, V )

Kφφ(u)Φ = q(u, V )(4.1)

dove Kφφ e la matrice di rigidezza elettrica, Kmechuu la matrice di rigidezza

meccanica, Φ e il vettore dei potenziali nodali, u il vettore degli spostamenti

meccanici nodali, q il vettore delle cariche nodali, V e il potenziale imposto

64

4.1 Introduzione 65

Figura 4.1: Approccio co-simulation (a sinistra) e fully coupled-simulation

(a destra).

sul bordo con condizioni di Dirichlet e f elec sono le forze elettrostatiche ap-

plicate alla parte meccanica.1

Nelle equazioni (4.1) e evidente la gia citata non-linearita del problema

elettro-meccanico accoppiato: innanzitutto le forze elettrostatiche applicate

alla struttura dipendono non-linearmente dal potenziale imposto V (equa-

zione (3.49)) e dalla configurazione del dominio elettrico; in secondo luogo,

sussite la dipendenza del dominio elettrostatico dal campo di spostamenti

meccanici.

La soluzione del sistema (4.1) puo essere cercata seguendo due diversi

tipi di approcci:

- partizionati, o co-simulations;

- monolitici, o fully coupled-simulations.

Nei metodi co-simulation i problemi elettrico e meccanico vengono

risolti separatamente, utilizzando codici differenti; l’accoppiamento fra gli

stessi e garantito da uno scambio di informazioni che avviene in determi-

nate fasi del calcolo, consentendo di determinare la soluzione di uno dei

due domini in un determinato passo o iterazione ricevendo come input la

soluzione dell’altro. In particolare, gli spostamenti del dominio meccanico

vengono trasmessi all’altra parte per l’aggiornamento delle coordinate nodali

della mesh elettrica; le forze elettriche invece, calcolate sul dominio elettrico

deformato, vengono successivamente applicate alla struttura, della quale si

determina la risposta (figura 4.1).

1La notazione e stata modificata rispetto a quella del capitolo 2 in funzione della

trattazione seguente.

4.2 Metodo Monolitico 66

I vantaggi dei metodi co-simulation sono molteplici. Innanzitutto per-

mettono di adottare differenti discretizzazioni nei due domini elettrico e mec-

canico: cio consente sia l’accoppiamento di differenti tipologie di discretiz-

zazione per le varie fisiche (FEM-FEM, FEM-BEM ([42], [43], [44])), sia

l’utilizzo di mesh diversamente rifinite. Nei casi piu complessi si possono

presentare alcune complicazioni nelle fasi di trasmissione dei dati da un ri-

solutore all’altro ([39], [40], [41]). In secondo luogo, come gia accennato, e

possibile utilizzare software specifici per la risoluzione di ciascuno dei due

problemi elettrico o meccanico.

Nelle fully-coupled simulations l’intero problema accoppiato viene ri-

solto in maniera monolitica da un unico codice ed entrambi i domini sono

discretizzati con la medesima mesh. Queste procedure comportano in genere

la linearizzazione delle equazioni del problema accoppiato.

4.2 Metodo Monolitico

Come detto, nell’approccio monolitico ([24]) i diversi problemi fisici vengo-

no risolti simultaneamente con una formulazione completamente accoppiata.

Il problema elettro-meccanico, intrinsecamente non-lineare, viene ricondot-

to ad una successione di passi lineari e affrontato in termini incrementali,

provvedendo a correggere la soluzione ad ogni iterazione.

Considerando il sistema (4.1), si definiscono le forze e le cariche interne

come:

f int = Kmechuu u (4.2)

qint = Kφφ(u)Φ (4.3)

Le equazioni di equilibrio del problema accoppiato diventano di conseguenza:

Felec(U, V )− Fint(U) = 0 (4.4)

dove

U =

[

u

Φ

]

F =

[

f

q

]

(4.5)


Attraverso uno sviluppo in serie di Taylor e possibile linearizzare attorno

alla generica configurazione di equilibrio (u,Φ) le forze e le cariche elettro-

statiche:

f elec(u,Φ, V ) = f elec(u,Φ)+∂f

∂u

∣∣∣∣(u,Φ)

du+∂f

∂Φ

∣∣∣∣(u,Φ)

dΦ+ϑ(du, dΦ) (4.6)

qelec(u,Φ, V ) = qelec(u,Φ)+∂q

∂u

∣∣∣∣(u,Φ)

du+∂q

∂Φ

∣∣∣∣(u,Φ)

dΦ+ϑ(du, dΦ) (4.7)

in cui si e considerato il caso generale in cui siano presenti delle cariche con-

centrate all’interno del dominio elettrico. Arrestando gli sviluppi al primo

ordine e sostituendo le espressioni (4.6) e (4.7) nel problema (4.1), o analo-

gamente nel (4.4), si ottiene la formulazione linearizzata dell’accoppiamento

elettro-meccanico:

KT∆U = ∆Felec (4.8)

dove ∆U rappresenta il vettore degli incrementi di spostamento e di poten-

ziale. Esplicitando i vari termini:

[

Kuu(u,Φ) Kuφ(u,Φ)

KTuφ(u,Φ) Kφφ(u)

]

︸︷︷︸

KT

[

∆u

∆Φ

]

︸︷︷︸

∆U

=

[

∆f(u,Φ, V )

∆q(u, V )

]

︸︷︷︸

∆Felec

(4.9)

dove:

-

Kuu = Kmechuu − ∂f elec(u.Φ)

∂u(4.10)

e la sottomatrice di KT della sola parte meccanica e Kmechuu e la matrice

di rigidezza meccanica;

Kuφ = −∂f elec(u,Φ)

∂Φ(4.11)

e il termine di accoppiamento. In particolare risulta Kφu = KTuφ, per

cui la matrice di rigidezza tangente K e simmetrica;


- Kφφ e la matrice di rigidezza elettrica;

- ∆f e ∆q sono rispettivamente gli incrementi delle forze e delle cariche

nodali equivalenti.

L’accoppiamento elettro-meccanico richiede l’utilizzo di un risolutore di tipo

iterativo per la risoluzione del sistema (4.8). Uno degli algoritmi piu diffusi

e senza dubbio quello di Newton-Raphson (NR), il cui schema e riportato

in figura 4.2. Ad ogni iterazione k, il sistema risolvente viene linearizzato

rispetto alla configurazione equilibrata dell’iterazione precedente k − 1, de-

terminando il valore dei residui ∆Fk:

∆Fk = Felec(Uk−1, V )− Fint(Uk−1) (4.12)

Utilizzando il sistema di equilibrio statico linearizzato (4.9), si calcola la

correzione da apportare alla soluzione dell’iterazione precedente k − 1, at-

traverso l’inversione della matrice di rigidezza tangente:

∆Uk = K−1T ∆Fk (4.13)

Si osserva che la matrice tangente KT deve essere aggiornata ad ogni ite-

razione, data la sua dipendenza dal vettore soluzione U. Da un punto di

vista computazionale il calcolo di KT rappresenta l’operazione piu onerosa

dell’intero approccio monolitico. La soluzione aggiornata:

Uk = Uk−1 +∆Uk (4.14)

deve essere sottoposta ad un controllo di convergenza del tipo:

‖Uk −Uk−1‖‖Uk−1‖

< ε (4.15)

dove ε e la tolleranza prefissata. Se la condizione (4.15) e soddisfatta il

processo iterativo si arresta e Uk da la soluzione del problema.

In prossimita della configurazione di equilibrio l’algoritmo di Newton-

Raphson presenta un ordine di convergenza quadratico. La buona velocita di

convergenza e piu che compensata dall’onere computazionale derivante dal


Figura 4.2: Diagramma di flusso: metodo monolitico.


calcolo ed inversione della matrice di rigidezza tangente ad ogni iterazione.

Per questa ragione e stata proposta una versione modificata del precedente

algoritmo (Newton-Raphson Modificato, o NRM, [37]): tale procedura

mantiene costante la matrice di rigidezza tangente durante il ciclo iterativo

e pari al valore iniziale di prima iterazione (nel caso dinamico quindi tale

matrice verra aggiornata una volta per ogni istante temporale).

Naturalmente l’algoritmo modificato ha il vantaggio di diminuire l’onere

computazionale garantendo la stessa soluzione fornita dalla versione origi-

naria; tuttavia il metodo NRM compie un numero maggiore di iterazioni

rispetto al classico NR e cio determina una minore velocita di convergenza.

Nella trattazione svolta si e considerato il problema meccanico stati-

co. Il passaggio al caso dinamico e comunque immediato: aggiungendo alla

prima delle (4.1) il termine dinamico, il sistema (4.1) si modifica nel modo

seguente:

Muuu+Kmechuu u = f elec(u,Φ, V )

Kφφ(u)Φ = q(u, V )(4.16)

Introducendo la linearizzazione, le equazioni di equilibrio dinamico possono

essere scritte nella seguente forma compatta:

M∆U+KT∆U = ∆F, (4.17)

dove

M =

[

Muu 0

0 0

]

. (4.18)

e la matrice delle masse tangente in cui l’unico sottoblocco non nullo e quello

relativo alla parte meccanica.

4.2.1 Esempio: metodo monolitico applicato ad un caso 1D

Si consideri il modello monodimensionale di figura 4.3 in cui i due domini

meccanico ed elettrico sono rappresentati da due bielle di lunghezza rispet-

tivamente Lm e Le. Nel caso piu generale il sistema globale e dotato di due

gradi di liberta: spostamento u del nodo comune e potenziale φ fra i due


Figura 4.3: Modello di riferimento monodimensionale.

nodi di estremita della biella elettrica. Si vuole studiare l’accoppiamento

elettro-meccanico in dinamica utilizzando l’approccio di tipo monolitico ap-

pena illustrato.

Le equazioni governanti sono fornite dal sistema (3.51), scritto per il generico

istante n: 2

Muuu = felec(u, φ) + F

qelec(u, φ) + qext = 0(4.19)

Si osservi che:

- si e implicitamente assunto che la parte meccanica sia gia stata discre-

tizzata temporalmente;

- la parte meccanica e soggetta sia a forzanti elettriche sia ai carichi

dinamici (F );

- la seconda equazione, che esprime l’equilibrio elettrico del sistema, puo

essere esplicitata:3

εφ

Le − u= qext (4.20)

dove qext sono le cariche esterne applicate sul contorno oppure dovute

2Onde evitare di appesantire la trattazione e stato omesso il pedice temporale n nelle

equazioni che seguono.3Si pone qelec = −εφ/(Le − u).


al voltaggio imposto mentre ε = ε0A, con ε0 permittivita dell’aria

(materiale di cui e composta l’ipotetica biella elettrica) e A area della

sezione delle bielle;

- la forza elettrica puo essere valutata utilizzando la formula canonica

valida per il caso monodimensionale:4

felec(u, φ) =1

2ε

φ2

(Le − u)2(4.21)

- il modello permette di evidenziare non solo l’accoppiamento fra i due

problemi meccanico ed elettrico (sia lo spostamento u che il potenziale

φ sono presenti in entrambe le equazioni del sistema (4.19)), ma anche

la non-linearita di cui si e discusso al capitolo 2. Infatti, sulla base

delle espressioni introdotte, il sistema (4.19) viene riscritto nel modo

seguente:

Muuu = 12ε

φ2

(Le−u)2+ F

ε φLe−u + qext = 0

(4.22)

in cui e evidente la dipendenza della forza elettrica da termini al

quadrato dipendenti da spostamento e potenziale.

Per far fronte a tale non-linearita e possibile applicare il metodo di

Newton-Raphson per il caso dinamico, che consente la riduzione del pro-

blema non-lineare accoppiato (4.22) ad un sistema lineare scritto in termini

incrementali. Si introduce quindi la linearizzazione delle quantita felec e

qelec intorno alla configurazione individuata dalla coppia di valori (u, φ), che

si riferiscono all’iterazione precedente nell’ambito dell’algoritmo NR:

felec(u, φ) = felec(u, φ) +∂f

∂u

∣∣∣∣(u,φ)

du+∂f

∂φ

∣∣∣∣(u,φ)

dφ+ ϑ(du, dφ) (4.23)

qelec(u, φ) = qelec(u, φ) +∂q

∂u

∣∣∣∣(u,φ)

du+∂q

∂φ

∣∣∣∣(u,φ)

dφ+ ϑ(du, dφ). (4.24)

4La formula e del tutto analoga a quella dei condensatori piani; il modello analizzato

infatti rappresenta il caso ideale in cui i piatti paralleli si concentrano nei due nodi di

estremita della biella elettrica.


Arrestando gli sviluppi in serie al primo ordine, si ottengono le espressioni

seguenti:

felec(u, φ) = felec(u, φ) +εφ

2

(Le − u)3du+

εφ

(Le − u)2dφ (4.25)

qelec(u, φ) = qelec(u, φ) +εφ

(Le − u)2du− ε

(Le − u)dφ (4.26)

Sostituendo le (4.25) e (4.26) nelle equazioni di equilibrio dinamico (4.19) si

ottiene:

M − εφ

2

(Le−u)3εφ

(Le−u)2

εφ(Le−u)2

− ε(Le−u)

[

du

dφ

]

=

[

0

−dqext

]

(4.27)

dove

KT =

M − εφ

2

(Le−u)3εφ

(Le−u)2

εφ(Le−u)2

− ε(Le−u)

(4.28)

e la matrice di rigidezza tangente linearizzata intorno alla configurazione

(u, φ). La presenza di termini extra-diagonali evidenzia l’accoppiamento

elettro-meccanico: il campo elettrico che si genera nel gap, induce una forza

di attrazione, che viene equilibrata dalla forza di richiamo elastica; mentre

il movimento del punto meccanico modifica la capacita e quindi la quantita

di cariche presente.

Si noti come l’approssimazione, inevitabilmente introdotta dalla lineariz-

zazione dei termini di accoppiamento, consente di scrivere un sistema di

equazioni algebriche lineari.

La soluzione del sistema (4.27) scritto in forma incrementale permette

di ricavare le correzioni da apportare alla configurazione di riferimento:

U = U+ dU =

[

u

φ

]

+

[

du

dφ

]

(4.29)

Il ciclo iterativo si arresta quando la condizione di convergenza risulta sod-

disfatta. Nel caso il test di convergenza abbia esito negativo, si applica

4.3 Metodi partizionati 74

una ulteriore linearizzazione intorno alla configurazione (u, φ) appena tro-

vata; la sequenza viene ripetuta fino a raggiungimento della convergenza.

La procedura iterativa viene ripetuta in ogni istante temporale dell’analisi

dinamica.

4.3 Metodi partizionati

4.3.1 Metodo staggered e metodo delle iterazioni simultanee

L’approccio monolitico illustrato nel paragrafo 4.2 prevede che ad ogni ite-

razione si debba non solo calcolare, ma anche invertire la matrice tangente;

questa operazione e, dal punto di vista computazionale, molto onerosa. E

anche per ovviare a questo problema che si sono sviluppati metodi alternati-

vi, denominati partizionati, caratterizzati dalla soluzione disaccoppiata delle

fisiche con approccio iterativo; essi consentono fra l’altro la discretizzazio-

ne indipendente dei domini elettrico e meccanico, i cui problemi associati

vengono risolti separatamente grazie ad appositi codici.

Una classe di metodi partizionati, detta staggered, prevede che la

matrice di rigidezza tangente del problema accoppiato (K), data dalla rela-

zione (4.9), venga suddivisa in una parte implicita ed una esplicita:

K = KI+KE =

[

Kuu(u,Φ) Kuφ(u,Φ)

0 Kφφ(u)

]

+

[

0 0

−KTu,φ(u,Φ) 0

]

(4.30)

Data la partizione di KT , il problema scritto in forma compatta risulta:

KIdU = dP−KEdU (4.31)

Esplicitando i vari termini si ottiene:

Kuu(u,Φ)du+Kuφ(u,Φ)dΦ = dPu

Kφφ(u)dΦ = dQφ +KTuφ(u,Φ)du

(4.32)

Per risolvere il sistema (4.32) e necessario impostare un procedimento

iterativo. Di seguito si riportano le operazioni principali dell’algoritmo alla

generica iterazione i-esima.


1. Predizione dupi sugli incrementi del campo di spostamenti rispetto al-

la configurazione (ui−1,Φi−1), attorno alla quale viene linearizzato il

sistema:

dupi = dui−1 (4.33)

2. Calcolo degli incrementi di voltaggio dΦi utilizzando la seconda equa-

zione del sistema (4.32):

dΦi = K−1φφ(ui−1)

[dQφ +KT

uφ(ui−1,Φi−1)dupi

](4.34)

3. Determinazione del campo di spostamenti all’iterazione corrente, so-

stituendo nella prima equazione del sistema (4.32) i valori dΦi:

dui = K−1uu (ui−1,Φi−1) [dPu −Kuφ(ui−1,Φi−1)dΦi] (4.35)

4. Controllo di convergenza:

‖ui − ui−1‖‖ui−1‖

< ε (4.36)

dove ε e la tolleranza fissata. Se il controllo ha esito positivo, (ui,Φi) =

(ui−1 + dui,Φi−1 + dΦi) e la soluzione del problema elettro-meccanico. In

caso contrario bisogna linearizzare nuovamente le equazioni rispetto alla

configurazione di equilibrio (ui,Φi) e si ripetono le operazioni dal punto

1.

Nella procedura staggered i problemi vengono risolti in modo disac-

coppiato. I passi di calcolo prevedono l’inversione delle matrici Kφφ e Kuu,

che hanno le stesse dimensioni dei singoli problemi elettrico e meccanico.

Un’altra possibile partizione della matrice di rigidezza tangente, conduce

alla formulazione del metodo delle iterazioni simultanee:

KT =

[

Kuu(u,Φ) 0

0 Kφφ(u)

]

+

[

0 Kuφ(u,Φ)

−KTuφ(u,Φ) 0

]

. (4.37)

Portando a secondo membro la parte esplicita della matrice KT , il problema


Figura 4.4: Procedura staggered ([22]).

diventa:

Kuu(u,Φ)du = dPu −Kuφ(u,Φ)dΦ

Kφφ(u)dΦ = dQφ +KTuφ(u,Φ)du

(4.38)

Il metodo delle iterazioni simultanee si articola nei punti riportati nel seguito.

1. Doppia predizione sia sugli incrementi di spostamento che su quelli di

voltaggio, rispetto alla configurazione (ui−1,Φi−1) attorno alla quale

viene linearizzato il sistema:

dupi = dui−1

dΦpi = dΦi−1

(4.39)

2. Stima simultanea degli incrementi dui e dΦi dal sistema (4.38):

dui = K−1uu (ui−1) [dPu −Kuφ(ui−1,Φi−1)dΦ

pi ]

dΦi = K−1φφ(uui−1,Φi−1)

[

dQφ +KTuφ(ui−1,Φi−1)du

pi

] (4.40)

3. Controllo di convergenza sulla soluzione totale:

‖Ui −Ui−1‖‖Ui−1‖

< ε (4.41)

Come nel caso del metodo staggered classico, il processo iterativo ter-

mina quando il test di convergenza ha esito positivo; in questa situazione Ui

e la soluzione del problema elettro-meccanico In caso contrario, si linearizza

nuovamente il sistema e si ritorna al punto 1.

In figura 4.4 e rappresentata la sequenza di operazioni di un semplice

solutore staggered. All’iterazione k si calcolano i valori nodali di potenzia-

le utilizzando la soluzione meccanica dell’iterazione k − 1. Successivamente


Figura 4.5: Diagramma di flusso: metodo staggered.

si determinano le forze elettrostatiche f elec, che vengono trasmesse all’al-

goritmo meccanico; risolto il problema meccanico, si ricavano i valori degli

spostamenti nodali uk, che vengono passati al risolutore elettrico, in modo

da determinare i nuovi potenziali nodali Φk. Queste operazioni devono es-

sere ripetute fino a convergenza. In caso la condizione non sia rispettata, la

geometria del dominio elettrico viene modificata imponendo gli spostamenti

meccanici nei nodi di interfaccia, e si ripetono le operazioni precedenti. La

convergenza viene raggiunta quando la correzione del vettore spostamenti

nodali e molto piccola. Si osserva infine che non viene assemblata alcuna

matrice di accoppiamento, ne la matrice di rigidezza tangente globale KT .

Il diagramma di flusso in figura 4.5 mostra in dettaglio le fasi dell’algoritmo

staggered per il problema elettro-meccanico accoppiato.

Naturalmente i vantaggi del metodo staggeded sono evidenti nel ca-

so si utilizzi un approccio co-simulation; i problemi elettrico e meccanico

sono risolti in modo disaccoppiato attraverso algoritmi particolari, con un

vantaggio computazionale evidente, poiche i due problemi sono lineari.

4.4 Voltaggio di pull-in 78

Figura 4.6: Ricerca del pull-in per potenziale imposto ([22]).

4.4 Voltaggio di pull-in

Nella progettazione dei microsistemi e fondamentale mantenersi lontani dalla

condizione di instabilita elettro-meccanica (pull-in), poiche una volta che le

armature di un condensatore vengono a contatto il funzionamento del dispo-

sitivo e irrimediabilmente compromesso. Per questo e necessario conoscere il

voltaggio di pull-in in corrispondenza del quale il sistema diventa instabile.

L’algoritmo iterativo di tipo staggered converge nel tratto stabile del-

la curva in figura 3.7, mentre non puo piu essere utilizzato per trovare la

soluzione del problema elettrostatico, quando il sistema diventa instabile.

Se l’algoritmo raggiunge la convergenza, questo significa che il valore di

pull-in e maggiore del voltaggio imposto e quindi la situazione di equilibrio e

stabile. Al contrario, nel caso l’algoritmo non converga, si puo concludere che

il voltaggio di pull-in e inferiore al valore di potenziale (Vi) imposto al sistema

(figura 4.6). La procedura per la stima del valore di voltaggio di instabilita

ha una forte dipendenza dalla robustezza dell’algoritmo. In prossimita della

situazione di pull-in, il problema diventa mal-condizionato; l’algoritmo diver-

ge prima che venga raggiunta la condizione di instabilita elettro-meccanica.

Per questa ragione, il metodo staggered classico sottostima il voltaggio di

pull-in.

Per la corretta valutazione del voltaggio di pull-in esistono altri algo-

4.5 Esempio: problema 1D 79

Figura 4.7: Ricerca del pull-in con il metodo di Riks-Crisfield.

ritmi tra cui quello di Riks-Crisfield ([47],[48]), che consentono di seguire il

ramo instabile della curva. La tecnica di Riks-Crisfield consiste nel fissare

la distanza tra due punti nel grafico (d,V) in figura 4.7. Si cerca un nuovo

punto di equilibrio imponendo il seguente vincolo aggiuntivo:

sk · sk = ∆S2i (4.42)

dove

sk =

[

∆di +∆dk

∆Vi −∆Vk

]

(4.43)

L’ aggiunta di un’equazione di vincolo permette di trattare i ∆V come in-

cognite, riuscendo ad oltrepassare il punto di pull-in e a descrivere anche il

ramo instabile della curva (figura 4.7).

4.5 Esempio: problema 1D

Si analizza ora l’esempio di figura 4.8, proposto in [25], utilizzando le due

procedure monolitica e staggered, opportunamente implementate in un co-


Parte elettrica (aria)

Geometria / proprieta valore

Le [µm] 1

ε0 [C2/N ·m2] 8.8542 · 10−12

V [Volt] 15

Parte meccanica

Lm [µm] 0.2

A [m2] 4 · 10−16

ρ [kg/m3] 2330

E [N/m2] 5000

ν 0

Tabella 4.1: Geometria e proprieta dei domini elettrico e meccanico del caso

1D.

dice Matlab; l’obiettivo e dimostrare che entrambe conducono alla stessa

soluzione e confrontare i tempi di analisi.

Il sistema e composto da un elemento meccanico e tre elettrici. La par-

te meccanica e idealmente concentrata nel punto 4, che rappresenta anche

l’interfaccia tra i domini; per questo motivo l’unico grado di liberta mec-

canico della struttura e dato dallo spostamento verticale di quel punto. In

seguito all’applicazione della differenza di potenziale, la massa viene pertur-

bata e, a causa dello spostamento u dello stesso nodo di interfaccia, anche il

dominio elettrico varia la propria configurazione.

La parte elettrica e discretizzata in tre elementi di trave di lunghezza

le = Le/3. Ai nodi 2 e 3 sono associati i due gradi di liberta elettrici del

sistema, mentre nei punti 1 e 4 viene imposto il voltaggio. Si ipotizza inoltre

che la deformazione del dominio elettrico riguardi esclusivamente l’elemento

compreso fra i nodi 3 e 4, che si configura quindi come elemento di interfaccia.

Con riferimento all’approccio monolitico, trattato nel paragrafo 4.2.1,

le relazioni che governano il problema sono date dal sistema (4.19). Esplici-

tando i termini, nel caso in esame si ottiene:

Muuu = 12ε

V 2

(Le−u)2+ F

KφΦ = qext

(4.44)

dove la matrice M , derivante dall’introduzione della discretizzazione tempo-


Figura 4.8: Esempio monodimensionale.

rale (Newmark), ha l’espressione seguente:

M =1

β∆t2M +K =

1

β∆t2ρALm

2+

EA

Lm(4.45)

M e la massa meccanica concentrata nel punto 4, K la rigidezza e β un

parametro dello schema di integrazione; F e il termine noto aggiuntivo,

dovuto anch’esso alla discretizzazione temporale.5

La risposta transitoria del sistema (4.44) puo essere trovata attra-

verso l’algoritmo di Newton-Raphson nel caso dinamico. Per comprende-

re meglio il procedimento all’istante generico n, si riporta la figura 4.9,

che rappresenta gli andamenti della forza elastica f int(u) (curva blu) ed

elettrica felec(u, V ) (curva verde) nell’esempio analizzato, al variare dello

spostamento u. L’algoritmo di NR si articola nei passi seguenti:

1. Partendo dalla configurazione equilibrata dell’iterazione precedente k−1, si calcola la forza elettrostatica e la sua derivata rispetto allo sposta-

mento, in modo da poter linearizzare, in quel punto, la curva della felec:

felec(u, V ) = feleci−1 (ui−1, V ) +

εV 2

(Le − ui−1)3(ui − ui−1) (4.46)

Si noti come il termine della linearizzazione (4.23), contenente le de-

5Si rimanda all’appendice B per la trattazione riguardante Newmark.


Figura 4.9: Algoritmo di Newton-Raphson (istante n).

rivate rispetto alle incognite di potenziale, sia nullo; infatti le forze

elettriche dipendono unicamente dal gap tra gli elettrodi (Le−ui−1) e

dal voltaggio esterno imposto (V ).

2. Si scrive la formulazione del problema linearizzato:

(

Muu − εV 2

(Le − ui−1)3

)

ui = feleci−1 − εV 2

(Le − ui−1)3ui−1 (4.47)

La soluzione ui vale:

ui =

(

Muu − εV 2

(Le − ui−1)3

)−1(

feleci−1 − εV 2

(Le − ui−1)3ui−1

)

(4.48)

3. Controllo di convergenza:

‖ui − ui−1‖‖ui−1‖

< ε (4.49)

Se il test ha esito positivo, ui e soluzione del sistema e si ripete il ciclo

di operazioni per l’istante successivo n+1 ; altrimenti si ritorna al passo

1.


Figura 4.10: Spostamento verticale nodo 4: soluzioni monolitica e staggered.

Per quanto concerne la soluzione del problema con approccio stag-

gered, si faccia riferimento al diagramma di flusso in figura 4.5. In primo

luogo viene risolto il problema elettrico ricavando i valori nodali di potenziale

(φ2, φ3). Essi dipendono dalla soluzione meccanica all’iterazione precedente,

in quanto la configurazione del dominio elettrico e cambiata in seguito allo

spostamento u imposto nel nodo 4 di interfaccia. Successivamente si calcola

la forza elettrostatica f elec con la relazione (3.27) e la si applica al nodo 4;

si risolve quindi il problema meccanico. Queste operazioni devono essere

ripetute fino a convergenza.

Nella tabella 4.2 vengono riportati i parametri utilizzati per la discre-

tizzazione temporale con il metodo di Newmark implicito. La figura 4.10

mostra l’andamento nel tempo dello spostamento del nodo 4 di interfaccia.

∆t [sec] 10−8

ttot [sec] 10−5

β 1/4

γ 1/2

Tabella 4.2: Discretizzazione temporale

Si osserva la perfetta sovrapposizione delle curve ottenute con i due

approcci monolitico e staggered, che quindi risultano equivalenti in termini


di soluzione restituita.

Si propongono alcune semplici verifiche per controllare la correttezza

dei risultati ottenuti con gli algoritmi implementati. Innanzitutto e possibile

confrontare la soluzione con lo spostamento statico della parte meccanica,

sollecitata da una forza pari al valore della forza elettrica iniziale applicata

nel punto di interfaccia 4. L’equazione di equilibrio statico risulta:

kustatico = felec −→ ustatico =

12(ε)

V 2

L2e

EALm

= 3.9825 · 10−8m (4.50)

dove

ε = ε0A = 8.85 · 10−12 · 4 · 10−16 = 3.54 · 10−27 C2

N(4.51)

e il valore di permittivita elettrica associato alla sezione di area A, ideal-

mente concentrata nel punto 4. Il valor medio intorno al quale la risposta

dinamica oscilla e ricavabile dalla figura 4.10:

umedio =max(u)−min(u)

2= 4.3605 · 10−8m (4.52)

La differenza tra i valori ottenuti dalle (4.50) e (4.52) e del 8.66%.

Si verifica ora che il periodo proprio della struttura T ricavato dall’anali-

si sia comparabile con il valore teorico atteso. Infatti per il semplice esempio

in esame e possibile calcolare T utilizzando la formula dell’oscillatore sem-

plice:

T =2π

ω(4.53)

dove ω e la pulsazione del sistema, definita come:

ω =

√

K

M(4.54)

Sostituendo la (4.54) nella (4.53), si ottiene:


Coordinate dei punti di picco

Picchi 1 2 3 4

t [10−6 · sec] 0.33 0.97 1.61 2.24

u [10−8 ·m] 8.7187 8.7212 8.7195 8.7178

Periodo di oscillazione

Intervalli 1-2 2-3 3-4

Periodo (T) [10−7 · sec] 6.400 6.300 6.300

Tabella 4.3: Periodo di oscillazione ottenuto dalla simulazione numerica.

Algoritmo Tempo di analisi

Monolitico 0.6298 sec

Staggered 0.3781 sec

Guadagno [-] −39.96%

Tabella 4.4: Confronto dei tempi analisi monolitico-staggered.

T =2π√

Km

=2π

√EALm

ρALm2

= 6.066 · 10−7 s (4.55)

La tabella 4.3 riporta i valori del periodo ricavati dalla simulazione numerica.

Si osserva che i valori ottenuti sono del tutto confrontabili con quello teorico,

trovato con il modello di oscillatore semplice.

Il grafico 4.11 mostra come, anche per quanto riguarda i valori di

potenziale elettrico, i codici monolitico e staggered restituiscano risultati

perfettamente sovrapponibili.

E interessante sottolineare come, gia nel semplice modello monodi-

mensionale analizzato, con il codice staggered si abbia un notevole risparmio

computazionale (tabella 4.4).

L’ultima analisi svolta riguarda la stima del voltaggio di instabilita

elettro-meccanica. Da un punto di vista analitico, il pull-in dinamico puo

essere calcolato con riferimento alla formula 3.82:

VDPI = 0.92 ·√

8L3e · k27ε

= 0.92 ·

√

8L3e · EA

Lm

27ε= 26.61V (4.56)

Inoltre, nella situazione di pull-in, ci si attende che il gap tra gli elettrodi


Figura 4.11: Potenziale nel nodo 2.

valga:

uDPI =1

2Le = 0.5µm. (4.57)

Nel grafico 4.12 si presentano i risultati numerici ottenuti attraverso l’analisi

di pull-in con l’algoritmo monolitico. Si noti come il codice riesca a cogliere

molto bene l’andamento della curva spostamento-voltaggio, fino al valore di

pull-in (differenza del 0.3%). In corrrispondenza del voltaggio di pull-in la

tangente alla curva tende all’orizzontale.


Figura 4.12: Pull-in dinamico nel caso monodimensionale.

Capitolo 5

La decomposizione dei

domini applicata al problema

elettro-meccanico

5.1 Introduzione

L’utilizzo dei metodi di decomposizione dei domini trattati nel capitolo 2

applicati al sistema oggetto di una simulazione (ad elementi finiti per esem-

pio) comporta dei vantaggi in termini di riduzione dei tempi di calcolo, so-

prattutto nel caso in cui l’analisi riguardi sistemi complessi e di notevoli

dimensioni; in questo modo e possibile ottenere la soluzione localmente su

ciascun sottodominio e, successivamente, ricostruire quella globale, risolven-

do un particolare problema di interfaccia. I problemi a piu campi (multi-field

problems), nei quali il sistema presenta diversi comportamenti fisici governa-

ti da equazioni differenti, rappresentano una categoria per la quale i tempi

di analisi sono in genere molto elevati: tipici esempi sono il problema termo-

elastico, l’interazione fluido-struttura ed il problema elettro-meccanico di cui

si e discusso nel capitolo 3. La loro complessita risiede nell’accoppiamento

o interazione fra i fenomeni fisici che si presentano, per cui in genere non

e possibile determinare in modo indipendente la soluzione associata ad una

delle fisiche senza dover determinare anche quella delle altre.

Nel capitolo 3 si e gia discusso dei metodi che consentono di determi-

nare la soluzione di un problema elettro-meccanico accoppiato con domini

non sovrapposti. In particolare sono emersi, gia dal semplice esempio mono-

88

5.2 Primo livello di decomposizione: contributifree e link 89

dimensionale, come peraltro era gia noto, i vantaggi della procedura di tipo

staggered nei confronti di un approccio di tipo monolitico al problema. Il

metodo staggered, infatti, comporta una automatica partizione del sistema

accoppiato originale in sottoproblemi, l’uno meccanico l’altro elettrico, che

vengono risolti alternativamente all’interno di un processo iterativo fino al

raggiungimento della convergenza. A partire da tale metodo, ci si propone

nel seguito di aumentare il livello di decomposizione del problema elettro-

meccanico: l’idea e quella di utilizzare gli algoritmi proposti in letteratura,

di cui si e gia trattato nel capitolo 2, per i quali sono stati dimostrati i van-

taggi computazionali nell’ambito dei problemi meccanici.

5.2 Primo livello di decomposizione: contributi

free e link

5.2.1 Considerazioni iniziali

Uno degli aspetti negativi della procedura di tipo staggered consiste nel-

la bassa velocita di convergenza perfino nei confronti del metodo monoliti-

co, che pure richiede un tempo maggiore per la risoluzione dei sistemi che

governano il problema.

E per far fronte a questa mancanza che il metodo staggered puo essere

riformulato secondo la filosofia dell’algoritmo di Gravouil e Combescure (GC

nel seguito) esposto nel paragrafo 2.4 (si veda [3]). Si consideri un generico

dominio di analisi Ω, composto per semplicita da due parti, l’una meccanica

Ωm e l’altra elettrica Ωe; sia inoltre Γ la frontiera comune alle due parti e

sia V il potenziale esterno applicato, costante durante l’analisi (figura 5.1.a).

Cosı come per il problema puramente meccanico, anche in questo caso la so-

luzione associata a ciascuno dei due domini puo essere espressa come somma

di due contributi free e link, il primo calcolato trascurando l’interazione fra

le due parti e il secondo, inteso come correzione da apportare al precedente,

determinato come effetto dell’altra parte su quella considerata. Introdotta

una discretizzazione spaziale per entrambe le parti meccanica ed elettrica, la

soluzione per ciascuno dei due sottodomini viene quindi scomposta nei due

contributi free e link :

um = ufreem + ulink

m (5.1)


Figura 5.1: Decomposizione elettro-meccanica (a) e puramente meccanica

(b).

Φ = Φfree +Φlink (5.2)

in cui um e Φ rappresentano i vettori delle incognite nodali, spostamenti per

la parte meccanica e potenziali per quella elettrica.

Per meglio comprendere il significato fisico dei contributi che compa-

iono nelle (5.1) e (5.2) si consideri l’analogia col problema puramente mecca-

nico risolto nella formulazione originale di Gravouil e Combescure; dato un

dominio Ω′

, esso si immagina suddiviso nelle due parti Ω′

1 e Ω′

2 che, diver-

samente dai due sottodomini Ωm e Ωe, sono caratterizzati dalla stessa fisica

(sono entrambi solidi deformabili)(figura 5.1). In questo caso, il calcolo del

contributo free dello spostamento avviene considerando le due parti come

indipendenti l’una dall’altra e soggette quindi ai soli carichi esterni variabili

nel tempo; la correzione link viene, invece, determinata per ciascuno dei due

sottodomini come effetto dell’interazione con l’altro, che nel caso meccanico

corrisponde all’applicazione delle forze di interfaccia, uguali ed opposte per il

principio di azione-reazione. Il calcolo del secondo contributo puo effettuarsi

solo dopo aver risolto il problema di interfaccia che, nell’esempio meccanico,

deriva dalla condizione di continuita degli spostamenti (o velocita, secondo

la trattazione di Gravouil e Combescure) all’interfaccia stessa.

Allo stesso modo si procede per l’esempio elettro-meccanico che diffe-

risce dal precedente per la diversa fisica dei due sottodomini; in particolare

occorre identificare le due situazioni free e link per ciascuna delle due parti.

Valgono le seguenti osservazioni:


- dal punto di vista del dominio meccanico, la situazione free (assenza

della parte elettrica) corrisponde al caso in cui esso e soggetto ai soli

carichi esterni o a quelli dovuti alla dinamica; d’altronde l’azione della

parte elettrica su quella meccanica si esplica attraverso le forze di Cou-

lomb calcolate sulla base della distribuzione di potenziale nel dominio

Ωe;

- dal punto di vista del dominio elettrico, sussiste una forte dipendenza

della configurazione (cioe forma, dimensioni e mesh) dello stesso Ωe

dagli spostamenti della parte meccanica, imposti lungo l’interfaccia Γ;

questi ultimi a loro volta dipendono dalle forze applicate, comprese

quelle elettriche (questa dipendenza reciproca rende necessario l’uti-

lizzo di un metodo di risoluzione di tipo iterativo, come peraltro gia

evidenziato in precedenza). Ne deriva quindi che il contributo free

della (5.2) coincide con la distribuzione di potenziale nella situazione

indeformata del dominio elettrico, cioe di spostamenti meccanici im-

posti nulli. Il secondo termine link corrisponde invece ad una sorta

di assestamento del potenziale dovuto alle deformazioni del dominio

elettrico;

- mantenendo l’analogia col caso puramente meccanico, e necessario in

questo caso risolvere un doppio problema di interfaccia: da una parte il

calcolo delle forze elettriche in funzione dell’andamento del potenziale e

dall’altra la risoluzione del cosiddetto problema pseudo-meccanico che

consente di determinare le nuove coordinate nodali della mesh elettrica

per effetto degli spostamenti meccanici imposti.

5.2.2 Formulazione del nuovo algoritmo

Sottodominio meccanico

Si consideri il sistema di equazioni di equilibrio dinamico ottenuto a valle

della discretizzazione spaziale della parte meccanica Ωm del dominio di figu-

ra (5.1.a), scritto nell’istante tn+1:

Mumn+1 +Kumn+1 = f extn+1 + f elecn+1(umn+1 ,Φ, V ) (5.3)

dove M, K e f extn+1 rappresentano rispettivamente la matrice delle masse, la


matrice di rigidezza ed il vettore dei carichi esterni all’istante tn+1, come e

noto; ad essi si aggiunge il contributo delle forze elettriche che traducono

l’azione del sottodominio Ωe su Ωm; in particolare, e stata messa in evidenza

la dipendenza di f elecn+1 dalla configurazione del dominio elettrico (dipendente

a sua volta dallo spostamento meccanico imposto um1), dalla distribuzione

di potenziale e, indirettamente, dall’attuazione V . Si introduca una discre-

tizzazione temporale e si adotti il metodo di Newmark per l’integrazione nel

tempo delle equazioni semi-discretizzate di equilibrio. Precisamente, combi-

nando le (2.55) e (2.56), si puo scrivere:

umn+1 = pumn +1

β∆t2(umn+1 − pumn) (5.4)

Quindi, sostituendo nel sistema (5.3):

( 1

β∆t2M+K

)

︸︷︷︸

M

umn+1 = f elecn+1(um,Φ, V ) + f extn+1 −Mpumn +Mumn︸︷︷︸

Fn+1

(5.5)

Mumn+1 = f elecn+1(um,Φ, V ) + Fn+1 (5.6)

Si introduce nel sistema ottenuto la decomposizione free-link attraverso la

(5.1) (si omette da questo punto in poi l’indice temporale n+ 1):

Mufreem + Mulink

m = f elec(um,Φ, V ) + F (5.7)

Sulla base delle considerazioni riportate nel paragrafo 5.2.1, il problema mec-

canico (5.7) puo essere scomposto nei due sottosistemi seguenti:

Mufreem = F problema free meccanico (in Ωm)

Mulinkm = f elec(um,Φ, V ) problema link meccanico (in Ωm)

(5.8)

Occorre osservare che il primo dei due sottosistemi (5.8) viene calcolato una

sola volta nel generico istante tn+1 in quanto il vettore dei termini noti

F dipende dal vettore dei carichi esterni all’istante tn+1 e dalla soluzione

1Il pedice m indica gli spostamenti nodali del dominio meccanico.


all’istante precedente tn. Viceversa, il secondo sottosistema viene risolto ad

ogni iterazione poiche l’accoppiamento rientra proprio nel vettore f elec.

Sottodominio elettrico

Si consideri il sistema 3.16 risolvente il problema elettro-statico nel dominio

Ωe:

Kφ(um)Φ = Q(um, V ) in Ωe (5.9)

in cui la matrice Kφ dipende dalla configurazione del dominio elettrico e

quindi dai cedimenti meccanici impostigli; il vettore dei termini noti a sua

volta dipende, nel caso piu generale, sia da um che dal potenziale imposto V .

Sulla base delle considerazioni del paragrafo 5.2.1 il problema free elettrico

e immediatamente scritto come:

K0φΦ

free = Q0(V ) problema free elettrico (in Ωe) (5.10)

Nella (5.10) la matrice di rigidezza elettrica K0φ e associata alla configurazio-

ne indeformata del dominio elettrico, per cui puo essere calcolata una volta

nota la geometria iniziale dello stesso. Infatti:

K0φ = Kφ(um = 0) (5.11)

Nel corso dell’analisi dinamica, il sistema (5.10) viene risolto un’unica volta,

essendo il problema elettro-statico stazionario.

Si consideri ora la generica situazione deformata del dominio elettrico.

La matrice di rigidezza elettrica puo essere espressa come segue:

Kφ(um) = K0φ +∆Kφ(um) (5.12)

Allo stesso modo le cariche nodali:

Q(um, V ) = Q0(V ) + ∆Q(um, V ) (5.13)


in cui ∆Kφ(um) e ∆Q(um, V ) sono le variazioni rispettivamente della ma-

trice di rigidezza e delle cariche nodali dovute alla deformazione del dominio

elettrico rispetto alla configurazione indeformata (0). Sostituendo le (5.2),

(5.12), (5.13) nella (5.9), si ottiene (omettendo le dipendenze delle grandezze

in gioco):

(

K0φ +∆Kφ

)(

Φfree +Φlink)

= Q0 +∆Q

Svolgendo il prodotto al primo membro:

(

K0φ +∆Kφ

)

Φlink +K0φΦ

free = Q0

︸︷︷︸

free

+∆Q−∆KφΦfree (5.14)

Tenendo conto del problema free (5.10), la (5.14) diventa:

(

K0φ +∆Kφ

︸︷︷︸

Kφ(um)

)

Φlink = ∆Q−∆KφΦfree

︸︷︷︸

∆Q

Kφ(um)Φlink = ∆Q in Ωe (5.15)

Nella (5.15) l’accoppiamento col problema meccanico e contenuto nella di-

pendenza della matrice Kφ dagli spostamenti meccanici imposti su Γ; nella

pratica tale matrice viene assemblata ad ogni iterazione sulla base degli spo-

stamenti nodali della mesh elettrica ue, ottenuti dalla risoluzione del pro-

blema pseudo-meccanico; il sistema (5.15) puo essere quindi scritto in modo

piu esplicito, come peraltro gia evidenziato al capitolo 3:

Kφ(ue(um))Φlink = ∆Q problema link elettrico (5.16)

La matrice ∆Kφ viene invece determinata come semplice differenza:

∆Kφ = Kφ −K0φ

Le (5.10) e (5.16) sono equivalenti al sistema (5.9).


Figura 5.2: Decomposizione di 1 livello: problemi free e link.

La figura (5.2) mostra in maniera sintetica i problemi che devono essere

risolti con il nuovo algoritmo.

In figura (5.3), invece, si riporta il diagramma di flusso del nuovo al-

goritmo proposto. I pedici n e n + 1 si riferiscono all’istante temporale

rispettivamente precedente e attuale; i e l’indice di iterazione. I pedici m ed

e applicati ai vettori u sono necessari per distinguere gli spostamenti nodali

meccanici di Ωm da quelli elettrici di Ωe, questi ultimi determinati risolvendo

il problema pseudo-meccanico. Il vantaggio di un algoritmo di questo tipo

non risiede tanto nella riduzione della dimensione dei sistemi da risolvere

(nessuno dei due domini elettrico e meccanico sono infatti stati fisicamen-

te suddivisi; i sistemi che compaiono non solo sono dello stesso ordine di

quelli del metodo staggered, ma addirittura aumentano in numero), quanto

nella velocita di convergenza: adottando un ciclo iterativo riferito alle sole

correzioni link, anche di piu ordini di grandezza inferiori rispetto alle cor-

rispondenti quantita free, il numero di iterazioni necessarie per giungere a

convergenza diminuisce rispetto alla procedura di tipo staggered.


Figura 5.3: Algoritmo di 1 livello: decomposizione secondo i due contributi

free e link.

5.3 Secondo livello di decomposizione: algoritmo GC applicato allarisoluzione del problema meccanico 97

Figura 5.4: Decomposizione della parte meccanica del problema di

riferimento.

5.3 Secondo livello di decomposizione: algoritmo

GC applicato alla risoluzione del problema mec-

canico

L’algoritmo di Gravouil e Combescure per la dinamica strutturale ([3]) di-

scusso nel capitolo 2 puo essere applicato alle due fasi di risoluzione del pro-

blema meccanico free e link (rispettivamente operazioni 2 e 6 dello schema

di figura (5.3)) della procedura proposta nel paragrafo 5.2.

Si consideri il dominio globale Ω costituito dai due domini rispetti-

vamente meccanico Ωm ed elettrico Ωe. Si immagini ora di suddividere

ulteriormente la parte meccanica in due sottodomini Ωm1 e Ωm2 (figura 5.4):

il problema meccanico e quindi ricondotto a due sottoproblemi associati ai

due sottodomini, a cui si aggiunge un problema di interfaccia.

5.3.1 Decomposizione del problema meccanico free

Con riferimento al passo 2 dello schema (5.3), il problema meccanico free e

dato da:

Mufreem = F in Ωm (5.17)

valido per l’intero dominio meccanico Ωm (le equazioni sono da intendersi


riferite all’istante tn+1). La matrice M e il vettore F hanno le espressioni

ottenute nella (5.5). Introducendo la decomposizione secondo l’algoritmo

GC, il sistema (5.17) viene sostituito dalle seguenti:

M1ufreem1 = F1 + LT

1 λfree in Ωm1

M2ufreem2 = F2 + LT

2 λfree in Ωm2

(5.18)

L’apice free e da riferirsi alla decomposizione di livello 1 (secondo le due fisi-

che). Si osserva nelle due equazioni di (5.18) la presenza a secondo membro

delle forze di interazione; LTj costituiscono gli operatori di interfaccia, cosı

come definiti nel paragrafo 2.4. Seguendo la strategia dell’algoritmo GC, la

soluzione puo essere espressa come somma di due contributi free e link :

ufreemj = u

free/freemj + u

free/linkmj j = 1, 2 (5.19)

I due termini della (5.19) si determinano risolvendo i corrispondenti sistemi

free e link :

Mjufree/freemj = Fj j = 1, 2

Mjufree/linkmj = LT

j λfree j = 1, 2

(5.20)

Per poter ottenere il secondo contributo dai sistemi costituiti dalla seconda

delle (5.20) e necessario determinare preliminarmente le forze di interfaccia

(moltiplicatori di Lagrange) risolvendo il corrispondente problema conden-

sato:

Hλfree = bfree (5.21)

dove H rappresenta l’operatore di condensazione definito nel modo seguente:

H =2∑

j=1

Lj(ajMj + bjKj)−1LT

j (5.22)

dove aj e bj sono due coefficienti dipendenti dallo schema di integrazione

adottato. Il termine noto bfree e dato da:

bfree = −2∑

j=1

Ljufree/freej (5.23)


Ovviamente la formulazione puo essere estesa al caso di decomposizione della

parte meccanica in Nm sottodomini; in questo caso, accanto a Nm sistemi

del tipo (5.20), occorre aggiugere nI problemi condensati all’interfaccia del

tipo (5.21), con nI numero delle interfacce fra i sottodomini meccanici.

I sistemi (5.20) e (5.21) sostituiscono il sistema monolitico (5.17). Lo

schema di figura 5.5 mostra come il passo 2 del diagramma di flusso (5.3)

possa essere risolto in maniera alternativa utilizzando appunto la decompo-

sizione della parte meccanica secondo l’algoritmo GC. E evidente in questo

caso che l’applicazione di tale algoritmo comporta la risoluzione di sistemi

di dimensione inferiore con conseguente risparmio di risorse computazionali.

Negli esempi del capitolo 5 verranno mostrati i risultati nel caso in cui quel

tipo di decomposizione sia applicata ad un problema che invece risulta ac-

coppiato.

5.3.2 Decomposizione del problema meccanico link

La decomposizione intesa alla maniera di Gravouil e Combescure puo essere

applicata anche per l’altro passo dell’algoritmo (5.3) che comporta la risolu-

zione di un problema meccanico, cioe il settimo. Il problema meccanico link

differisce dal corrispondente free affrontato al paragrafo 5.3.1 esclusivamente

per i carichi applicati (forze elettriche in luogo dei carichi esterni e dinamici)

e per il fatto che viene risolto ad ogni iterazione; nonostante queste diffe-

renze, la decomposizione introduce una serie di operazioni (alternative alla

risoluzione del sistema totale) che sono del tutto analoghe al caso free.

Si consideri il sistema risolto al passo 7 dell’algoritmo di 1 livello:

Mulinkm = f elec in Ωm (5.24)

Anche in questo caso il sistema (5.24) e riferito all’intero dominio meccanico

Ωm (si veda la figura 5.1). Considerando una suddivisione in N sottodomi-

ni, l’equilibrio puo essere scritto per ciascuno di essi introducendo anche il

contributo delle forze di interazione tra sottodomini meccanici, agenti all’in-

terfaccia:

Mjulinkmj

= f elecj +Λlinkj in Ωmj

, j = 1, . . . , Nm (5.25)


Figura 5.5: Decomposizione di 2 livello: schema di risoluzione del problema

free meccanico.

5.4 Terzo livello di decomposizione: metodo FETI applicato alproblema elettrico e pseudo-meccanico 101

I vettori Λlinkj , rappresentanti le forze di interazione del sottodominio j-

esimo, tengono conto dei contributi provienienti da tutte le interfacce tra

il sottodominio stesso e quelli adiacenti; nel caso piu generale infatti esse

possono essere piu d’una. La soluzione e considerata come somma di due

contributi:

ulinkmj

= ulink/freemj + u

link/linkmj j = 1, . . . , Nm (5.26)

I due termini vengono calcolati risolvendo, in ordine, i tre sistemi di seguito

riportati:

Problemi link/free : Mjulink/freemj = f elecj j = 1, . . . , Nm

Problemi di interfaccia link : Hiλlinki = blink

i i = 1, . . . , nI

Problemi link/link : Mjulink/linkmj = Λlink

j j = 1, . . . , Nm

(5.27)

La soluzione puo essere infine ricostruita tramite la (5.26). La figura 5.6

mostra le operazioni previste dalla decomposizione del problema meccanico

link, alternative alla risoluzione del sistema monolitico del passo 7 dell’algo-

ritmo di 1 livello.

5.4 Terzo livello di decomposizione: metodo FE-

TI applicato al problema elettrico e pseudo-

meccanico

Nei precedenti paragrafi e stato sviluppato, a partire dal metodo staggered

di cui si e discusso al capitolo 4, un algoritmo ottenuto introducendo diversi

livelli di decomposizione: secondo le fisiche (contributi free e link delle par-

ti elettrica e parte meccanica) e suddivisione della parte meccanica. Dopo

l’applicazione dell’algoritmo di Gravouil e Combescure al dominio Ωm, un

ulteriore passo avanti consiste nell’introduzione della decomposizione anche

per quei sottoproblemi dell’algoritmo (5.3) che devono essere risolti sul do-



link meccanico.


Figura 5.7: Decomposizione della parte elettrica del problema di riferimento.

minio elettrico Ωe: problema elettrico free, problema pseudo-meccanico e

problema elettrico link (rispettivamente passi 1, 4 e 6 dell’algoritmo di de-

composizione di 1 livello); il primo e risolto una sola volta, i restanti ad

ogni iterazione. Come gia emerso nel capitolo 3, questi tre problemi sono

stazionari: e quindi naturale fare riferimento al metodo FETI ([1]) esposto

al capitolo 2 nel momento in cui si procede ad una partizione del dominio

elettrico.

Con riferimento al consueto dominio Ω costituito da una parte mec-

canica Ωm ed una elettrica Ωe, si supponga di suddividere la parte elettrica

in due sottodomini Ωe1 e Ωe2 (figura 5.7); l’estensione al caso di Ne parti

e comunque immediata. Si analizzano separatamente i due casi di proble-

ma elettrico e pseudo-meccanico, ai quali il metodo FETI viene applicato in

maniera leggermente diversa.

5.4.1 Metodo FETI applicato al problema pseudo-meccanico

Introduzione

Come discusso nel capitolo 3, l’accoppiamento elettro-meccanico si presenta

anche in termini di deformazione del dominio elettrico per gli spostamenti in-

dotti dalla parte meccanica attraverso l’interfaccia Γ tra le due parti. Si e vi-

sto che i modi per tenerne conto sono essenzialmente due: re-meshing e mesh

deformation. Nel primo caso si tratta di definire una nuova mesh ogni qual-

volta si verifichi uno spostamento dei nodi di interfaccia; il secondo metodo,

invece, tiene conto della deformazione della mesh del dominio elettrico. Con


riferimento a quest’ultimo approccio il problema pseudo-meccanico consente

proprio di determinare gli spostamenti dei nodi della suddetta mesh.

Le caratteristiche del problema pseudo-meccanico sono gia state di-

scusse nel capitolo 3: si tratta di fatto di un problema meccanico fittizio,

risolto sul dominio elettrico, in cui sono assenti carichi esterni e nel quale

sono applicati esclusivamente dei cedimenti all’interfaccia Γ tra i domini Ωm

e Ωe. La risoluzione del problema pseudo-meccanico consente di determi-

nare gli spostamenti dei nodi della mesh elettrica e quindi di aggiornare le

coordinate nodali della stessa mesh per il calcolo della nuova Kφ.

Occorre sottolineare che il problema pseudo-meccanico viene risolto

in statica e soprattutto con riferimento alla configurazione indeformata del

dominio elettrico, per cui la corrispondente matrice di rigidezza Kelec viene

assemblata un’unica volta al di fuori del ciclo temporale dell’analisi dinami-

ca. Il problema pseudo-meccanico in forma monolitica, come si presenta al

passo 4 dell’algoritmo di 1 livello, risulta:

Kelecue = Fps in Ωe (5.28)

in cuiKelec e la matrice di rigidezza meccanica del dominio pseudo-meccanico,

Fps e un vettore di termini nulli, non essendo applicato alcun carico al domi-

nio Ωe2; infine, ue e il vettore degli spostamenti nodali della mesh elettrica.

Si ribadisce ancora una volta che la deformazione della parte elettrica non e

dovuta a carichi esterni ma a cedimenti che traducono l’azione degli sposta-

menti meccanici.

Formulazione del problema

Ora si consideri la suddivisione di Ωe nei due sottodomini Ω(1)e e Ω

(2)e . Si

immagini di aver introdotto preliminarmente una discretizzazione spaziale

per entrambi i sottodomini; la formulazione del metodo FETI, cosı come in-

trodotta nel paragrafo 2.3.3 consente di scrivere il seguente sistema in luogo

del (5.28):

2Ovviamente in seguito all’imposizione dei vincoli, costituiti dai cedimenti imposti, il

vettore avra dei termini diversi da zero; in questo modo dalla risoluzione del sistema (5.28)

si ottiene una soluzione diversa da quella nulla.


K(1)elecu

(1)e = F

(1)ps +B

(1)T

ps Λps in Ω(1)e

K(2)elecu

(2)e = F

(2)ps +B

(2)T

ps Λps in Ω(2)e

B(1)ps u

(1)e +B

(2)ps u

(2)e = 0 su Γe

(5.29)

dove con Γe si e indicata l’interfaccia fra i due sottodomini pseudo-meccanici.

Il precedente e del tutto analogo al (2.36). Alcune osservazioni:

- i primi due sistemi di (5.29) esprimono l’equilibrio meccanico per i due

sottodomini Ω(1)e e Ω

(2)e pensati come separati;

- i moltiplicatori di Lagrange Λps rappresentano le forze di interazione

fra le due parti; il principio di azione-reazione garantisce che tali forze,

applicate ai due sottodomini, siano uguali ed opposte;

- il terzo sottosistema di (5.29) esprime la continuita degli spostamenti

all’interfaccia Γe fra i due sottodomini pseudo-meccanici;

- infine, le matrici booleane B(j)ps , con j = 1, 2, sono del tutto analoghe

a quelle definite al capitolo 2, consentono cioe il passaggio dai gradi

di liberta del singolo sottodominio a quelli di interfaccia e viceversa;

si sottolinea di nuovo che, nonostante sia risolto su dominio elettrico,

il problema e di tipo meccanico a cedimento imposto, pertanto a cia-

scun nodo corrisponde un numero di incognite pari alla dimensione del

problema (2D o 3D).

Nel capitolo 2 e stato introdotto il concetto di matrice pseudo-inversa per

quei casi in cui la particolare decomposizione adottata determini la forma-

zione di domini cosiddetti fluttuanti, cioe non sufficientemente o per niente

vincolati. Immaginando, come nell’esempio a due sottodomini del paragrafo

2.3.3, che il dominio Ω(2)e sia fluttuante, si invertono le prime due delle (5.29)

e si sostituiscono nella terza, avendo cura di sostituire l’inversa di K(2)elec con

la sua pseudo-inversa e di aggiungere il contributo dei modi rigidi:

u(1)e = K

(2)−1

elec (F(1)ps +B

(1)T

ps Λps) in Ω(1)e

u(2)e = K

(2)+

elec (F(2)ps +B

(2)T

ps Λps) +R(2)ps αps in Ω

(2)e

(B(1)ps K

(1)−1

elec B(1)T

ps +B(2)ps K

(2)+

elec B(2)T

ps )Λps =

= −(

B(1)ps K

(1)−1

elec F(1)ps +B

(2)ps (K

(2)+

elec F(2)ps +R

(2)ps αps)

) su Γe

(5.30)


Accanto al sistema (5.30) occorre considerare l’ulteriore condizione di orto-

gonalita (2.43), adattata al problema pseudo-meccanico:

R(2)T

ps (F(2)ps +B(2)T

ps Λps) = 0 (5.31)

Dalla combinazione delle (5.30) e (5.31) si ottiene il seguente sistema, valido

per il caso di due sottodomini pseudo-meccanici:

u(1)e = K

(1)−1

elec (F(1)ps +B

(1)T

ps Λps)

u(2)e = K

(2)+

elec (F(2)ps +B

(2)T

ps Λps) +R(2)ps αps

[

(B(1)ps K

(1)−1

elec B(1)T

ps +B(2)ps K

(2)+

elec B(2)T

ps ) B(2)ps R

(2)ps

B(2)T

ps R(2)T

ps 0

][

Λps

αps

]

=

=

[

−(B(1)ps K

(1)−1

elec F(1)ps +B

(2)ps K

(2)+

elec F(2)ps )

−R(2)T

ps F(2)ps

]

(5.32)

Il sistema (5.32) e del tutto analogo al (2.44). Estendendo al caso generale

di N sottodomini di cui Nf fluttuanti, si ottiene il sistema analogo al (2.45):

u(s)e = K

(s)+

elec (F(s)ps +B

(s)T

ps Λps) +R(s)ps α

(s)ps s = 1, . . . , Ne

[

Sdps −Gps

−GTps 0

][

Λps

α3ps

]

=

[

dps

−e3ps

] (5.33)

in cui le quantita che compaiono nel (5.33) sono definite in modo analogo

alle (2.46-2.50), adattate al caso pseudo-meccanico.

La soluzione del problema pseudo-meccanico e ottenuta dapprima ri-

solvendo il problema di interfaccia (costituito dal secondo sottosistema delle

(5.33) e successivamente ricavando gli spostamenti nodali u(s)e attraverso le

prime equazioni dello stesso sistema.

Modi rigidi

La determinazione della matrice Rpss dei modi rigidi del sottodominio s puo

avvenire attraverso l’utilizzo del cosiddetto metodo geometrico-algebrico; es-


so e stato proposto da Farhat e Geradin in [16] ed e valido per problemi di

meccanica strutturale, quale e appunto il problema pseudo-meccanico, pur

essendo risolto su dominio elettrico.

Il metodo distingue sottodomini privi di vincoli e arbitrariamente vincolati.

STRUTTURA PRIVA DI VINCOLI

Si consideri il generico sottodominio Ω(s)e della partizione e si immagini che

esso sia completamente fluttuante. Poiche si vogliono determinare i mo-

di rigidi, tale sottostruttura viene considerata rigida. Come e noto dalla

Meccanica Razionale, il numero di moti rigidi consentiti ad un solido in-

deformabile e pari rispettivamente a 3 o 6, a seconda che il problema sia

bidimensionale o tridimensionale. La Meccanica Razionale suggerisce anche

il cosiddetto teorema di composizione delle velocita che, per il caso 2D, ri-

sulta (si consideri contemporaneamente la figura 5.8):

uM

vM

ωM

=

1 0 −(yM − yO)

0 1 (xM − xO)

0 0 1

α1

α2

α3

(5.34)

in cui:

- uM , vM , ωM rappresentano rispettivamente le velocita in direzione x e

y e la velocita angolare nel generico punto M(xM , yM ) di Ωs;

- α1, α2, α3 sono invece le velocita in direzione x e y e la velocita angolare

di un punto O(xO, yO) di riferimento della sottostruttura.

In forma compatta la (5.34) puo scriversi:

vM = RMα(s) (5.35)

Le relazioni (5.35) possono essere scritte per tutti i nodi N(s)n della mesh

della sottostruttura Ω(s)e ; si definisce quindi il vettore globale delle velocita

dei nodi stessi:

v3(s)ps = R3(s)

ps α(s)ps (5.36)

dove:


Figura 5.8: Teorema di composizione della velocita.

R3(s)ps = [R1, . . . , RN

(s)n

]T (5.37)

v3(s)ps = [v1, . . . ,vN

(s)n

]T (5.38)

α(s)ps = [α1, α2, α3]

T (5.39)

STRUTTURA VINCOLATA IN MODO ARBITRARIO

Si consideri la generica sottostruttura Ω(s)e vincolata in modo arbitrario; in-

dicando con nT e nV rispettivamente il numero totale dei gradi di liberta

associati agli N(s)n nodi della mesh ed il numero di quelli vincolati, si defini-

sce la seguente matrice di dimensioni nT × nV :

E(s)ps = [e1, . . . , enV

] (5.40)

in cui ei e il vettore booleano associato al grado di liberta vincolato i-esimo,

costituito da zeri e da un termine unitario in corrispondenza della posizione

i-esima:

ei = [0, . . . , 0, 1, 0, . . . , 0]T (5.41)

Occorre valutare nel caso generale se alla struttura sono consentiti moti ri-

gidi oppure se e sufficientemente vincolata da evitarli; in pratica, si tratta

di determinare l’eventuale presenza di moti ad energia nulla. Si definisce la


seguente matrice:

Z(s)ps = E(s)T

ps R3(s)ps (5.42)

di dimensioni nV × 3 nel caso di problema bidimensionale (come si assume

da questo momento in poi). L’operazione (5.42) estrae le righe di R3(s)ps cor-

rispondenti ai gradi di liberta vincolati. Si considera quindi il sistema:

Z(s)ps α

(s)ps = 0 (5.43)

Possono verificarsi due casi distinti:

- sottostruttura ben vincolata: la matrice Z(s)ps e non-singolare, di conse-

guenza l’unica soluzione del sistema (5.43) e α(s)ps = 0; alla sottostrut-

tura non sono consentiti moti rigidi;

- sottostruttura non vincolata: data la singolarita della matrice Z(s)ps ,

essendo il sistema (5.43) omogeneo, esso ammette sempre almeno una

soluzione; in particolare esistono ∞l(s) soluzioni, in cui l(s) rappresenta

il rango di Z(s)ps e, quindi, il numero di modi rigidi consentiti.

Come suggerisce l’Algebra Lineare, attraverso la decomposizione ai valori

singolari della matrice Z(s)ps e possibile ricavare il numero di modi rigidi:

Z(s)ps = USVT (5.44)

in cui:

- U e una matrice ortogonale di dimensione nV × nV ;

- V e una matrice ortogonale di dimensione 3× 3 (caso 2D);

- S e una matrice diagonale di dimensione nV × 3 avente lo stesso rango

di Z(s)ps ; in particolare essa ha tutti valori nulli a meno di quelli sulla

diagonale, maggiori o uguali a zero (autovalori).

Il passo successivo consiste nella determinazione degli l(s) vettori linearmente

indipendenti α(s)k (autovettori), con k = 1, . . . , l(s), che soddisfano il sistema


(5.43) nel caso di struttura mal vincolata. In particolar modo, essi rappre-

sentano il nucleo della matrice Z(s), secondo la definizione data nel capitolo

2. Una volta noti questi vettori, e possibile risalire agli spostamenti dei nodi

dell’intera mesh del sottodominio attraverso il teorema (5.34).

Si definisce la seguente matrice, di dimensioni 3× l(s) che raccoglie i vettori

α(s)k :

Q(s)ps = [α

(s)1 , . . . ,α

(s)

l(s)] (5.45)

Si risale infine al vettore dei modi rigidi effettivo della sottostruttura, a

partire dalla matrice definita dalla (5.37) e ottenuta ipotizzando in via pre-

liminare che lo stesso sottodominio fosse non vincolato:

R(s)ps = R3(s)

ps Q(s)ps (5.46)

Risoluzione del problema di interfaccia

Una volta determinate le matrici dei modi rigidi, e possibile risolvere il pro-

blema di interfaccia costituito dal secondo sottosistema delle (5.33). Come

gia esposto al paragrafo 2.3.3, le potenzialita dei metodi di decomposizione

dei domini sono esaltate dall’utilizzo di procedure di tipo iterativo per la

risoluzione del suddetto problema di interfaccia, che consentono di evitare

l’inversione del complemento duale di Schur. Nei codici implementati (ap-

pendice D) si e utilizzato il metodo del gradiente coniugato precondizionato;

si rimanda ancora una volta all’appendice A per approfondimento.3

In figura 5.9 e riportato la schema delle operazioni da effettuare per la riso-

luzione del problema pseudo-meccanico con il metodo FETI.

3In ambiente Matlab e disponibile una funzione che risolve sistemi lineari classici del tipo

Ax = b proprio attraverso il metodo del gradiente coniugato precondizionato: pcg.mat.



pseudo-meccanico.


5.4.2 Metodo FETI applicato al problema elettrico

Introduzione

Nell’algoritmo di 1 livello (5.2) i problemi elettrici free e link vengono risolti

rispettivamente nei passi 1 e 6; essi consentono di ricavare la distribuzione di

potenziale nella situazione di dominio elettrico indeformato e la correzione

dovuta alla deformazione del dominio elettrico stesso:

K0φΦ

free = Q0 problema free elettrico (in Ωe)

KφΦlink = ∆Q problema link elettrico (in Ωe)

(5.47)

Nel seguito, per la formulazione dell’algoritmo, si fara riferimento, in luogo

dei due (5.47), al generico problema elettrico KφΦ = Q; ovviamente la

procedura ottenuta potra essere utilizzata anche per i sistemi (5.47), che a

quel caso generale si riconducono.

Formulazione del problema

Si consideri di nuovo la partizione del dominio elettrico di figura 5.7. Es-

sendo stato proposto originariamente per problemi strettamente meccanici

([1]), occorre innanzitutto comprendere come il metodo FETI possa essere

esteso al problema elettrico; a tale scopo e necessario tenere in considerazio-

ne l’analogia fra le grandezze in gioco nei due tipi di problemi. Tralasciando

anche in questo caso la derivazione delle equazioni che governano il problema

decomposto, gia trattata nel paragrafo 2.3.3, il sistema a cui si giunge a valle

dell’introduzione di una discretizzazione spaziale per entrambi i sottodomini

elettrici e il seguente:

K(1)φ Φ(1) = Q(1) +B

(1)T

el Λel in Ω(1)e

K(2)φ Φ(2) = Q(2) +B

(2)T

el Λel in Ω(2)e

B(1)el Φ

(1) +B(2)el Φ

(2) = 0 su Γe

(5.48)

Alcune osservazioni:

- i primi due sistemi di (5.48) esprimono l’equilibrio elettrico per i due

sottodomini Ω(1)e e Ω

(2)e pensati come a se stanti;

- i moltiplicatori di Lagrange Λel, diversamente dal problema pseudo-

meccanico, hanno qui il significato fisico di cariche di interfaccia fra


i due sottodomini; per analogia col problema meccanico, un principio

di azione-reazione di tipo elettrico impone che quelle cariche debbano

essere complessivamente nulle su Γe (condizione che anche le forze di

interazione dovevano soddisfare nella versione originaria del FETI);

- un altro aspetto, di cui occorre tener conto, e la continuita del poten-

ziale, grandezza analoga agli spostamenti, attraverso l’interfaccia Γe;

tale condizione e espressa dalla terza di (5.48);

- le matrici booleane B(j)el , con j = 1, 2, pur mantenendo la loro con-

sueta funzione legata al passaggio dalla numerazione del sottodominio

a quella dell’interfaccia e viceversa, differiscono dalle B(j)ps del pseudo-

meccanico. Nel problema elettrico infatti ad ogni nodo e associata

un’unica incognita (il potenziale appunto), in luogo dei due sposta-

menti del caso meccanico 2D.

Il concetto di pseudo-inversa introdotto nel paragrafo 2.3.3 e estendi-

bile anche alla matrice di rigidezza elettrica di un sottodominio fluttuante.

Quest’ultimo concetto puo essere esportato, forse in modo un po’ astratto,

anche al problema elettrico: cosı come nel caso meccanico l’assenza di vin-

colo coincide con la mancata assegnazione di spostamenti sul contorno del

sottodominio, allo stesso modo nel caso elettrico si considerera fluttuante

quella parte sul cui contorno non e imposto alcun potenziale. Immaginando

quindi, come nell’esempio a due sottodomini del paragrafo 2.3.3, che il domi-

nio Ω(2)e sia fluttuante, si invertono le prime due delle (5.48) e si sostituiscono

nella terza:

Φ(1) = K(1)−1

φ (Q(1) +B(1)T

el Λe) in Ω(1)e

Φ(2) = K(2)+

φ (Q(2) +B(2)T

el Λel) +R(2)el αel in Ω

(2)e

(B(1)el K

(1)−1

φ B(1)T

el +B(2)el K

(2)+

φ B(2)T

el )Λel =

= −(

B(1)el K

(1)−1

φ Q(1) +B(2)el (K

(2)+

φ Q(2) +R(2)el αel)

) su Γe

(5.49)

Nel sistema (5.49) e stata introdotta la matrice pseudo-inversa K(2)+

φ ed il

contributo dei modi rigidi, racchiusi nel vettore R(2)el e combinati attraverso

il vettore αel, analogamente al caso meccanico4.

Il sistema risolvente il problema elettrico decomposto in Ne sottodomini di

4Tale modifica e stata aggiunta per semplice analogia con il problema meccanico; si

esploreranno nei paragrafi seguenti i veri significati delle quantita in gioco.


cui alcuni fluttuanti e quindi:

Φ(s) = K(s)+

φ (Q(s) +B(s)T

el Λel) +R(s)el α

(s)el s = 1, . . . , Ne

[

Sdel −Gel

−GTel 0

][

Λel

α3

el

]

=

[

del

−e3el

] (5.50)

Formalmente il problema elettrico decomposto tramite FETI e del tutto

analogo al caso meccanico. La differenza sostanziale fra i due casi risiede

pero nel significato assunto dai modi rigidi.

Modi rigidi

La ricerca della matrice R(s)el dei modi rigidi elettrici puo avvenire tenendo

conto di alcune considerazioni derivanti dall’analogia fra problema meccanico

e problema elettrico:

- cosı come i modi rigidi meccanici sono caratterizzati da energia di de-

formazione nulla, tale deve essere anche l’energia elettrica; questa con-

dizione traduce l’assenza di deformazione nel caso meccanico mentre

nel caso elettrico l’annullarsi del campo elettrostatico;

- una seconda analogia riguarda le incognite dei due problemi: se da

una parte i modi rigidi meccanici consentiti (nel caso 2D) sono le due

traslazioni in x e y e la rotazione ω, allo stesso modo nel caso elettrico

la situazione di modo rigido coincide con quella di potenziale ovunque

costante sul sottodominio (il potenziale costante non produce alcun

campo elettrico, che quindi rimane nullo);

- se nel problema meccanico bidimensionale sono possibili tutti casi in-

termedi tra quelli limite, costituiti da numero massimo di modi rigidi

consentiti, pari a 3, e da struttura ben vincolata con nessuna possibi-

lita di muoversi rigidamente, nel caso elettrico si possono presentare

unicamente le due situazioni di sottodominio vincolato (quando su una

parte del suo contorno e assegnato il potenziale) o fluttuante (non e

imposto alcun vincolo sul potenziale); non sussistono in questo caso

situazioni intermedie.

Guardando di nuovo al problema meccanico, si cerca per il caso elettri-

co una relazione analoga al teorema (5.34); si consideri un nodo O(xO, yO)


della mesh del sottodominio elettrico fluttuante Ω(s)e ; tale nodo e assunto

come riferimento e ad esso e associato un potenziale αO. Considerando l’i-

potesi di sottodominio elettricamente rigido (definizione data in base alle

considerazioni precedenti), il potenziale nel generico nodo M(xM , yM ) della

stessa mesh si puo scrivere come:

φM = αO (5.51)

cioe il potenziale, unica incognita nodale, e lo stesso per tutti i nodi della

mesh. Allora la (5.51) e analoga alla scrittura compatta (5.35), con RM = 1.

Si definisce il vettore globale dei potenziali degli N(s)n della mesh, come fatto

nella (5.36):

Φ3(s) = R3(s)el α(s) (5.52)

in cui

Φ3(s) = [φ1, . . . , φN(s)n

]T (5.53)

R3(s)el = [1, 1, . . . , 1, 1]T (5.54)

α(s) = αO (5.55)

L’analogia con il problema puramente meccanico ha consentito di giungere

alle precedenti relazioni. Valgono in particolar modo le considerazioni che

seguono:

- nella situazione in cui il sottodominio e vincolato (potenziale imposto

ad almeno uno dei nodi della mesh) la soluzione viene determinata

tramite la prima delle (5.50) una volta noti i moltiplicatori di Lagrange;

in questo caso non e presente il termine legato ai modi rigidi;

- se invece il dominio e fluttuante (da un punto di vista elettrico) si

osserva, sempre nella (5.50), che il secondo termine e dato da un vettore

di costanti uguali; cio significa che la soluzione in termini di potenziale

e ricostruita sommando ad una distribuzione di potenziale ‘locale’, data

dal primo contributo, un valore costante.



elettrico free.

Le figure (5.10) e (5.11) riportano i diagrammi di flusso riguardanti

rispettivamente la risoluzione del problema elettrico free e link con l’intro-

duzione del metodo FETI.

Nel presente capitolo si e giunti ad una sintesi dei due macro-argomenti

trattati nei primi due capitoli, cioe la decomposizione dei domini ed il proble-

ma elettro-meccanico. In particolare e stato ottenuto per gradi un algoritmo

che prevedesse l’applicazione di quei metodi nella risoluzione dei sistemi del

problema elettro-meccanico.

Il prossimo capitolo e dedicato ad analisi effettuate con i codici realizzati

ad hoc: verranno innanzitutto considerati alcuni casi prova caratterizzati

da geometria semplice e da maggiore governabilita, dopodiche si passera ad

analizzare problemi piu complessi; verranno pero sempre messi in eviden-

za i guadagni in termini di tempo di calcolo che comporta l’introduzione



elettrico link.


della decomposizione nei confronti della piu tradizionale procedura di tipo

staggered.

La figura 5.12 riporta uno schema riassuntivo dei livelli di decomposi-

zione che sono stati introdotti a partire dall’algoritmo staggered. I problemi

free e link sono da intendersi risolti su ogni sottodominio della suddivisione.


Figura 5.12: Schema riassuntivo dei livelli di decomposizione introdotti.

Capitolo 6

Applicazioni

6.1 Introduzione

Nel presente capitolo si studiano semplici esempi di problema elettro-meccani-

co, utilizzando gli algoritmi discussi nei capitoli precedenti: in particolare,

si fara riferimento da una parte alla procedura staggered, di dimostrata vali-

dita, e dall’altra ai nuovi algoritmi proposti nel capitolo 4 (decomposizione).

Per poter effettuare queste analisi sono stati realizzati dei codici in Matlabr

dei diversi algoritmi in versione 2D, al fine di ottenere in modo immediato un

confronto, sia in termini di soluzione che di tempi di calcolo, fra i vari meto-

di (l’attenzione e volta soprattutto al secondo aspetto, oltre che ovviamente

alla validazione degli algoritmi di decomposizione).

L’implementazione dei codici e stata necessaria poiche gli algoritmi

proposti sono di nuova concezione; i programmi realizzati ad hoc hanno

rappresentato quindi un primo strumento di convalida di cio che e stato

affermato e proposto nei capitoli precedenti, nonche un punto di partenza

per eventuali studi futuri.

Sono stati realizzati i seguenti tipi di codice:

- codice 0: staggered;

- codice 1: decomposizione di 1 livello (o semplice);

- codice 2: decomposizione di 2 livello (o doppia, con l’aggiunta della

decomposizione meccanica);

- codice 3: decomposizione di 3 livello (o tripla, con l’ulteriore e defini-

tiva introduzione della decomposizione della parte elettrica).

120

6.2 Casi prova 121

Figura 6.1: Microbridge (a) e Comb Drive (b).

La realizzazione del codice 3 ha preso spunto dal codice FETI di secondo

ordine realizzato in [17] per l’analisi di piastre: le subroutines sviluppate in

[17] sono state adattate sia al problema pseudo-meccanico che, soprattutto,

a quello elettrico.

Nelle analisi che seguono verranno considerati elementi finiti di tipo

CST: la scelta e ricaduta su questo tipo di elementi per la loro semplicita

e praticita, con la consapevolezza che la loro poverta richiede un maggior

grado di raffinamenti della mesh. I codici implementati sono comunque im-

postati in modo tale da consentire l’aggiunta di altre tipologie di elementi.

6.2 Casi prova

Per validare i codici implementati si considerano innanzitutto due problemi

a geometria semplice. La scelta di questi due esempi non e casuale. Come gia

affermato nei capitoli precedenti, il problema elettro-meccanico trova appli-

cazione nei cosiddetti MEMS (Micro Sistemi Elettro-Meccanici); nonostante

questi siano in genere caratterizzati da geometrie piu o meno complesse,

gli elementi costruttivi base di quei dispositivi possono essere realisticamen-

te ricondotti agli schemi elementari di trave doppiamente incastrata (figura

6.1.a) e di mensola (figura 6.1.b).

6.2.1 Trave a doppio incastro

Il primo modello analizzato e quello della trave doppiamente incastrata. In

figura 6.2 si riporta la geometria dell’esempio, mentre nelle tabelle 6.1 e

6.2 Casi prova 122

Figura 6.2: Geometria del primo caso prova.

6.2 le caratteristiche dei materiali e i parametri dell’analisi dinamica. Ana-

logamente agli esempi proposti nei capitoli precedenti per la formulazione

dell’accoppiamento elettro-meccanico e degli algoritmi di risoluzione, anche

in questo caso viene applicata una differenza di potenziale V fra un elettrodo

fisso ed uno deformabile (la trave) separati da un gap costituito da aria. Il

voltaggio imposto V e costante nel corso di tutta l’analisi dinamica (anda-

mento secondo la funzione di Heaviside) e pari a 30 Volt; e inoltre assente

lo smorzamento: si studiano quindi le oscillazioni della struttura, tenendo

conto dell’analogia con il caso di forzante meccanica applicata all’istante ini-

ziale t = 0 e mantenuta costante, in cui la risposta del sistema e ben nota.

Data la simmetria del problema, l’analisi si concentra su meta trave, con

Parte meccanica (elettrodo deformabile)

Proprieta Simbolo Valore Unita

Modulo elastico E 1011 Pa

Densita ρ 2330 kg/m3

Coefficiente di Poisson ν 0.3 -

Parte elettrica (gap di aria)

Permittivita elettrica ε0 8.854 · 10−12 C2/N ·m2

Tabella 6.1: Proprieta dei materiali delle parti meccanica ed elettrica.

l’introduzione di un vincolo pattino a pattino per tener conto dell’altra meta

mancante; la mesh utilizzata nelle analisi che seguono, realizzata con il pro-

gramma GiDr ([52]), e riportata in figura 6.3.

6.2 Casi prova 123

Parametri di Newmark β 0.25 -

(schema di integrazione implicito) γ 0.50 -

Passo temporale ∆t 10−9 sec

Tempo totale ttot 5 · 10−6 sec

Tabella 6.2: Parametri dell’analisi dinamica.

Figura 6.3: Mesh di riferimento.

Soluzione staggered

Innanzitutto e necessario confrontare la soluzione ottenuta con il codice stag-

gered, che rappresenta il riferimento per le soluzioni degli algoritmi successivi,

con quella ottenuta tramite codice commerciale Comsolr ([53]). Si assume

come parametro di confronto la massima freccia della trave (spostamento

verticale del nodo 1 di figura 6.3): se ne riporta l’andamento nella figura

6.4, in cui vengono messi a confronto i risultati ottenuti col codice 0 e con

Comsol. I due andamenti sono praticamente quasi sovrapposti.

Una seconda verifica consiste nel confrontare il valore medio delle oscillazio-

ni con la soluzione statica; si prevede infatti che il nodo considerato oscilli

all’incirca intorno alla sua posizione di equilibrio statico. I due valori sono

riportati nella tabella 6.3, in cui si e indicato con umedio il valore medio del-

le oscillazioni e con ustatico la soluzione statica: i due valori sono del tutto

confrontabili (a meno di effetti dinamici).

6.2 Casi prova 124

Figura 6.4: Confronto della soluzione staggered con quella ottenuta con il

programma Comsol.

Soluzione statica ustatico −1.72309 · 10−3 µm

Valore medio delle oscillazioni umedio −1.70357 · 10−3 µm

Errore 1.13%

Tabella 6.3: Confronto con la soluzione statica.

Un ulteriore confronto puo essere effettuato fra il periodo proprio della

struttura e i risultati numerici dedotti dalla figura 6.4, ottenuti calcolando la

distanza fra picchi adiacenti. La figura 6.5 riporta il primo modo di vibrare

della trave, di tipo flessionale; il periodo corrispondente puo essere valuta-

to a valle della risoluzione di un problema di autovalori (qui effettuata con

pdetool di Matlab). Esso risulta:

ω1 =√

λ1 = 6.066 · 106rad/sec

T1 =2π

ω1= 1.036 · 10−6sec

La tabella 6.4 riporta i valori del periodo ricavati dall’analisi effettuata col

6.2 Casi prova 125

Figura 6.5: 1 modo di vibrare della trave doppiamente incastrata.


Picchi 1 2 3 4 5

t [10−6 · sec] 0.503 1.478 2.463 3.471 4.440

u [10−3 · µm] −3.4637 −3.4914 −3.4470 −3.4781 −3.4756


Intervalli 1-2 2-3 3-4 4-5

Periodo [10−6 · sec] 0.975 0.985 1.008 0.969


codice 0. Si osserva che i valori ottenuti sono del tutto confrontabili con

quello teorico.

I confronti precedenti di fatto provano che la soluzione ottenuta con il

codice 0 corrisponde a quella dell’algoritmo (staggered). Essa costituisce la

soluzione di riferimento per i risultati ottenuti con gli algoritmi di decompo-

sizione che verranno mostrati nei paragrafi seguenti.

6.2 Casi prova 126

Decomposizione di 1 livello: confronto con la soluzione staggered

Il problema della trave a doppio incastro viene ora analizzato utilizzando il

codice 1 (decomposizione semplice). I dati riguardanti geometria, materiali e

dinamica sono quelli di riferimento per l’intero esempio della trave a doppio

incastro.

Il confronto viene qui effettuato in termini sia di freccia del nodo 1

del dominio meccanico, che di potenziale del nodo 2 del dominio elettrico

(figura 6.3), di coordinate rispettivamente (60µm; 4µm) e (60µm; 2µm).

Le figure 6.6 e 6.7 riportano gli andamenti dei due parametri di confronto

valutati con i due diversi codici. Si osserva che le soluzioni ottenute, sia in

termini di spostamento che di potenziale, sono perfettamente sovrapposte:

l’algoritmo di decomposizione di 1 livello fornisce quindi la stessa soluzione

dell’algoritmo staggered. La perfetta corrispondenza fra i due algoritmi puo

essere dimostrata anche attraverso il confronto di altre grandezze in gioco,

come il campo elettrico nella parte elettrica o gli sforzi in quella meccanica.

Si riportano i confronti fra le seguenti grandezze, i cui andamenti sono stati

ottenuti coi due codici 0 e 1:

- spostamenti in direzione x (figura 6.9);

- spostamenti in direzione y (figura 6.10);

- sforzi σx (figura 6.11);

- sforzi σy (figura 6.12);

- sforzi τxy (figura 6.13);

- andamento del potenziale elettrico (figura 6.14);

- andamento del campo elettrico in direzione x (figura 6.15);

- andamento del campo elettrico in direzione y (figura 6.16);

Queste quantita sono riferite all’istante t = 0.5 · 10−6 sec, in corrispondenza

del primo picco. Data la simmetria del problema, le deformate e i contour

plot sono relativi solo alla meta sinistra della struttura (in figura 6.8 si mo-

stra comunque l’andamento qualitativo della deformata dell’intera trave).

Dalle figure 6.9-6.16 emerge che le grandezze calcolate coi due diversi algo-

ritmi sono del tutto coincidenti.

6.2 Casi prova 127

Figura 6.6: Confronto della freccia massima ottenuta con codici 0 e 1.

Figura 6.7: Confronto del potenziale di riferimento ottenuto con codici 0 e

1.

6.2 Casi prova 128

Figura 6.8: Deformata qualitativa della trave a doppio incastro.

Avendo utilizzato elementi CST, le deformazioni (e quindi gli sforzi) e

i campi elettrici risultano costanti sul singolo elemento finito; i valori ottenu-

ti sono stati interpolati linearmente (tramite la funzione pdeprtni di Matlab

che consente di ottenere un andamento continuo e non a scalini di quelle

quantita).

Si confrontano a questo punto i tempi di analisi; essi sono stati valutati

sulla stessa unita di calcolo e sotto le medesime condizioni. La tabella 6.5

mostra i tempi che sono stati necessari per l’analisi di tipo staggered e con al-

goritmo di 1 livello, ricordando che il tempo totale e pari a ttot = 5 ·10−6sec

ed il passo adottato ∆t = 10−9sec.


Staggered 100160 sec

Decomp. 1livello 77616 sec


Tabella 6.5: Confronto dei tempi analisi staggered -decomposizione 1livello.

A parita di condizioni (geometria, materiali, parametri dinamici, elaborato-

re, . . . ) e di soluzioni, l’algoritmo di decomposizione di 1 livello presenta

un vantaggio computazionale consistente nei confronti della procedura

staggered e, a maggior ragione, di quella monolitica (come messo in evidenza

nell’esempio 1D del capitolo 3). Come accennato nel capitolo 4, la ragione

di tale risparmio di tempo risiede essenzialmente nella maggiore velocita di

convergenza del nuovo algoritmo: infatti, il ciclo iterativo riguarda le sole

6.2 Casi prova 129

Figura 6.9: Confronto degli spostamenti in direzione x (spostamenti in µm).

6.2 Casi prova 130

Figura 6.10: Confronto degli spostamenti in direzione y (spostamenti in µm).

6.2 Casi prova 131

Figura 6.11: Confronto degli sforzi σx (in N/m2).

6.2 Casi prova 132

Figura 6.12: Confronto degli sforzi σy (in N/m2).

6.2 Casi prova 133

Figura 6.13: Confronto degli sforzi τxy (in N/m2).

6.2 Casi prova 134

Figura 6.14: Confronto del potenziale (in V olt).

6.2 Casi prova 135

Figura 6.15: Confronto del campo elettrico in direzione x (in V/m).

6.2 Casi prova 136

Figura 6.16: Confronto del campo elettrico in direzione y (in V/m).

6.2 Casi prova 137

Figura 6.17: Suddivisione del dominio meccanico in 2 sottodomini.

correzioni link degli spostamenti e dei potenziali, inferiori anche di piu ordini

di grandezza rispetto alle corrispondenti quantita free.


Nel presente paragrafo si mostrano i risultati ottenuti introducendo la de-

composizione della parte meccanica. Si considera, in primo luogo, una di-

visione in due sottodomini 1 (figura 6.17). Si riportano nelle figure 6.18 e

6.19 i confronti effettuati sui due valori di freccia massima e di potenziale

intermedio (associati ai consueti nodi 1 e 2, figura 6.17), ottenuti con l’al-

goritmo di decomposizione di 2 livello e con la procedura staggered. Dai

grafici si osserva che gli andamenti di quei due parametri sono coincidenti.

L’algoritmo con decomposizione di 2 livello (dec.semplice + dec.meccanica)

fornisce quindi la stessa soluzione della procedura staggered di riferimento.

I grafici di figura 6.20 riportano il confronto fra lo spostamento e la

velocita verticali di un nodo di interfaccia (nodo 3 di figura 6.17), conside-

rato come appartenente ai sottodomini 1 e 2 rispettivamente. Si osserva che

i due andamenti sono pressoche sovrapposti, per cui la continuita degli spo-

stamenti e delle velocita e verificata. In particolare, si osserva che la velocita

si annulla in corrispondenza dei picchi di spostamento, quando cioe cambia

il verso dello stesso spostamento; risulta, invece, massima quando lo spo-

1Si ricorda che il problema della trave doppiamente incastrata e stato studiato per

meta, data la simmetria della struttura; quindi la suddivisione riguarda solo meta trave.

6.2 Casi prova 138



2.

6.2 Casi prova 139

stamento e nullo, cioe quando il dispositivo ritorna nella sua configurazione

indeformata durante il suo moto oscillatorio.

Il confronto fra i tempi di analisi e riportato in tabella 6.6. Si os-

serva che l’introduzione della decomposizione meccanica (due sottodomini)

comporta un vantaggio in termini di tempo di analisi leggermente maggiore

rispetto alla sola decomposizione di 1 livello (2% circa in piu di guadagno).



Decomp. 2livello (2 Sd) 76024 sec


Tabella 6.6: Confronto dei tempi analisi staggered -decomposizione 2livello

(2 sottodomini).

Si considerino ora i casi di suddivisione del dominio meccanico in 3

e 4 sottodomini uguali, ottenute attraverso una partizione verticale dello

stesso. Tralasciando gli andamenti della freccia massima e del potenziale

di riferimento, che risultano anche per questi casi sovrapposti alla soluzio-

ne staggered, si riportano in figura 6.21 i tempi di analisi corrispondenti ai

diversi casi. L’analisi e stata effettuata su un tempo ridotto (10−6 sec in

luogo di 5 · 10−6 sec). Si osserva dalla figura 6.21 che il tempo di analisi non

subisce modifiche sostanziali all’aumentare del numero di sottodomini mec-

canici. L’introduzione della decomposizione meccanica in questo esempio,

quindi, non produce vantaggi apprezzabili a livello computazionale: i risul-

tati ottenuti sono, infatti, riferiti alla mesh di figura 6.3 uniformemente fitta,

per la quale la risoluzione del problema elettrico ha un peso maggiore in ter-

mini di onere di calcolo, essendo appunto la parte elettrica piu estesa. Ci si

aspetta quindi che i tempi di analisi subiscano una riduzione consistente con

l’introduzione del 3 livello di decomposizione (della parte elettrica), come

peraltro verra messo in evidenza al paragrafo successivo. Nel caso analizza-

to al paragrafo 6.3.2 si mettera in evidenza, inoltre, come la decomposizione

meccanica sia vantaggiosa nel caso, non raro, in cui la mesh elettrica sia piu

rada rispetto a quella meccanica.

6.2 Casi prova 140

Figura 6.20: Spostamenti e velocita del nodo 3 di interfaccia fra i domini 1

e 2.

6.2 Casi prova 141

Figura 6.21: Confronto dei tempi di analisi con diverso numero di

sottodomini meccanici.


Si introduce infine la decomposizione del dominio elettrico, che si aggiunge

alla precedente suddivisione della parte meccanica in due sottodomini, se-

condo l’algoritmo di decomposizione di 3 livello presentato al capitolo 5. Si

consideri inizialmente la suddivisione del dominio elettrico in quattro par-

ti, mostrata in figura 6.22. Ancora una volta si riporta, nelle figure 6.23 e

6.24, il confronto della freccia massima e del potenziale di riferimento con

la procedura staggered ; gli andamenti sono ancora una volta sovrapposti.

Si riporta nella tabella 6.7 il confronto dei tempi di analisi. A parita di

soluzione, l’algoritmo di 3 livello risulta molto vantaggioso anche solo con

una suddivisione base come quella di figura 6.22.

Si considerino ora i casi di partizione verticale del dominio elettrico

rispettivamente in 3 e 4 sottodomini (per un totale rispettivmaente di 6 e

8 sottodomini elettrici, avendo mantenuto una partizione in orizzontale di

2). La figura 6.25 riporta un confronto dei tempi di analisi rispettivamente

con 4, 6 e 8 sottodomini elettrici (e, corrispondentemente, 2, 3 e 4 sotto-

domini meccanici), ottenuti su un tempo pari a ttot = 10−6sec. I tempi

6.2 Casi prova 142

Figura 6.22: Suddivisione della parte elettrica in 4 sottodomini e della parte

meccanica in 2 sottodomini.



Decomp. 3livello (2 sd Mecc. + 4 sd El.) 69605 sec



(2 sd meccanici, 4 sd elettrici).

rilevati dimostrano che una maggiore suddivisione della parte elettrica com-

porta vantaggi ancor piu evidenti rispetto alla partizione di figura 6.22; cio

conferma quanto affermato alla fine del paragrafo precedente: il vantaggio

e apprezzabile per la decomposizione della parte elettrica proprio perche al

problema elettrico e associato un onere di calcolo maggiore, nel caso in cui

si utilizzi una mesh come quella di figura 6.3, per cui l’introduzione della

decomposizione di 3 livello ha un piu evidente effetto benefico.

6.2 Casi prova 143

Figura 6.23: Confronto della freccia massima ottenuta con codici 0 e 3 (2

sottodomini meccanici, 4 sottodomini elettrici).

Riassunto dei tempi di analisi

Si riassumono nella tabella 6.8 i tempi di analisi e i guadagni percentuali

ottenuti con i diversi livelli di decomposizione, per l’esempio della trave a

doppio incastro (ttot = 5 · 10−6sec; ∆t = 10−9sec). I tempi relativi al caso

di decomposizione di 2 e 3 livello si riferiscono rispettivamente ad una

suddivisione del dominio meccanico in 2 parti e di quello elettrico in 4 parti.

Il guadagno percentuale e relativo all’algoritmo staggered.

Algoritmo Tempo di analisi Guadagno %

Staggered 100160 sec -

Decomp. 1livello 77616 sec −22.51%



Tabella 6.8: Riassunto dei tempi analisi ottenuti con i vari algoritmi (trave

a doppio incastro).

6.2 Casi prova 144


3 (2 sottodomini meccanici, 4 sottodomini elettrici).


sottodomini elettrici e meccanici.

6.2 Casi prova 145

6.2.2 Trave a mensola

Si studia in questa sezione il secondo caso prova dato dalla trave a menso-

la di figura 6.26. Geometria, materiali e parametri dinamici sono gli stessi

del caso prova precedente (si vedano le tabelle 6.1 e 6.2) cosı come la mesh

adottata (figura 6.3). La differenza fra i due casi risiede nelle condizioni di

vincolo: in questo secondo caso l’estremo di destra della trave risulta libero.

Il potenziale imposto e costante e pari a 30 Volt. Si riportano nel seguito i

Figura 6.26: Geometria del secondo caso prova.

risultati ottenuti, procedendo in modo del tutto analogo al caso della trave

doppiamente incastrata.

Soluzione staggered

Si assume come parametro di confronto lo spostamento verticale del nodo 1

di figura 6.26. La figura 6.27 mostra l’andamento delle curve ottenute con

il codice 0 (staggered) e con Comsol ; si osserva che i due andamenti sono

pressoche sovrapposti. La tabella 6.9 riporta il confronto fra valore medio

delle oscillazioni e soluzione statica.

Soluzione statica ustatico −5.1723 · 10−3 µm

Valore medio delle oscillazioni umedio −4.9344 · 10−3 µm

Errore 4.59%

Tabella 6.9: Confronto con la soluzione statica.

6.2 Casi prova 146

Figura 6.27: Confronto della soluzione staggered con quella ottenuta con il

programma Comsol.

Il primo modo di vibrare della mensola e riportato in figura 6.28. Vie-

ne calcolato il primo periodo modale:

ω1 =√

λ1 = 3.811 · 106rad/sec

T1 =2π

ω1= 1.649 · 10−6sec

Tale valore deve essere confrontato con le stime riportate nella tabella 6.10.

Si osserva che i valori risultano del tutto confrontabili.


Picchi 1 2 3

t [10−6 · sec] 0.775 2.335 3.952

u [10−3 · µm] −10.4795 −10.3841 −10.3726


Intervalli 1-2 2-3

Periodo [10−6 · sec] 1.560 1.617


6.2 Casi prova 147

Figura 6.28: 1 modo di vibrare della trave a mensola.


Le figure 6.29 e 6.30 riportano il confronto fra spostamenti e potenziali di

estremita2 nei due casi staggered e decomposizione di 1 livello. Ancora una

volta le soluzioni risultano sovrapposte. Le figure 6.31-6.34 mostrano l’an-

damento degli spostamenti nelle due direzioni x e y, degli sforzi σx, σy, τxy,

del potenziale φ e dei campi elettrici Ex e Ey. Gli andamenti e i contour

plot riportati sono validi per entrambi gli algoritmi con cui sono stati de-

terminati. L’istante a cui essi si riferiscono e in prossimita del primo picco

(a t = 0.7 · 10−6sec). Infine, la tabella 6.11 riporta il confronto fra i

tempi delle analisi staggered con l’algoritmo di decomposizione di 1 livello

(ttot = 5 · 10−6sec; ∆t = 10−9sec). Il secondo caso prova conferma l’entita

del guadagno di tempo percentuale nel passaggio dal primo al secondo algo-

ritmo, che si attesta intorno al 20%.

2Oltre allo spostamento del nodo 1, si assume come parametro di confronto il potenziale

a meta altezza del gap in corrispondenza del nodo 2 di figura 6.26.

6.2 Casi prova 148



1.

6.2 Casi prova 149

Figura 6.31: Spostamenti in direzione x e y (in µm).

6.2 Casi prova 150

Figura 6.32: Sforzi normali in direzione x e y (in N/m2).

6.2 Casi prova 151

Figura 6.33: Sforzi di taglio (in N/m2) e potenziale (in V olt).

6.2 Casi prova 152

Figura 6.34: Campi elettrici in direzione x e y (in V/m).

6.2 Casi prova 153


Staggered 95378 sec



Tabella 6.11: Confronto dei tempi di analisi staggered -decomposizione

1livello.


Si immagini di suddividere la mensola in due parti uguali attraverso una par-

tizione verticale, analogamente a quanto fatto nel primo caso prova. Dando

per assodata la sovrapposizione fra gli andamenti della freccia massima e

del potenziale intermedio, e quindi l’equivalenza fra gli algoritmi in termini

di soluzione restituita, ci si concentra sui tempi di analisi. La tabella 6.12

rappresenta il confronto fra i tempi di analisi staggered e decomposizione di

2 livello (con 2 sottodomini meccanici). Ancora una volta il vantaggio deri-

vante dalla decomposizione meccanica non e apprezzabile con 2 sottodomini

meccanici (in questo caso il tempo di analisi e addirittura maggiore, seppur

di poco, a quello della decomposizione di 1 livello); cio e dovuto, come gia

anticipato, al maggior peso della risoluzione del problema elettrico all’inter-

no dell’algoritmo. La figura 6.35 mostra invece i tempi di analisi in funzione

del numero di sottodomini, per la quale valgono le stesse osservazioni del

primo caso prova.


Staggered 95378 sec

Decomp. 2livello (2 Sd) 75667 sec



(2 sottodomini).


Si introduce infine la decomposizione del dominio elettrico. Anche per questo

caso, si consideri innanzitutto la partizione del dominio elettrico in (2 × 2)

parti. A parita di soluzione (non si riportano gli andamenti di freccia mas-

6.2 Casi prova 154


sottodomini meccanici.

sima e potenziale intermedio, che risultano del tutto analoghi a quelli delle

figure 6.29 e 6.30) i tempi di analisi sono riportati in tabella 6.13. Il gua-


Staggered 95378 sec

Decomp. 3livello (4 Sd Elettr.) 69399 sec



(2 sd meccanici, 4 sd elettrici).

dagno percentuale e analogo a quello ottenuto nel primo caso prova. Infine,

la figura 6.36 riporta l’andamento dei tempi di analisi (su intervallo di tem-

po ridotto e pari a 10−6sec) ottenuti aumentando il numero di sottodomini

elettrici secondo una partizione analoga al caso della trave a doppio incastro

(2× 3, 2× 4).

6.2 Casi prova 155


sottodomini elettrici e meccanici.

Riassunto dei tempi di analisi

La tabella 6.14 riporta un riassunto dei tempi di analisi ottenuti nel caso

di trave a mensola (ttot = 5 · 10−6sec; ∆t = 10−9sec). Anche in questo

caso i tempi relativi al caso di decomposizione di 2 e 3 livello si riferiscono

rispettivamente ad una suddivisione del dominio meccanico in 2 parti e di

quello elettrico in 4 parti. Il guadagno percentuale e relativo all’algoritmo

staggered.

Algoritmo Tempo di analisi Guadagno %

Staggered 95378 sec -




Tabella 6.14: Riassunto dei tempi analisi ottenuti con i vari algoritmi (trave

a mensola).

6.2 Casi prova 156

6.2.3 Considerazioni finali

I due casi prova analizzati hanno consentito da una parte di convalidare i co-

dici implementati e dall’altra di mettere in evidenza alcuni vantaggi e limiti

dei diversi livelli di decomposizione. In particolare e emerso che l’introdu-

zione della decomposizione di 1 livello comporta, in ogni caso, un netto

abbattimento (20− 30%) del tempo di analisi rispetto al metodo staggered,

garantendo allo stesso tempo che la soluzione ottenuta sia la stessa degli

algoritmi gia convalidati.

I vantaggi derivanti dall’introduzione dell’algoritmo GC e del metodo

FETI, applicati rispettivamente alla parte meccanica ed elettrica, dipendo-

no fortemente dall’onere di calcolo associato alla risoluzione dei problemi

meccanico ed elettrico rispettivamente; in generale si possono avanzare le

seguenti osservazioni:

- qualora sia richiesto un certo livello di precisione nella soluzione mec-

canica, e possibile adottare una mesh molto fitta per il solido e una

mesh altrettanto fitta per il dominio elettrico in corrispondenza del-

l’interfaccia, in modo tale da garantire una buona approssimazione nel

calcolo delle forze elettrostatiche; quest’ultima mesh va pero diradan-

dosi nelle zone piu lontane dal dominio meccanico, riducendo l’onere

computazionale associato al calcolo dei potenziali. E in questo caso

che emergono con piu evidenza le potenzialita della decomposizione

meccanica (si veda il caso notevole al paragrafo 6.3.2);

- viceversa, nel caso in cui sia maggiore l’onere di calcolo associato alla

risoluzione del problema elettrico, i vantaggi della decomposizione sono

piu evidenti dall’introduzione del metodo FETI (caso notevole 6.3.3).

Le considerazioni precedenti sottolineano da una parte la versatilita

dei metodi di decomposizione dei domini, applicati al problema elettro-

meccanico, ma allo stesso tempo la necessita di un loro corretto utilizzo

a seconda del problema che ci si trova ad affrontare (dimensioni dei domini

meccanico ed elettrico, tipo di mesh, . . . ); in caso contrario, l’introduzione

della decomposizione non solo non produrrebbe alcun ulteriore vantaggio ma

addirittura un incremento dell’onere di calcolo.

6.3 Casi notevoli 157

6.3 Casi notevoli

In questa sezione si intende mostrare le potenzialita della decomposizio-

ne dei domini nell’ambito di problemi particolari, confrontando risultati e

tempi di analisi, ricavati con diversi algoritmi fra quelli precedentemente

implementati:

1. trave doppiamente incastrata soggetta ad attuazione sinusoidale: con-

fronto staggered -decomposizione 1 livello;

2. portale semplice con attuazione laterale: confronto staggered -decompo-

sizione 2 livello;

3. trave doppiamente incastrata smorzata: confronto staggered -decompo-

sizione 3 livello;

4. trave doppiamente incastrata con doppia attuazione: confronto stag-

gered -decomposizione 1 livello;

6.3.1 Esempio 1: trave a doppio incastro soggetta ad attua-

zione sinusoidale

Si consideri la trave di figura 6.37 analoga per caratteristiche geometriche e

meccaniche a quella del primo caso prova; si immagini pero in questo caso

che il potenziale imposto non sia costante ma che vari con legge sinusoidale,

come mostrato in figura 6.38. Per il caso in esame si ha ampiezza massima

V0 = 30 Volt e periodo T = 10−5 sec. In tabella 6.15 si riportano i parametri

assunti per l’analisi dinamica. La mesh adottata e quella di figura 6.3. La

Figura 6.37: Trave doppiamente incastrata soggetta a differenza di

potenziale variabile nel tempo.


Figura 6.38: Attuazione sinusoidale.






figura 6.39 mostra l’andamento della forza elettrica nodale applicata al nodo

di mezzeria inferiore della trave nei due casi (staggered e decomposizione 1

livello): l’andamento risulta sovrapposto a conferma dell’equivalenza dei ri-

sultati dei due algoritmi. Si osserva che il periodo di oscillazione della forza

elettrica e diverso da quello del potenziale attuato; infatti la forza elettrica

e proporzionale al quadrato del voltaggio, che nel caso in esame risulta sinu-

soidale:

Fel ∝ sin2(2πt

T

)

(6.1)

Utilizzando una nota formula di bisezione, il quadrato del seno puo essere

riscritto come:

Fel ∝1

2

[

1− cos(4πt

T

)]

(6.2)

Il periodo della funzione che compare nella 6.2 puo essere calcolato come


Coordinate dei punti di picco - forza elettrica

Picchi 1 2 3 4

t [10−5 · sec] 0.25 0.75 1.25 1.75

Fel [µN ] −51.9281 −51.9282 −51.9284 −51.9287

Soluzione statica - forza elettrica

Fel(statica) [µN ] −51.9281

Tabella 6.16: Stima dei valori di picco della forza elettrica.

Coordinate dei punti di picco - freccia massima

Picchi 1 2 3 4

t [10−5 · sec] 0.25 0.75 1.25 1.75

u [10−3 · µm] −1.7276 −1.7396 −1.7520 −1.7635

Soluzione statica - freccia massima

u(statica) [10−3 · µm] −1.7231

Tabella 6.17: Stima dei valori di picco della freccia massima.

segue:

T ∗ =2π

(4πT )=

T

2= 0.5 · 10−5 sec (6.3)

Il periodo della forzante e quindi pari alla meta di T , coerentemente con

quanto si evince dalla figura 6.39.

Si osserva inoltre che il segno della forzante e sempre negativo (e rivolta verso

il basso) e che i valori massimi in corrispondenza dei picchi sono praticamente

coincidenti con quello ottenuto dall’analisi statica (il confronto e riportato

nella tabella 6.16).

Si riporta inoltre l’andamento della freccia massima della trave (spo-

stamento verticale del nodo di mezzeria) determinato coi due diversi algorit-

mi e caratterizzato dallo stesso periodo di oscillazione della forzante elettrica.

Dalla figura 6.40 si evince che entrambi gli algoritmi confrontati restituisco-

no lo stesso andamento. La tabella 6.17 riporta i valori in corrispondenza

dei picchi e il loro confronto con la soluzione statica.

Infine, si riportano nella tabella 6.18 i tempi di calcolo per i due casi

analizzati. Ancora una volta si ha la conferma che l’introduzione della de-

composizione (di 1 livello in questo caso) garantisce un vantaggio in termini

computazionali.


Figura 6.39: Forzante elettrica nel nodo di mezzeria della trave.

Figura 6.40: Spostamento verticale del nodo di mezzeria della trave.


Staggered 26386 sec




1livello (attuazione sinusoidale).


6.3.2 Esempio 2: portale semplice con attuazione laterale

L’esempio presentato nel paragrafo seguente riguarda il portale di figura

6.41, soggetto a differenza di potenziale imposta lungo il montante sinistro.

Con questo esempio si vogliono evidenziare i vantaggi computazionali che

la decomposizione meccanica di 2 livello comporta nel caso in cui si adotti

una mesh molto fitta in corrispondenza del dominio meccanico e piu rada

in quello elettrico: si mostrera come, aumentando il numero di sottodomini,

si ottenga un abbattimento dei tempi di calcolo, mantenendo inalterata la

qualita della soluzione. Le caratterstiche dei materiali sono riportate in

tabella 6.1.

Figura 6.41: Geometria del portale.

Analogamente agli esempi precedenti, viene applicata una differenza

di potenziale V fra un elettrodo fisso ed uno deformabile; quest’ultimo e rap-

presentato proprio dal montante sinistro del portale, il cui contorno esterno,

rivolto verso l’elettrodo fisso, identifica l’interfaccia sulla quale avviene lo

scambio di informazioni fra dominio meccanico ed elettrico (forze elettriche

applicate al portale e cedimenti imposti al dominio elettrico). Il voltaggio

imposto V = 500V e costante nel corso di tutta l’analisi dinamica (anda-


mento secondo la funzione di Heaviside) ed e assente lo smorzamento: si

studiano quindi le oscillazioni forzate della struttura.

L’esempio considerato puo essere assimilato approssimativamente al

caso di portale sollecitato da un carico distribuito non uniformemente lun-

go il montante di sinistra e variabile nel tempo; si ottiene quindi una de-

formata analoga qualitativamente a quella che si ricaverebbe in quel caso,

rappresentata in figura 6.42.

Figura 6.42: Deformata del portale sollecitato a carico distribuito variabile

sul montante di sinistra.

Le analisi sono state eseguite con i seguenti codici

- codice 0: staggered ;

- codice 2: decomposizione 2 livello con

– 2 sottodomini meccanici;




Per le analisi che seguono e stata utilizzata una mesh di elementi trian-

golari CST, realizzata con GiD, e suddivisa, di volta in volta, in un diverso

numero di sottodomini, a seconda dell’analisi effettuata (figura 6.43). Nelle

tabelle 6.19, 6.20 e 6.21 vengono indicati il numero di nodi ed elementi che

compongono i sottodomini nei vari casi.

E stata utilizzata una mesh variabile, molto fitta nel dominio mecca-

nico (17 elementi nello spessore) e via via piu rada nel dominio elettrico.


Figura 6.43: Mesh e suddivisione in sottodomini delle casistiche

implementate.


Algoritmo n nodi

staggared 2472

2 sottodomini 1242 1240

3 sottodomini 931 631 930

4 sottodomini 931 321 320 930

6 sottodomini 470 471 321 320 469 471

Tabella 6.19: Numero di nodi dei sottodomini meccanici.

Algoritmo n elementi

staggared 4464

2 sottodomini 2234 4230

3 sottodomini 1670 1126 1668

4 sottodomini 1670 564 562 1668

6 sottodomini 834 836 564 562 832 836

Tabella 6.20: Numero di elementi dei sottodomini meccanici.

La scelta scaturisce da una duplice necessita: mettere in evidenza i vantaggi

computazionali derivanti dalla decomposizione del dominio meccanico, evi-

tando pero di dover eseguire analisi che richiedano tempi non compatibili con

il lavoro svolto. Nella tabella 6.22 sono riportati i parametri di Newmark

dell’analisi dinamica.

I grafici 6.44 e 6.45 mostrano gli andamenti degli spostamenti e poten-

ziali nel nodo A (figura 6.43), valutati con i due diversi algoritmi e nei casi di

diverso numero di sottodomini. Si osservi come le curve siano perfettamente

sovrapposte alla soluzione staggered di riferimento.

La correttezza della soluzione ottenuta puo essere dimostrata anche

atttraverso il confronto di altre grandezze in gioco, quali il regime tensio-

nale nella parte meccanica oppure il campo elettrico nella parte elettrica.

Nella tabella 6.23 si riportano i confronti tra i valori di output delle varie

analisi, valutati nel punto A. Essi si riferiscono all’istante t = 2 · 10−7 sec, in

prossimita del primo picco di spostamento.

n. nodi n. elementi

dominio elettrico 1 279 418

Tabella 6.21: Numero di nodi ed elementi del dominio elettrico.







Figura 6.44: Confronto degli spostamenti in direzione x del nodo A.

Figura 6.45: Confronto dei potenziali nel nodo B.


Grandezza Algoritmo

staggered 2 sd 3 sd 4 sd 6 sd

ux,max [µm] 0 0 0 0 0

ux,min [µm] −0.012485 −0.012485 −0.012485 −0.012485 −0.012485

uy,max [µm] 0.000888 0.000888 0.000888 0.000888 0.000888

uy,min [µm] −0.000998 −0.000998 −0.000998 −0.000998 −0.000998

σx,max [N/m2] 3122388 3122388 3122388 3123715 3123780

σx,min [N/m2] −2682908 −2682908 −2682909 −2682211 −2682174

σy,max [N/m2] 8357410 8357410 8357411 8357897 8357913

σy,min [N/m2] −8880781 −8880781 −8880781 −8881786 −8881772

τxy,max [N/m2] 997531 997530 997530 997588 997586

τxy,min [N/m2] −2166801 −2166801 −2166801 −2167204 −2167199

Ex,max [V/m] 0 0 0 0 0

Ex,min [V/m] −1 · 108 −1 · 108 −1 · 108 −1 · 108 −1 · 108

Ey,max [V/m] 4457 4457 4457 4457 4457

Ey,min [V/m] −30515 −30515 −30515 −30515 −30515

Tabella 6.23: Valori massimi e minimi delle varie grandezze, valutate nel

nodo A, nell’istante t = 2 · 10−7 sec.

Si riportano, a titolo di esempio, i grafici delle grandezze meccaniche

ed elettriche in gioco, valutate con il codice staggered e con la decomposizione

di 2 livello, nel caso di dominio meccanico diviso in sei sottodomini.

Un’ulteriore verifica della correttezza del codice di decomposizione ap-

plicato alla parte meccanica e la continuita di spostamento in corrisponden-

za dei nodi di interfaccia. I grafici 6.54-6.57 verificano tale condizione nelle

diverse suddivisioni.

E, infine, interessante confrontare i tempi di calcolo registrati per le

varie analisi (in tabella 6.24). Anche l’esempio del portale dimostra come,

a parita di condizioni di analisi, l’algoritmo di decomposizione di 1 livello

presenta un vantaggio computazionale evidente nei confronti della procedura

staggered, garantendo in tutti i casi la stessa soluzione.


Figura 6.46: Confronto degli spostamenti in direzione x (in µm).


Figura 6.47: Confronto degli spostamenti in direzione y (in µm).


Figura 6.48: Confronto degli sforzi σx (in N/m2).


Figura 6.49: Confronto degli sforzi σy (in N/m2).


Figura 6.50: Confronto degli sforzi τxy (in N/m2).


Figura 6.51: Confronto del campo elettrico in direzione x (in V/m).


Figura 6.52: Confronto del campo elettrico in direzione y (in V/m).


Figura 6.53: Confronto del potenziale (in V olt).


Figura 6.54: Spostamenti in direzione x del nodo C (2 sottodomini).

Figura 6.55: Spostamenti in direzione x del nodo E (3 sottodomini).


Figura 6.56: Spostamenti in direzione x del nodo C (4 sottodomini).

Figura 6.57: Spostamenti in direzione x del nodo E (6 sottodomini).


Algoritmo Tempo di analisi Guadagno [-]

Staggered 2498 sec –

2 sottodomini 2121 sec −15.09%




Tabella 6.24: Confronto tra i tempi di analisi e relativo guadagno rispetto

all’algoritmo staggered.

6.3.3 Esempio 3: trave a doppio incastro con smorzamento

Si consideri l’ormai consueto problema della trave a doppio incastro di figura

6.3. Si introduca uno smorzamento di tipo proporzionale, in cui la corrispon-

dente matrice e definita come combinazione lineare delle matrici di massa e

di rigidezza:

[C] = αR[M] + βR[K] (6.4)

in cui αR e βR rappresentano i cosiddetti coefficienti di Rayleigh. Questi ul-

timi sono legati ai rapporti di smorzamento relativi al critico ζ1 e ζ2 associati

ai primi due modi di vibrare della trave attraverso le seguenti equazioni (si

veda Abaqus, Theory Manual):

ζ1 =1

2ω1αR + ω1

2 βR

ζ2 =1

2ω2αR + ω2

2 βR

(6.5)

dove ω1 e ω2 sono le pulsazioni corrispondenti ai primi due modi. Ipotiz-

zando uno stesso rapporto di smorzamento (ζ = ζ1 = ζ2), si esplicitano i

coefficienti di Raylegh dalle 6.5:

αR = ζ 2ω1ω2ω1+ω2

βR = ζ 2ω1+ω2

(6.6)

Dall’analisi modale svolta nel paragrafo 6.2.1 si ricavano le due pulsazioni

ω1 e ω2 (figura 6.58). I due coefficienti quindi risultano pari a:


Figura 6.58: Pulsazioni dei primi due modi della trave a doppio incastro.

αR = 0.8897µsec−1

βR = 0.8792 · 10−2 µsec(6.7)

avendo assunto un rapporto di smorzamento ζ = 10%.

Si vuole valutare la risposta della struttura che inizialmente si trova in

quiete (spostamento, velocita e accelerazione nulle), essendo essa soggetta ad

un potenziale costante durante tutta l’analisi (V = 30 Volt); questa situazio-

ne e del tutto analoga a quella delle vibrazioni forzate del problema meccani-

co. Le figure 6.59 e 6.60 riportano rispettivamente l’andamento della freccia

massima e della corrispondente forza elettrostatica nodale in funzione del

tempo, nel caso dei due algoritmi staggered e decomposizione 3 livello (per

quest’ultimo caso e stata adottata una suddivisione in 4 sottodomini elettrici,

analogamente ai casi prova). Si osserva che, come prevedibile, le oscillazioni

tendono a smorzarsi dopo un transitorio iniziale e la soluzione dinamica a

stabilizzarsi nell’intorno di quella statica (ustatica = −1.72309 · 10−3 µsec).

Il periodo di oscillazione e ovviamente lo stesso per la forzante. La tabella

6.25 riporta le stime grafiche dei periodi di oscillazione, da confrontarsi con

il valore teorico (primo periodo modale: T = 1.036µsec).

Infine, si riportano i tempi di calcolo registrati nelle due analisi rispettiva-

mente con algoritmo staggered e decomposizione 3 livello nella tabella 6.26.


Figura 6.59: Andamento della freccia massima nel tempo.

Figura 6.60: Andamento della forza elettrostatica nodale di mezzeria nel

tempo.



Picchi 1 2 3 4 5

t [10−6 · sec] 0.499 1.493 2.487 3.482 4.477

u [10−3 · µm] −2.9931 −2.4061 −2.0887 −1.9176 −1.8255


Intervalli 1-2 2-3 3-4 4-5

Periodo [10−6 · sec] 0.994 0.994 0.995 0.995

Tabella 6.25: Periodo di oscillazione ottenuto numericamente.


Staggered 98898 sec

Decomp. 3livello (4 sd elettrici) 69258 sec



3livello (trave smorzata).

6.3.4 Esempio 4: trave a doppio incastro con doppia attua-

zione

Si consideri il problema di figura 6.61, dato da una trave doppiamente inca-

strata soggetta a due voltaggi imposti V1(t) e V2(t), in generale funzioni del

tempo; sono quindi presenti due domini elettrici separati, indicati rispetti-

vamente con 1 (sottostante) e 2 (sovrastante). Si vuole studiare la risposta

della trave soggetta a doppia attuazione, con voltaggio rispettivamente si-

nusoidale e costante (in Volt):

V1(t) = 30 sin(2πtTV

)

(TV = 10−5 sec)

V2(t) = 30(6.8)

Si riportano gli andamenti in figura 6.62.3 Considerando separatamente le

due attuazioni e tenendo conto dei risultati ottenuti nei casi analoghi dei

paragrafi 6.2.1 e 6.3.1 (primo caso prova e attuazione sinusoidale), si puo

tentare una previsione sull’andamento della risposta, espressa ancora una

volta in termini di freccia massima. Si fanno le seguenti considerazioni:

3La scelta non e casuale, in quanto i casi di attuazione singola con V1 e V2 sono stati in

precedenza studiati separatamente; sara quindi possibile stabilire un confronto coi risultati

ottenuti in quelle occasioni.


Figura 6.61: Trave a doppio incastro soggetta a doppia attuazione.

- effetto di V1: l’attuazione sinusoidale, caratterizzata da un periodo

pari a TV = 10−5 sec, provoca delle oscillazioni che, come gia osser-

vato al paragrafo 6.3.1, hanno periodo dimezzato rispetto a quello del

voltaggio;

- effetto di V2: l’attuazione costante produce oscillazioni secondo il perio-

do proprio della trave, assunto coincidente con quello del modo prin-

cipale (primo). Dall’analisi modale effettuata al paragrafo 6.2.1 era

emerso un primo periodo modale pari a T = 1.036 · 10−6 sec;

- i due effetti devono essere combinati.

Assumendo quindi come parametro di risposta la freccia massima, si ripor-

ta in figura 6.63 il suo andamento nel tempo, calcolato coi due algoritmi

che si e deciso di utilizzare per questo esempio: staggered e decomposizione

1livello. Il risultato conferma quanto previsto: l’attuazione costante supe-

riore provoca delle oscillazioni aventi periodo prossimo a quello della trave

(circa 10−6sec), analoghe a quelle di figura 6.6 del primo caso prova; queste

pero risentono dell’attuazione sinusoidale inferiore, per cui subiscono una

traslazione in direzione y la cui entita varia essa stessa nel tempo secon-

do un andamento sinusoidale. In figura 6.63 la linea tratteggiata in verde

rappresenta proprio il valore medio, non costante ma funzione sinusoidale,

attorno a cui avvengono le oscillazioni. Infine, il confronto fra i tempi di

calcolo (tabella 6.27) evidenzia ancora una volta il vantaggio derivante dalla


Figura 6.62: Andamento nel tempo dei due potenziali di attuazione.


Staggered 2987 sec




1livello (doppia attuazione).

decomposizione dei domini.


Figura 6.63: Andamento della freccia massima nel tempo.

Capitolo 7

Conclusioni

Nel presente lavoro di tesi sono state applicate varie tecniche di decomposi-

zione dei domini al problema elettro-meccanico accoppiato; questo appartie-

ne alla categoria dei problemi a piu campi, nei quali sono presenti fenomeni

che possono interessare piu ambiti della fisica e che interagiscono fra di loro.

Il problema elettro-meccanico e tipico dei micro-sistemi elettro-meccanici

(MEMS): le dimensioni di tali dispositivi sono tali (ordine di grandezza del

micrometro) che le forze gravitazionali risultano trascurabili nei confronti di

forze elettriche e adesive. Le simulazioni numeriche che vengono effettuate

per determinare le forze elettrostatiche, che e necessario determinare con

sufficiente precisione soprattutto in ambito progettuale, comportano la riso-

luzione di un sistema di equazioni non solo accoppiato, ma anche non-lineare,

che puo richiedere un consistente dispendio di risorse computazionali; il la-

voro svolto ha voluto dimostrare che l’introduzione della decomposizione dei

domini (a piu livelli) puo comportare una notevole riduzione del tempo di

calcolo, sfruttando gli algoritmi originariamente proposti per problemi pu-

ramente meccanici ed estesi al particolare problema oggetto di studio.

Le analisi volte a dimostrare l’effettivo vantaggio che le tecniche di decom-

posizione possono comportare sono state effettuate con l’ausilio di Matlab,

in cui sono stati implementati ex novo dei codici di calcolo associati ai vari

algoritmi proposti nel capitolo 5; questi sono stati utilizzati per analizzare

dei casi prova la cui soluzione e stata confrontata con quella fornita dalla

procedura staggered e dal codice commerciale Comsol. Accanto alla valida-

zione dei nuovi algoritmi presentati, particolare attenzione e stata dedicata

al confronto dei tempi di calcolo dei vari algoritmi, soprattutto con riferi-

mento al metodo staggered. I risultati ottenuti hanno messo in evidenza

184

185

che l’introduzione del 1 livello di decomposizione, che sfrutta la naturale

suddivisione fra dominio meccanico ed elettrico e che comporta, nello spi-

rito dell’algoritmo di Gravouil e Combescure ([8]), una scomposizione delle

soluzioni - spostamenti e potenziali - in due contributi free e link, garantisce

un risparmio computazionale che si attesta intorno al 20-25% del tempo di

calcolo necessario per un’analisi staggered.

Per quanto riguarda la decomposizione della parte meccanica (2 livello) e

di quella elettrica (3 livello), i casi prova hanno mostrato la possibilita di

raggiungere un abbattimento dei tempi di analisi anche del 30%. Cio signi-

fica che, rispetto alla procedura monolitica, il vantaggio computazionale si

eleva al 60%.

Lo studio svolto evidenzia da una parte i benefici dei metodi di decomposizio-

ne dei domini applicati al problema elettro-meccanico, dall’altra la necessita

di un loro corretto utilizzo a seconda del problema che ci si trova ad af-

frontare. Infatti i vantaggi derivanti dall’introduzione della decomposizione

meccanica (2 livello) sono piu evidenti nel caso in cui la maggior parte delle

risorse computazionali sia spesa per la risoluzione del problema meccanico:

e il caso del secondo caso notevole (portale attuato lateralmente), nel quale

la mesh meccanica risulta molto piu fitta di quella elettrica, che invece di-

venta piu rada allontanandosi dall’interfaccia con la parte meccanica. Allo

stesso modo, la decomposizione di 3 livello (FETI, [1]) applicata al caso di

mesh elettrica rada non produce vantaggi apprezzabili. Si deduce quindi la

forte dipendenza del risparmio dal tipo di mesh utilizzata; a seconda della

modellazione adottata sara necessario scegliere il livello di decomposizione

piu adatto, in grado di garantire il massimo guadagno possibile.

L’adozione delle tecniche di decomposizione dei domini per la risolu-

zione del problema elettro-meccanico ricopre un interesse pratico soprattutto

a livello progettuale, in quegli ambiti dove e richiesto effettuare un gran nu-

mero di simulazioni: la diminuzione dei tempi di calcolo, sfruttando il giu-

sto livello di decomposizione, puo portare ad un incremento considerevole

dell’efficienza.

Nella progettazione dei micro-sistemi e di grande importanza riuscire

a determinare la cosiddetta curva di pull-in, dalla quale e possibile dedur-

re il massimo voltaggio applicabile al singolo dispositivo senza che i piatti

del condensatore che lo costituiscono collassino l’uno sull’altro. Per ottenere

una migliore approssimazione del voltaggio di pull-in, in corrispondenza del

quale si ha quel fenomeno di instabilita, e necessario eliminare l’ipotesi, im-

186

plicitamente assunta a monte del lavoro, di piccoli spostamenti della parte

meccanica. Cio comporta la necessita di tener conto non solo della defor-

mazione del dominio elettrico, ma anche di quello meccanico, aggiornandone

le coordinate nodali ad ogni passo temporale. Cio comporta un aggravio

in termini computazionali, dovendo aggiornare istante per istante la matri-

ce di rigidezza e delle masse del dominio meccanico, dipendenti dalla sua

configurazione (variabile).

I risultati ottenuti nel presente lavoro di tesi valgono per il proble-

ma elettro-meccanico; questo appartiene, come accennato, alla piu ampia

categoria dei problemi multi fisici ed e caratterizzato da fisiche non sovrap-

poste (contrariamente, per esempio, al problema termo-meccanico nel quale

dominio elastico e termico sono, in genere, coincidenti). Gli algoritmi propo-

sti in questo lavoro possono essere estesi, con opportuni arrangiamenti che

dipendono dalle fisiche in gioco, a tutti quei problemi multi fisici caratteriz-

zati da domini non sovrapposti e governati da equazioni analoghe a quelle

del problema elettro-meccanico. L’interazione fluido-struttura ne e il tipico

esempio.

Appendice A

Requisiti matematici degli

algoritmi di decomposizione

dei domini

A.1 Metodo del gradiente coniugato

Il metodo del gradiente coniugato permette di risolvere, attraverso una pro-

cedura iterativa, un sistema di equazioni nella forma:

Ax = b (A.1)

dove A e la matrice dei coefficienti, per ipotesi simmetrica e definita positiva,

x il vettore delle incognite e b il vettore dei termini noti.

Il metodo del gradiente coniugato si basa sull’ipotesi che la soluzione

del sistema (A.1) renda minima la forma quadratica:

Q(x) =1

2xTAx− bx (A.2)

Si procede quindi in modo iterativo per trovare il punto di minimo sulla

superficie Q(x). Si costruisce una successione di approssimazioni xk, con-

vergente alla soluzione esatta, a partire da un vettore iniziale arbitrario x0,

secondo la formula ricorsiva:

xk+1 = xk + αkpk (A.3)

187

A.1 Metodo del gradiente coniugato 188

dove pk identifica la direzione che consente al vettore dato dalla (A.3) di

minimizzare la forma quadratica Q(x) lungo la direzione stessa; i vettori pk

sono ottenuti mediante A-ortogonalizzazione dei residui. 1 Si definisce il

coefficiente βk tale da soddisfare:

pk+1 = rk+1 + βkpk (A.4)

dove rk = b −Axk e il residuo associato all’iterazione k. Il nuovo residuo,

relativo all’iterazione k+1, si ricava facilmente utilizzando la relazione (A.3).

Imponendo la condizione di A-ortogonalita fra i vettori pk e pk+1:

pTkApk+1 = pT

kA(rk+1 + βkpk) = 0 (A.5)

si ottiene l’espressione del coefficiente βk:

βk = −pTkArk+1

pTkApk

(A.6)

Sostituendo la (A.3) all’interno della forma quadratica (A.2), si ottiene:

Q(x) =1

2(xk + αkpk)

TA(xk + αkpk)− b(xk + αkpk) (A.7)

Si cerca quindi il punto di minimo della (A.7), imponendo che sia nulla la

derivata prima:

∂Q(x)

∂αk= pT

k (Axk − b) + αkpTkApk = 0 (A.8)

Esplicitando l’equazione (A.8) rispetto al coefficiente αk:

αk =pTk rk

pTkApk

(A.9)

L’algoritmo del gradiente coniugato puo essere sintetizzato nei seguenti

passi:

1Due vettori xi e xj si dicono A-ortogonali se xTi Axj = 0 per i 6= j.

A.2 Gradiente coniugato precondizionato (PGC) 189

1. inizializzazione:

dato x0

r0 = b−Ax0

p0 = r0

(A.10)

2. utilizzando la (A.6) e la (A.9) si calcolano, in corrispondenza della k-

esima iterazione (con k = 0, 1, 2, 3, ..., n):

xk+1 = xk + αkpk

rk+1 = rk − αkApk

(A.11)

3. controllo di convergenza:

‖rk+1 − rk‖‖rk‖

< ε, (A.12)

dove ε e una tolleranza prefissata. In corrispondenza dell’iterazione in

cui il test ha esito positivo, si assume che l’ultima approssimazione del

vettore x sia soluzione del sistema (A.1).

Il metodo del gradiente coniugato trova particolare applicazione nel ca-

so di matrici dei coefficienti sparse, per le quali non e conveniente adottare

una triangolazione della matrice A; infatti in questo caso le matrici trian-

golari che si ottengono hanno un numero di elementi non nulli superiore

rispetto alla matrice originaria.

A.2 Gradiente coniugato precondizionato (PGC)

Il precondizionamento e una tecnica utilizzata per aumentare il numero di

condizionamento di una matrice. Si introduca il precondizionatore M, ma-

trice simmetrica e definita positiva, definito come:

M = ETE (A.13)


dove E e una matrice non singolare. Inoltre deve essere garantita la condi-

zione:

K(E−TAE−1) ≪ K(A) (A.14)

Il metodo del gradiente coniugato precondizionato e semplicemente il meto-

do GC (A.1) applicato al sistema:

E−TAE−1(Ex) = E−Tb (A.15)

L’algoritmo PGC puo essere sintetizzato secondo i passi seguenti:

1. inizializzazione

dato x0

r0 = b−Ax0

(A.16)

2. in corrispondenza della k-esima iterazione (con k = 0, 1, 2, 3, ..., n), si

eseguono in sequenza le seguenti operazioni:

zk = M−1rk (A.17)

βk+1 = rTk+1zk (A.18)

pk+1 = zk +βk+1

βkpkzk (A.19)

αk+1 =βk+1

pTk+1Apk+1

(A.20)

xk+1 = xk + αk+1pk+1 (A.21)


rk+1 = rk − αk+1Apk+1 (A.22)

3. controllo di convergenza:

‖rk+1 − rk‖‖rk‖

< ε (A.23)

Appendice B

Integrazione delle equazioni

del moto

Lo schema di integrazione temporale di Newmark appartiene alla categoria

dei metodi passo-passo; questi consentono di ricavare una soluzione dalle

equazioni di equilibrio dinamico (equazioni differenziali ordinarie nella va-

riabile tempo), che difficilmente sono risolvibili in forma chiusa. Essi preve-

dono l’introduzione di una discretizzazione temporale, con la quale il tempo

di riferimento viene suddiviso in intervalli finiti ∆t; l’equilibrio dinamico e

soddisfatto negli istanti che delimitano ciascuno di quegli intervallini (nodi

della mesh temporale). Al generico istante tk, l’equilibrio del sistema puo

essere scritte come:

Muk + f(uk,uk) = g(tk) (B.1)

Si concentri l’analisi nell’intervallo temporale ∆tk = tk+1 − tk, dove si sup-

pone siano note all’istante tk la configurazione del sistema, le velocita e le

accelerazioni. Si scrive quindi l’equazione di equilibrio dinamico in forma

incrementale sul passo:

M∆uk +∆fk(∆uk,∆uk) = ∆gk(tk) (B.2)

dove le quantita incrementali sono definite come:

192

193

∆uk = uk+1 − uk

∆fk = fk+1 − fk

∆gk = gk+1 − gk

(B.3)

Il vettore delle forze resistenti f e funzione solamente dell’atto di moto del

sistema; il suo incremento puo essere espresso come combinazione lineare

degli incrementi di spostamento e velocita:

∆fk = Kseck ∆uk +Csec

k ∆uk (B.4)

Nella formulazione (B.4), le matrici di rigidezza Kseck e smorzamento Csec

k

sono tipicamente incognite, in quanto dipendono dalla soluzione in termini di

spostamenti, velocita e accelerazioni in corrispondenza degli istanti discreti.

Sostituendo la (B.4) nel problema (B.2) si ottiene:

M∆uk +Kseck ∆uk +Csec

k ∆uk = ∆gk. (B.5)

Nel caso di comportamento lineare, le matrici di rigidezza e smorzamento

sono note, in quanto costanti nell’intervallo di analisi. Nel caso piu gene-

rale e necessario riscrivere la (B.5) in forma approssimata, sostituendo alle

matrici secanti le corrispondenti matrici tangenti, calcolate ad inizio passo

come sviluppo in serie di taylor arrestato al primo ordine:

Ctank = ∂f(u,u)

∂u

∣∣∣tk

Ktank = ∂f(u,u)

∂u

∣∣∣tk

(B.6)

Si ottiene l’equazione di riferimento per l’integrazione temporale con metodi

passo-passo:

M∆uk +Ktank ∆uk +Ctan

k ∆uk = ∆gk (B.7)

Il problema B.7 e composto da n equazioni nelle 3n incognite uk+1, uk+1,uk+1.

Affinche il problema sia ben posto e necessario considerare ulteriori 2n equa-

zioni, date dalle relazioni:

194

u(τ) = uk +∆tk

∫ τ

0u(τ)dτ (B.8)

u(τ) = uk +∆tkuk + (∆tk)2

∫ τ

0

(∫ τ

0u(τ)dτ

)

dτ (B.9)

dove

τ =(tk+1 − tk)

∆t(B.10)

e la coordinata temporale adimensionalizzata. Integrando per parti le rela-

zioni (B.9), si ottiene:

u(τ) = uk +∆tkuk + (∆tk)2

∫ τ

0(1− τ)u(τ) dτ (B.11)

Se fosse noto l’andamento delle accelerazioni nell’intervallo di anali-

si, sarebbe possibile risalire, attraverso le (B.11), ai valori di spostamento

e velocita in un generico istante appartenente all’intervallo. Per questa ra-

gione, i metodi di integrazione passo-passo introducono un’approssimazione

sull’andamento delle accelerazioni nell’intervallo ∆tk.

Un possibile modello per l’accelerazione puo essere definito come segue:

u(τ) = ωA(τ)ρA (B.12)

dove ωA e la matrice delle funzioni di interpolazione che dipende solamente

da τ ; ρA e il vettore contenente i valori di accelerazione nei punti interpolati.

Gli algoritmi di integrazione passo-passo si differenziano tra loro a se-

conda del modello assunto per l’accelerazione, definito sul generico passo di

calcolo. I metodi della famiglia di Newmark assumono:

ωA(τ) = [ω1(τ)− ω2(τ)]I2n (B.13)

ρA =

[

qk

qk+1

]

(B.14)

195

dove le relazioni (B.13) e (B.14) mostrano che l’andamento dell’accelerazione

sul singolo ∆t dipenda solamente dai valori ad inizio e fine passo; inoltre si

scelgono funzioni interpolanti uguali per ogni passo temporale e per ogni

coordinata libera del sistema.

E possibile stabilire un’analogia tra le funzioni di forma della discre-

tizzazione temporale e quelle utilizzate nel metodo degli elementi finiti;

con questa formulazione le ωi devono essere tali da garantire la continuita

dell’accelerazione nel passo.

In generale non e necessario definire espressamente le funzioni ω1 e ω2;

si introducono due parametri, β e γ, legati ad esse dalle relazioni:

∫ 1

0(1− τ)ω1(τ) dτ =

1

2− β (B.15)

∫ 1

0ω1 (τ)dτ = 1− γ (B.16)

β e γ sono due parametri liberi, i cui valori caratterizzano i diversi metodi

di integrazione di Newmark.

Introducendo le relazioni (B.15),(B.16) e l’espressione approssimata

dell’accelerazione sul passo, valutate per τ = 1, nel problema composto dal-

le (B.7), (B.8) e (B.9), si ottiene:

M∆uk +Ktank ∆uk +Ctan

k ∆uk = ∆gk

∆uk = ∆tkuk +∆tk γ∆uk

∆uk = ∆tkuk +12(∆tk)

2uk + (∆tk)2β∆uk

(B.17)

Il secondo e il terzo blocco di equazioni del sistema (B.17) forniscono l’incre-

mento di velocita (∆u) e di spostamento (∆u) in funzione di grandezze note,

quali le velocita (uk) e accelerazioni (uk) ad inizio passo, ed in funzione di

una grandezza incognita, l’incremento di accelerazione ∆u che le coordinate

libere della struttura subiscono durante il passo.

Le procedure di tipo passo-passo consentono pertanto, in ogni passo

temporale, di scrivere un sistema di equazioni algebriche del tipo:

D

[

∆uk

∆uk

]

= bk (B.18)

B.1 Algoritmi impliciti 196

La matrice dei coefficienti del sistema (B.18) dipende dalle caratteristiche

della struttura e dall’algoritmo di integrazione adottato. Il vettore dei ter-

mini noti bk dipende dai carichi dinamici applicati e dalle condizioni iniziali

del passo.

Possono presentarsi due diverse situazioni, che consentono di distin-

guere tra procedure esplicite e procedure implicite. Nelle procedure

esplicite la matrice dei coefficienti D e diagonale, ed e pertanto possibile la

determinazione immediata dell’incremento di velocita e accelerazione. Vice-

versa, nelle procedure implicite, la matrice D non e diagonale e per trovare la

soluzione del problema (B.18) e necessario risolvere il sistema di equazioni.

In generale, le procedure esplicite richiedono un basso onere compu-

tazionale per ottenere la soluzione, ma comportano problemi di precisione

e stabilita; percio esse richiedono passi di integrazione molto piccoli. Nelle

procedure implicite, al contrario, l’onere computazionale e decisamente piu

elevato ma le procedure sono solitamente piu stabili e accurate.

Le considerazioni precedenti dimostrano come non sia possibile stabili-

re a priori quale sia la migliore tipologia da adottare nelle analisi dinamiche;

la scelta dipende dal particolare tipo di problema.

Analizzando il problema (B.17), si osserva come i parametri β e γ

hanno l’effetto di pesare il contributo dato dall’incremento di accelerazione,

fra gli istanti di inizio e fine passo, rispettivamente alla variazione di velocita

e a quella di spostamento; tale incremento deve essere sommato ai termini

che si otterrebbero se, nel passo ∆tk, le accelerazioni si mantenessero costanti

e pari al valore di inizio passo.

B.1 Algoritmi impliciti

Gli algoritmi impliciti sono caratterizzati da β = 0: lo spostamento e la velo-

cita all’istante finale del passo dipendono dall’incremento dell’accelerazione

nel passo stesso:

∆uk = ∆tkuk +1

2(∆tk)

2uk + (∆tk)2β∆uk (B.19)

∆uk = ∆tkuk + γ∆t∆uk (B.20)


Sostituendo la (B.20) nella relazione (B.19), si ottiene:

∆uk =1

β∆t2kuk −

1

β∆tk, uk −

1

2βuk (B.21)

e inserendo le accelerazioni (B.21) nelle (B.20):

∆uk = ∆tkuk + γ1

β∆tkuk −

γ

βuk −

γ∆tk2β

uk. (B.22)

Le relazioni degli incrementi di accelerazione (B.22) e velocita (B.21), in

funzione dell’incremento di spostamento e delle grandezze ad inizio passo,

se sostituite all’interno dell’equazione di equilibrio dinamico, permettono di

giungere alla seguente formulazione:

[

M 1β∆t2

k

+C γβ∆tk

+K]

∆uk = ∆gk +M[

1β∆tk

uk +12β uk

]

+

+C[

∆tk

(γ2β − 1

)

uk +γβ uk

](B.23)

che scritta in forma sintetica risulta:

M∆uk = ∆Fk (B.24)

dove

M =

[

M1

β∆t2k+C

γ

β∆tk+K

]

, (B.25)

e la matrice delle masse, mentre

∆F = ∆gk+M

[1

β∆tkuk +

1

2βuk

]

+C

[

∆tk

(γ

2β− 1

)

uk +γ

βuk

]

(B.26)

e il vettore dei carichi esterni equivalenti.

B.1.1 Metodo di Newmark dell’accelerazione costante

Nelle analisi svolte in questo lavoro di tesi, si e utilizzato uno schema di

integrazione di Newmark con accelerazione costante nel passo temporale. In


questo caso, i valori dei parametri β e γ sono:

β = 14 γ = 1

2(B.27)

Come dice il nome stesso, il metodo assume che le accelerazioni si man-

tengano costanti nel generico ∆tk e sia pari alla media tra i valori ad inizio

e fine passo. Cio implica che le condizioni di continuita dell’accelerazione in

corrispondenza degli estremi del passo non siano soddisfatte. Tale metodo

risulta essere incondizionatamente stabile

B.1.2 Metodo di Newmark dell’accelerazione lineare

Il metodo di Newmark dell’accelerazione lineare assume che i parametri val-

gano:

β = 16 γ = 1

2(B.28)

Il metodo assume che l’accelerazione vari linearmente tra i valori iniziale

e finale del passo ∆tk. Anche questo algoritmo risulta condizionatamente

stabile.

Appendice C

Caratteristiche dei processori

Nella tabella C.1 vengono elencate le caratteristiche dei processori utiliz-

zati per le analisi e di conseguneza per il confronto dei tempi tra i diversi

algoritmi.

Processore Sistema Operativo

a Intel(R) Core(TM) 2 Duo

CPU E4600 @ 2.40GHz

2.39GHz

32 bit

b Intel(R) Core(TM) i5 CPU

M460 @ 2.53GHz 2.53GHz

64 bit

Tabella C.1: Caratteristiche dei processori.

199

Appendice D

Descrizione dei codici

Gli algoritmi di decomposizione dei domini applicati alla soluzione del pro-

blema elettro-meccanico, presentati nel capitolo 5, sono stati implementati

ex-novo in codici ad elementi finiti 2D in dinamica, realizzati con Matlab.

E stato adottato uno schema di integrazione implicito, con accelerazione co-

stante sul singolo passo temporale, della famiglia di Newmark (si veda B).

Nei paragrafi successivi si descrive la struttura dei codici implementati.

D.1 Decomposizione di 1 livello

Il codice di decomposizione di 1 livello risolve il problema elettro-meccanico

considerando le soluzioni associate ai due domini meccanico ed elettrico (spo-

stamenti e potenziali rispettivamente) come somma di due contributi free e

link.

Il codice e costituito da un blocco principale (main decomp 1livello),

nel quale si susseguono le fasi di input, di assemblaggio, le operazioni previste

dall’algoritmo, post-processing. In realta il codice e organizzato in maniera

tale da richiamare subroutines atte all’esecuzione di determinate operazioni.

Nella fase di input vengono caricati due file .txt, esportati dal pro-

gramma GiD, contenenti rispettivamente le coordinate e le incidenze nodali.

Le mesh degli esempi analizzati sono costituite da elementi CST.

Il codice, attraverso i file importati, identifica e distingue i domini

meccanico ed elettrico, in base a condizioni sulla geometria imposte dall’u-

tente; sia le quantita geometriche (coordinate, incidenze,. . . ) sia le grandezze

derivate (matrici, vettori,. . . ) vengono salvate all’interno di strutture dati,

identificate rispettivamente con il suffisso mech ed elec.

200

D.1 Decomposizione di 1 livello 201

Oltre alle grandezze citate, vengono salvate, per ogni problema fisi-

co, le condizioni di carico, le condizioni iniziali e di vincolo, queste ultime

espresse in termini di imposizione di potenziale per il problema elettrico, e di

spostamento per la parte meccanica. In aggiunta alle condizioni al contorno

elettriche e meccaniche, e necessario definire anche i vincoli del problema

pseudo-meccanico, imposti in termini di spostamenti.

Si definiscono le interfacce tra i due domini meccanico ed elettrico, sal-

vando le duplici numerazioni di nodi e gradi di liberta e le relative incidenze

sottodominio-interfaccia. Inoltre vengono inserite le proprieta dei materiali.

A questo punto grazie all’ausilio di funzioni particolari per assem-

blaggio vengono costruite, all’esterno del ciclo temporale, tutte le quantita

indipendenti dal tempo di analisi, necessarie per la risoluzione del problema.

Si riporta di seguito una breve descrizione delle subroutines utilizzate:

assemb M

Assembla la matrice delle masse di ogni sottodominio meccanico, salvandola

nella struttura dati M. All’interno del ciclo sugli elementi del singolo sottodo-

minio viene richiamata la funzione mas plane t3, che restituisce la matrice

delle masse locale del singolo elemento CST.

assemb K

Assembla la matrice di rigidezza di ogni sottodominio meccanico, salvandola

nella struttura dati Ku. All’interno dei cicli sugli elementi viene richiamata la

funzione rig plane t3, che calcola la matrice di rigidezza locale del singolo

elemento CST.

assemb f

Esegue l’assemblaggio del vettore delle forze nodali equivalenti meccaniche,

dovuta alla presenza di carichi esterni: forze di volume, di superficie o ca-

richi applicati direttamente nei nodi. All’interno del ciclo di assemblaggio

sugli elementi, vengono richiamate le functions bf space t3 e bf plane t3,

che calcolano per il singolo elemento il vettore dei carichi nodali equivalenti

rispettivamente ai carichi distribuiti sul volume (area nel caso 2D) e sulla

superficie (lati).

assemb Kfi0

La subroutine assembla, per ogni sottodominio elettrico, la matrice di ri-


gidezza elettrica associata alla configurazione indeformata, salvandola nella

struttura dati Kfi0. All’interno dei cicli sugli elementi di un sottodominio

viene richiamata la funzione rigelec plane t3, che calcola la matrice di

rigidezza locale del singolo elemento CST.

assemb Kelec

Assembla la matrice di rigidezza del dominio pseudo-meccanico, salvandola

nella struttura Kelec. Le operazioni sono analoghe a quelle eseguite nella

funzione assemb K. All’interno dei cicli sugli elementi viene utilizzata ancora

una volta la funzione rig plane t3.

Le subroutines precedenti consentono l’assemblaggio delle quantita che

rimangono costanti durante l’analisi. All’interno del corpo principale del

main vengono richiamate altre subroutines:

assemb Kfi

Formalmente analoga alla assemb Kfi0, differisce per le quantita fornite in

input, che in questo caso tengono conto della deformazione del dominio elet-

trico (alle coordinate iniziali sono aggiunti gli spostamenti nodali attuali).

La matrice calcolata si riferisce quindi alla situazione deformata e viene sal-

vata all’interno della struttura Kfi, che viene sovrascritta ad ogni iterazione.

emaxwell tensor t3

Richiamata a valle della risoluzione del problema elettrico link, all’interno

del ciclo iterativo, consente di determinare le forze di Coulomb nodali del

dominio elettrico attraverso l’utilizzo del tensore di Maxwell.

A valle del ciclo iterativo si colloca la fase di post-processing : le quantita

ottenute dalla risoluzione del problema elettro-meccanico vengono salvate in

corrispondenza di determinati istanti all’interno di files .vtk, che vengono

poi forniti in input al programma ParaView. I files .vtk vengono costruiti

dalle apposite funzioni, di seguito brevemente descritte insieme alla subrou-

tine per il calcolo degli sforzi:

stress plane t3

Calcola gli sforzi nel singolo elemento CST a partire dalle coordinate e spo-


stamenti nodali dell’elemento. Restituisce per ogni elemento un vettore con-

tenente le tre componenti di sforzo piano.

Writing Ucst t3

Costruisce un file Uplot.vtk, contenente tutte le informazioni, riguardanti

geometria e spostamenti, che devono essere passate a ParaView.

Writing elCST t3

Compila un file Vplot.vtk che, una volta caricato in ambiente ParaView,

consente di rappresentare la distribuzione di potenziale e il valore di campo

elettrico in direzione x e y sul dominio elettrico.

Writing stressCST t3

La funzione riceve in input la matrice degli sforzi di tutta la struttura (co-

stituita dall’insieme dei vettori restituiti da stress plane t3 ), che viene in-

terpolata linearmente sui nodi del dominio meccanico attraverso il comando

pdertni. Anche in questo caso la funzione compila un file stressplot.vtk,

che permette di visualizzare le distribuzioni delle tre componenti di sforzo

piano sul dominio meccanico.

Si riporta uno schema del codice di decomposizione 1 livello in figura D.1.


L’algoritmo di decomposizione 2 livello applica l’algoritmo di Gravouil e

Comberscure ([8]), proposto per problemi di dinamica strutturale, alle due

fasi di risoluzione del problema meccanico free e link (rispettivamente in

blu e in viola nello schema di figura D.2). L’algoritmo di GC e applicato

quindi due volte all’interno del codice, per cui le operazioni ad esso associa-

te sono state organizzate in subroutines che vengono richiamate dal main

(main decomp 2livello).

GC t3

E la funzione che gestisce le operazioni dell’algoritmo GC. Essa richiama al-

tre tre subroutines: GC free t3 per la risoluzione della parte free, GC int t3

per il problema di interfaccia e GC link t3 per la parte link. Si osserva che la

subroutine GC t3 viene richiamata un’unica volta in ogni istante temporale


Figura D.1: Schema del codice di decomposizione 1 livello.


per la risoluzione del problema meccanico free, ad ogni iterazione invece per

quanto riguarda il problema link.

La fase di input del programma di decomposizione di 2 livello rima-

ne pressoche inalterata; infatti il salvataggio delle grandezze all’interno di

strutture dati permette di adattare facilmente il codice alla presenza di sotto-

domini meccanici. La partizione del dominio meccanico comporta il calcolo

degli operatori di condensazione H nella fase di assemblaggio, necessari per

la risoluzione dei problemi di interfaccia fra domini meccanici.

La fase di output e analoga all’algoritmo di decomposizione 1 livello.


La decomposizione di 3 livello introduce, in aggiunta alle decomposizioni

dei livelli precedenti, la decomposizione del dominio elettrico. Al capitolo 5

si mostra che, poiche i problemi elettrico e pseudo-meccanico sono staziona-

ri, e necessario adottare il metodo FETI. In particolare si mette in evidenza

come quel metodo debba essere applicato in modo diverso per i due tipi di

problemi precedenti, trattandosi rispettivamente di un problema a potenzia-

le (elettrico) e di uno meccanico (su dominio elettrico).

Per l’implementazione delle subroutine di decomposizione secondo FETI si

e preso spunto dal codice implementato in [17], le cui suroutines sono state

adattate ai due problemi in esame, elettrico e pseudo-meccanico. Per tener

conto della diversa natura dei due problemi risolti su dominio elettrico, sono

quindi state realizzate due serie di funzioni, una per ciascuno di quei due

problemi. La figura D.3 mostra, all’interno del codice, le operazioni in cui

interviene la decomposizione di 3 livello. Di seguito si riporta una breve

descrizione delle subroutines principali.

Interf 2D elec

Assembla la matrice dei coefficienti e il vettore dei termini noti del problema

di interfaccia associato al problema elettrico.

Interf 2D psm

Assembla la matrice dei coefficienti e il vettore dei termini noti del problema

di interfaccia associato al problema pseudo-meccanico.


FETI 2D elec

Risolve il problema di interfaccia elettrico e, noto il vettore dei moltiplicatori

di Lagrange, risale al vettore soluzione (potenziali nodali) dei sottodomini

e, successivamente, dell’intero dominio.

FETI 2D psm

Risolve il corrispondente problema di interfaccia, determinando il vettore dei

moltiplicatori di Lagrange e risalendo agli spostamenti nodali del dominio

elettrico.

Bibliografia

[1] C. Farhat, F. Roux: A method of finite element tearing and intercon-

necting and its parallel solution algorithm, International Journal for

Numerical Methods in Engineering, Vol. 32, 1205-1227, 1991.

[2] A. Corigliano, A. Taliercio: Meccanica computazionale, soluzione del

problema elastico lineare, Ed. Esculapio, 2005.

[3] N. Mahjoubi, A. Gravouil, A. Combescure: Coupling subdomains

with heterogeneous time integrators and incompatible time steps,

Computational mechanics, Vol. 44, 825-843, 2009.

[4] A. Combescure, A. Gravouil, B. Herry: An algorithm to solve transient

structural non-linear problems for non-matching time-space domains,

Computers and Structures 81, 12111222, 2003.

[5] C. Farhat, P.S. Chen, J. Mandel: A scalable Lagrange multiplier based

domain decomposition method for time-dependent problems, Internatio-

nal Journal for Numerical Methods in Engineering, Vol. 38, 38313853,

1995.

[6] C. Farhat, M. Lesoinne, P. LeTallec, K. Pierson, D. Rixen: FETI-DP: a

dual-primal unified FETI method, International Journal for Numerical

Methods in Engineering, Vol. 50, 15441551, 2001.

[7] A. Combescure, A. Gravouil: Multi-time-step explicit-implicit method

for non-linear structural dynamics, International Journal for Numerical

Methods in Engineering, Vol. 50, 199-225, 2001.

[8] A. Combescure, A. Gravouil: Multi-time-step explicit-implicit method

for non-linear structural dynamics, International Journal for Numerical


209

BIBLIOGRAFIA 210

[9] J. Mandel: Balancing Domain Decomposition, Comm. on Numerical


[10] P. Gosselet, C. Rey: Non-overlapping domain decomposition me-

thods in structural mechanics, Archives of Computational Methods in

Engineering, Vol. 13, 515-572, 2006.

[11] F. Confalonieri: Metodi di decomposizione dei domini per problemi di

dinamica strutturale, Tesi di Laurea Specialistica, Politecnico di Milano,

2009.

[12] N. Fragakis, M. Papadrakakis: The mosaic of high performance domain

decomposition methods for structural mechanics: formulation, interre-

lation and numerical efficiency of primal and dual methods, Comput.

Methods Appl. Mech. Engrg., Vol. 191, 3799-3830, 2003.

[13] N. Fragakis, M. Papadrakakis: The mosaic of high performance do-

main decomposition methods for structural mechanics Part II: formu-

lation enhancements, multiple right-hand sides and implicit dynamics,

Comput. Methods Appl. Mech. Engrg., Vol. 19, 4611-4662, 2004.

[14] I.T. Carstea, D. Carstea, A.A. Carstea: A domain decomposition ap-

proach for coupled fields in induction heating devices, 6th WSEAS Inter-

national Conference on System Science and Simulation in Engineering,

Venice, Italy, November 21-23 2007.

[15] E.S. Shina, Y.U. Kima, S.N. Jung: Efficient solution algorithms for

thermo-mechanical contact problems based on a domain/boundary de-

composition method, International Journal of Non-Linear Mechanics,

Vol. 44, 610-622, 2009.

[16] C. Farhat and M. Geradin: On the general solution by a direct method

of a large-scale singular system of linear equations: application to the

analysis of floating structures, Int. Jnl. Num. Meth. Engng., Vol. 41,

675696, 1998.

[17] D. Roncoroni: Metodi di decomposizione di dominio applicati a pro-

blemi di meccanica di solidi e strutture, Tesi di Laurea Specialistica,

Politecnico di Milano, 2010.

BIBLIOGRAFIA 211

[18] C.A. Felippa, K.C. Park: Staggered transient analysis procedures for

coupled mechanical systems: formulation, Comput. Methods Appl.

Mech. Engrg, Vol. 24, 61-111, 1980.

[19] C. Farhat, K.C. Park, Y. Dubois-Pelerin: An unconditionally sta-

ble staggered algorithm for transient finite element analysis of coupled

thermo-elastic problems, Comput. Methods Appl. Mech. Engrg, Vol. 85,

365-439, 1991.

[20] A. Felippa, T. Geers: Partitioned analysis for coupled mechanical

systems, Eng. Comput, Vol. 5, giugno 1988.

[21] R. Ardito, L. Baldasarre, A. Corigliano: On the numerical evaluation

of capacitance and electrostatic forces in MEMS, 10th. Int. Conf. on

Thermal, Mechanical and Multiphysics Simulation and Experiments in

Micro-Electronics and Micro-Systems, EuroSimE, 2009.

[22] S.D.A. Hannot: Modeling strategies for electro-mechanical microsystems

with uncertainty quantification, Ph.D. Thesis, 2010.

[23] S.D.A. Hannot, D.J. Rixen: Pull-in curves determined with monolithic

FEM models, Congresso EuroSimE, Delft University of Technology, 26-

30 Aprile 2009.

[24] V. Rochus, D.J. Rixen, J.C. Golinval: Monolithic modelling of electro-

mechanical coupling in micro-structures, International Journal for

Numerical Methods in Engineering, Vol. 65, 461-493, 2009.

[25] V. Rochus: Finite element modelling of strong electro-mechanical

coupling in MEMS, Ph.D. Thesis, 2006.

[26] R. Ardito, C. Comi, A. Corigliano, A. Frangi: Solid damping in micro

electro-mechanical systems, Meccanica, Vol. 43, num. 4, 419-428, 2007.

[27] A. Andreykiv , D.J. Rixen: Numerical modelling of electro-mechanical

coupling using fictitious domain and level set methods, Int. J. Numer.

Meth. Engng, Vol. 80, 478-506, 2009.

[28] D. Rixen, V. Rochus: Extended finite element for electro-mechanical

coupling, 8th. Int. Conf. on Thermal, Mechanical and Multiphysics

Simulation and Experiments in Micro-Electronics and Micro-Systems,

EuroSimE 2007.

BIBLIOGRAFIA 212

[29] , L. Baldasarre: Metodi computazionali per la soluzione di problemi ac-

coppiati in microsistemi elettro-meccanici, Tesi di Laurea Specialistica,

Politecnico di Milano, 2007.

[30] L. Corradi Dell’Acqua: Meccanica delle strutture, McGraw-Hill Vol

I,II,III.

[31] S. Focardi, I. Massa, A. Uguzzoni: Fisica generale. Elettromagnetismo.,

Ed. Ambrosiana, Milano, 2004.

[32] J.G. Michopuolos, C. Farhat, J. Fish: Modeling and simulation of mul-

tiphysics systems, Journal of Computing and Information Science in

Engineering, Vol. 85, 213, 2005.

[33] A. Beccaluva, M.G. De Donno: Micro-accelerometri capacitivi e ri-

sonanti: problemi di progettazione e modellazione numerica, Tesi di

Laurea Specialistica, Politecnico di Milano, 2009.

[34] C.A. Felippa, K.C. Park, C. Farhat: Partioned analysis of coupled

mechanical systems, Comput. Methods Appl. Mech. Engrg, Vol. 190,

3247-3270, 2001.

[35] C.C. Zienkiewicz, R.L. Taylor, J.Z. Zhu: The finite element method: its

basis and fundamentals, Sixth edition.

[36] V. Rochus, D.J. Rixen, J.C. Golinval: Electrostatic coupling of MEMS

structures: transient simulations and dynamic pull-in, Nonlinear

Analysis, Vol. 6, 5-7, 2005.

[37] M.A. Crisfield: Non-linear finite element analysis of solids and struc-

tures. Volume 1: Essentials, John Wiley and Sons, New York,

1991.

[38] A. Corigliano, R. Ardito: Appunti delle lezioni di Meccanica

Computazionale II, Politecnico di Milano, 2009.

[39] C. Farhat, M. Lesoinne, P. Le Tallec: Load and motion transfer algori-

thms for fluid/structure interaction problems with non-matching discre-

te interfaces: momentum and energy conservation, optimal discretiza-

tion and application to aero-elasticity, Computer methods in applied

mechanics and engineering, 157, no. 1-2, 95114, 1998.

BIBLIOGRAFIA 213

[40] X. Jiao and M.T. Heath: Common-refinement-based data transfer bet-

ween non-matching meshes in multiphysics simulations, International

Journal Numerical Method Engng, 61, 24022427, 2004.

[41] A. Andreykiv and D.J. Rixen: Numerical modelling of electro-

mechanical coupling using fictitious domain and level set methods, In-

ternational Journal for Numerical Methods in Engineering, 80, no. 4,

2009.

[42] N.R. Aluru and J. White: An efficient numerical technique for electro-

chemical simulation of complicated microelectromechanical structures,

Sensors Actuators: A. Physical, 58, no. 1, 111, 1997.

[43] S. Mukherjee, Z. Bao, M. Roman, N. Aubry: Non.linear mechanics of

mems plates with a total lagrangian approach, Computers Structures,

83, 758768, 2005.

[44] R. Gilbert, R. Legtenberg, S.D. Senturia: 3D coupled electro-mechanics

for MEMS: applications of cosolve-EM, Micro Electro Mechanical

Systems, Proceedings IEEE, 1995.

[45] E. Brusa, M.Gh. Munteanu: Analisi agli elementi finiti del compor-

tamento dinamico nonlineare di microattuatori elettrostatici: soluzio-

ne del problema accoppiato, Politecnico di Milano, XXXIV convegno

nazionale, 14-17 settembre 2005.

[46] F. Shi, P. Ramesh, S. Mukherjee: Dynamic analysis of micro-electro-

mechanical systems, International Journal for Numerical Methods in

Engineering, Vol. 39, 4119-4139, 1996.

[47] E. Riks: An incremental approach to the solution of snapping and buc-

kling problems, International Journal of Solids and Structures 15 , no. 7,

529551, 1979.

[48] M.A. Crisfield: A fast incremental/iterative solution procedure that

handles snapthrough, Computers and Structures, Vol. 13, no. 1-3, 5562,

1981.

[49] N.R. Aluru, J. White: An efficent numerical technique for

electro-mechanical simulation of complicated microelectro-mechanical

structures, Sensor and Actuators, A 58, 1-11, 1997.

BIBLIOGRAFIA 214

[50] N.R. Aluru, J. White: Algorithms for coupled domain MEMS simu-

lation, Research Laboratory of Electronics, Massachusetts Istitute of

Tecnology, Cambridge.

[51] A. Quarteroni, R. Sacco, F. Saleri: Matematica numerica, 3a Edizione,

Springer Milano, 2008.

[52] GiD. User Manual.

[53] Comsol Multiphysics. Modeling Guide.

METODI DI DECOMPOSIZIONE DEI DOMINI APPLICATI ALLA ... · 3.8 In alto la curva statica di pull-in,...

Documents

Transcript of METODI DI DECOMPOSIZIONE DEI DOMINI APPLICATI ALLA ... · 3.8 In alto la curva statica di pull-in,...