Analisi Isogeometrica. Vantaggi derivanti dall'uso di funzioni di forma a continuità elevata per problemi monodimensionali

Introduzione
Ad oggi, nel panoramea dei metodi di maggior successo nell’ambito della meccanica
computazionale, al classico metodo degli elementi niti (FEM), nell’ultimo decen-
nio, si  e aancata l’Analisi Isogeometrica (IGA).
Il metodo degli elementi niti trae le sue origini negli anni ’50 da alcuni articoli di
Turner [1] e Argyris e Kelsey [2], che raccolsero una serie di articoli sull’argomento
pubblicati tra il 1954 ed il 1956. Esso si fonda sull’idea di studiare un problema
dierenziale partendo dalla sua formulazione debole; il concetto basilare del metodo
 e quello di approssimare la geometria del problema e la sua soluzione con un set
di funzioni note. Questo permise agli ingegneri di studiare strutture sempre pi u
complesse, potendo disporre di uno strumento di calcolo molto pi u potente rispetto
al tradizionale calcolo manuale. Questo ruolo cos  importante per la progettazione
port o il FEM ad un rapido progresso; in particolare, nella seconda met a degli anni
’60 vennero introdotti gli elementi isoparametrici [3, 4] che fecero fare un notevole
passo in avanti al metodo. Il problema era che questi elementi valevano solo per
funzioni di forma di classe C
0
; vennero introdotte nuove tecniche per estendere il
metodo a funzioni di forma con maggiore regolarit a, ma il risultato era troppo poco
eciente. Inoltre, rimaneva il problema della corretta denizione della geometria
mediante le sole funzioni di forma scelte.
A met a degli anni ’70 cominci o a svilupparsi l’idea di poter risolvere entrambi que-
sti problemi scegliendo come funzioni di forma le funzioni (altamente regolari) usate
per la rappresentazione CAD dell’elemento da studiare (di cui si esporr a in dettaglio
nel capitolo 2). Quest’idea non trov o piena realizzazione, poich e tali funzioni erano
molteplici e non tutte adatte alle tecniche di analisi numerica; tuttavia lo svilup-
parsi dell’IGA ha spinto anche gli sviluppatori delle rappresentazioni CAD a venire
incontro alle esigenze degli analisti strutturali cosicch e i due mondi stanno in questi
anni convergendo, secondo i piani iniziali.
In questa tesi considerer o due casi circa un’asta in stato tenso-deformativo monodi-
mensionale con lo scopo di:
 calcolare gli spostamenti assiali sotto un carico distribuito (per esempio sinu-
soidale) ed i relativi errori
11
Analisi Isogeometrica Luca Marioni
 calcolare le autofrequenze proprie dell’asta stessa ed i relativi errori
Una volta risolti con il classico FEM, e introdotta l’IGA, verranno confrontati i ri-
sultati.
Nelle prime due appendici verranno invece accennate alcune tecniche di migliora-
mento dell’IGA che sfruttano l’elevata continuit a delle funzioni di forma. Nell’ultima
appendice verranno presentati i codici MatLab con cui ho svolto le analisi. Questi
sono ovviamente da considerarsi a puro scopo didattico, atti solo ad evidenziare al-
cuni risultati importanti: sono stati quindi esclusi dalla loro stesura problemi quali
l’ecienza dei programmi stessi e l’inserimento dei parametri strutturali, seppure
ingegneristicamente rilevanti.
Con questa tesi si vuole mostrare come l’IGA tragga vantaggi rilevanti dalla re-
golarit a delle sue funzioni di forma.
 E per questo che sono stati scelti solo problemi
monodimensionali, in modo da escludere ogni contributo derivante dalla denizione
geometrica del problema (nonostante proprio questa sarebbe l’elemento basilare e
di maggior ecacia dell’Analisi Isogeometrica).
12
Capitolo 1
Metodo degli elementi niti
Il metodo degli elementi niti  e storicamente il metodo di analisi numerica pi u uti-
lizzato;  e stato oggetto di studi e, ad oggi, si presenta implementato per i pi u svariati
utilizzi, con tecniche varie e altamente strutturate.
In questo capitolo, per mantenere l’argomentazione pi u lineare possibile, senza di-
vagazioni eccessive, verranno trattati solo i princ pi base del metodo in funzione
della loro applicazione nei codici MatLab utilizzati. Inne, come gi a accennato
nell’Introduzione, verranno considerati solo problemi 1D.
1.1 Problema statico
Il problema statico considerato  e quello di un’asta in stato tensodeformativo mono-
dimensionale sollecitata assialmente; esso  e governato dall’equazione di Navier 1D
con due coppie di condizioni al contorno: entrambe sulla funzione soluzione agli
estremi (di Dirichlet) , oppure una sulla funzione soluzione e una sulla sua derivata
(di Neumann).
8
>
<
>
:
 u
00
(x) =p(x)
u(a) =c
u(b) =d
8
>
<
>
:
 u
00
(x) =p(x)
u(a) =c
u
0
(b) =d
(1.1)
Dove u(x) sono gli spostamenti assiali e p(x) =
q(x)
EA
, con q(x) il carico assiale di-
stribuito, E modulo di Young , A area della sezione dell’asta (considerando EA
costante lungo l’asta).
Dal problema nella sua forma forte (dierenziale) si passa quindi alla forma debole
(integrale), moltiplicando poi per una qualsiasi funzione tale da dare signicato al-
l’integrale, per esempio continua.
 Z
b
a
u
00
(x)v(x)dx =
Z
b
a
p(x)v(x)dx 8vregolari
13
Analisi Isogeometrica Luca Marioni
Le funzioni v usate vengono denominate funzioni test. Dal punto di vista sico
possono esser considerati spostamenti virtuali, seguendo l’interpretazione del metodo
come principio di minimo dell’energia. Mantenendo un approccio pi u strettamente
matematico, il passaggio  e giusticato dal fatto che non si considera una v specica,
ma un insieme indenito di funzioni test. Integrando per parti e ipotizzando per ora
una doppia condizione di Dirichlet omogenea, si ottiene
Z
b
a
u
0
(x)v
0
(x) =
Z
b
a
p(x)v(x) 8v2V eu2V
u
dove V
u
=fu2C
0
(a;b); u(0) =c; u(b) =dg e V =fu2C
0
(a;b); u(0)u(b) = 0g.
Il problema  e che la funzione soluzione andrebbe testata con innite funzioni di-
verse, dal momento che dim(V ) =1. Da qui l’idea di Galerkin di approssimare
lo spazio delle funzioni test con un sottospazio V
h
di dimensione nita, denominato
appunto spazio degli elementi niti . Questo sottospazio contiene tutte le funzioni di
una certa struttura (per esempio polinomi di grado 1 a tratti). Le funzioni scelte
come base dello spazioV
h
dovranno dare signicato all’integrale su tutto il dominio,
da cui deriva la non necessariet a di esser di classe C
0
in tutto l’intervallo (a;b), ma
solo di avere propriet a di integrabilit a. In questo modo posso sostituire le funzioni
v (di struttura generica) con un numero nito di funzioni base.
Sar a comodo costruire le funzioni base in modo tale che i coecienti di tale com-
binazione lineare siano anche la soluzione nodale approssimata del problema. Que-
sta richiesta, che verr a meno nel metodo isogeometrico, prende il nome di base
interpolante e s’identica con la struttura denita dall’equazione (1.2)
 i
(x
j
) = ij
(1.2)
dove  ij
 e l’operatore di Kronecker.
Per cui gli spostamenti incogniti saranno deniti come
u
h
(x) =
N
X
j
 j
(x) ^ u
j
Dove N  e il numero di nodi usati per la discretizzazione, u
h
sono gli spostamenti
nodali approssimati e (x) le funzioni scelte come base. Esprimendo le funzioni v(x) e
u(x) come combinazione lineare delle funzioni base, si ottiene quindi la formulazione
debole discreta del problema (1.1)
N
X
j=1
^ u
j
Z
b
a
 0
j
 0
i
dx =
Z
b
a
p(x) i
dx (1.3)
Oppure, in forma matriciale:
K ^ u=F
14
Capitolo 1. Metodo degli elementi niti
dove K  e detta matrice di rigidezza, ^ u  e il vettore degli spostamenti nodali appros-
simati e F  e il vettore proiezione dei carichi sulle funzioni base.
1.1.1 Applicazione delle condizioni al contorno
Nella pratica non risulta conveniente seguire il procedimento rigoroso per imporre
le condizioni al contorno. Il modo pi u comodo per imporre le condizioni di Dirichet
sulla funzione u degli spostamenti nodali  e modicare il sistema risolvente (equazione
(1.1)) in modo tale che la soluzione rispetti le condizioni al contorno. In particolare si
usa sostituire gli elementi della prima e dell’ultima riga della matrice di rigidezza con
zeri tranne, rispettivamente, il primo e l’ultimo elemento della riga (corrispondenti
al primo e ultimo nodo dell’asta); dopodich e si pongono il primo ed ultimo elemento
del vettore F uguali al valore degli spostamenti nodali negli estremi, come si vede
nell’equazione (1.4)
0
B
B
B
B
B
B
@
1 0 0 : :
k
2;1
k
2;2
: : k
2;n
: : : : :
: : : : :
k
n 1;1
k
n 1;2
: : k
n 1;n
: : 0 0 1
1
C
C
C
C
C
C
A
0
B
B
B
B
B
B
@
u
1
:
:
:
:
u
n
1
C
C
C
C
C
C
A
=
0
B
B
B
B
B
B
@
c
F (2)
:
:
F (n 1)
d
1
C
C
C
C
C
C
A
(1.4)
Per quanto riguarda le condizioni di Neumann (ipotizzandola, per semplicit a, sull’e-
stremo b ed aancata da una condizione di Dirichlet omogenea sul primo estremo)
si procede all’integrazione per parti dell’equazione (1.1)
Z
b
a
u
0
(x)v
0
(x)dx [u
0
(x)v(x)]
b
a
=
Z
b
a
p(x)v(x)dx
Z
b
a
u
0
(x)v
0
(x)dx =
Z
b
a
p(x)v(x)dx +dv(b)
N
X
j=1
^ u
j
Z
b
a
( 0
j
 0
i
)dx =
Z
b
a
p(x) i
dx +d
i
(b)
Dal momento che, per la condizione di base interpolante, i
(b)6= 0 solo per un valore
dell’indice i corrispondente all’estremo b (in cui si ha  (b) = 1), la condizione di
Neumann s’imporr a sommando il valore dela derivata di u all’estremo b all’ultimo
elemento del vettore dei carichi, come si vede nell’equazione (1.5)
0
B
B
B
B
@
1 0 0 : :
k
2;1
k
2;2
: : k
2;n
: : : : :
: : : : :
k
n;1
k
n;2
: : k
n;n
1
C
C
C
C
A
0
B
B
B
B
@
u
1
:
:
:
u
n
1
C
C
C
C
A
=
0
B
B
B
B
@
c
:
:
:
d +F (n)
1
C
C
C
C
A
(1.5)
15
Analisi Isogeometrica Luca Marioni
1.1.2 Funzioni di forma lineari
Le funzioni di forma lineari che rispettano la richiesta di base interpolante sono
denite, ragionando sull’elemento generico [ x
i
;xi + 1] di ampiezza h, secondo l’e-
quazione (1.6) e sono descritte visivamente in Fig.1.1. In particolare si noti che
sul singolo elemento vivono solo due funzioni di forma consecutive, come si vede in
gura 1.2.
(
 i
= x x
i+1
h
 i+1
=
x x
i
h
(1.6)
Figura 1.1: Funzioni di forma po-
linomiali di grado 1 su tutto il
dominio
Figura 1.2: Funzioni di forma po-
linomiali di grado 1 sul singolo
elemento
1.1.3 Funzioni di forma quadratiche
Le funzioni di forma polinomiali di grado due richiedono un grado di libert a in pi u
rispetto a quelle lineari, che viene determinato con la creazione di un sotto-nodo
all’interno dell’elemento (in genere al centro). Considerando di costruirle sull’in-
tervallo elementare [-1;1] e presupponendo la richiesta di base interpolante, esse
vengono denite dalle equazioni (1.7) e si rappresentano come in Figura 1.3.
8
>
>
>
<
>
>
>
:
 i
=
x
2
2
 x
2
 i+1
= x
2
+ 1
 i+2
=
x
2
2
+
x
2
(1.7)
Le funzioni di forma cos  costruite vengono poi trasportate sul particolare elemento
dopo l’integrazione moltiplicando per lo Jacobiano della trasformazione ane neces-
saria. Un altro metodo  e quello di eettuare il cambio di coordinate prima dell’inte-
grazione: in questo modo le funzioni di forma vengono rappresentate analiticamente
16
Capitolo 1. Metodo degli elementi niti
Figura 1.3: Funzioni di forma quadratiche
(Equazione (1.8)) da un polinomio di ordine 2.
8
>
>
>
>
>
<
>
>
>
>
>
:
 i
=
2(x x
i+1
)(x x
i+2
)
h
2
 i+1
=
 4(x x
i
)(x x
i+2
)
h
2
 i+2
=
2(x x
i
)(x x
i+1
)
h
2
(1.8)
dove si considera il generico elemento [x
i
;x
i+2
] [x
i
;x
i
+h], x
i+1
 e il nodo interno
centrale.
1.1.4 Studio dell’errore di consistenza del metodo agli ele-
menti niti
Uno schema, per essere convergente, deve essere consistente e stabile (Teorema di
Lax). In questa tesi useremo lo studio della consistenza (espressa analiticamente
in (1.9)). Elemento interessante di quest’analisi sar a la valutazione di p, ordine di
convergenza nei vari problemi trattati.
8k : lim
h!0
jy
k
 y(t
k
)j = 0 ^jy
k
 y(t
k
)j Ch
p
(1.9)
Valutazione dell’ordine di convergenza con funzioni di forma lineari
Nel caso di funzioni di forma lineari, la consistenza si trova come segue.
Si parta dal problema statico gi a analizzato e riassunto in (1.10).
(
 u
00
(x) =f(x) in[0; 1]
u(0) =u(1) = 0
(1.10a)
17
Analisi Isogeometrica Luca Marioni
Z
1
0
u
0
(x)v
0
(x)dx =
Z
1
0
f(x)v(x)dx8vV (1.10b)
Si introduca lo spazio V
h
denito come segue:
V
h
=
 v
h
2C
0
(0; 1) : v
h
I
2P
1
(I)8I; v
h
(0) =v
h
(1) = 0
 !V
h
 V (1.11)
dove I  e il generico intervallo della suddivisione.
Procedendo come visto in (1.1) si arriva alla scrittura delle relazioni (1.12a), solu-
zione debole discreta, e (1.12b), soluzione debole continua.
Z
1
0
u
0
h
(x)v
0
h
(x)dx =
Z
1
0
f(x)v
h
(x)dx 8v
h
2V
h
 V (1.12a)
Z
1
0
u
0
(x)v
0
h
(x)dx =
Z
1
0
f(x)v
h
(x)dx 8v
h
2V
h
 V (1.12b)
Quindi, sottraendo la (1.12b) alla (1.12a) si ottiene la relazione (1.13) che prende il
nome di ortogonalit a di Galerkin .
Z
1
0
(u(x) u
h
(x))
0
v
0
h
(x)dx = 0 8v
h
 V
h
(1.13)
Questa  e vera solo se non si commettono errori nel calcolo degli integrali. Si noti,
inne, che v
0
h
 e, nel caso considerato, una qualunque funzione costante a tratti (in
quanto derivata di funzioni di forma lineari).
Si denisca quindi
jju u
h
jj
2
L
2
:=
Z
1
0
(u(x) u
h
(x))
2
(1.14)
Si voglia quindi valutare l’errore dello schema numerico degli elementi niti, espresso
come segue:
      (u u
h
)
0
      2
L
2
=
Z
1
0
(u u
h
)
0
(u u
h
)
0
=
Z
1
0
(u u
h
)u
0
 Z
1
0
(u u
h
)
0
u
0
h
(1.15)
Essendo
R
1
0
(u u
h
)
0
u
0
h
= 0 per l’ortogonalit a di Galerkin (1.13), si pu o sostituire
u
0
h
con una qualunque funzione lineare, in particolare u
0
I
, funzione di interpolazio-
ne lineare degli spostamenti nodali. Si ottiene quindi la relazione (1.16) e, per la
disuguaglianza Cauchy-Schwartz, si trova la (1.18).
      (u u
h
)
0
      2
L
2
=
Z
1
0
(u u
h
)
0
(u u
I
)
0
(1.16)
18