Sul metodo degli elementi finiti applicato a problemi di elasticità piana
by user
Comments
Transcript
Sul metodo degli elementi finiti applicato a problemi di elasticità piana
MECH.1 - Modellistica e Approssimazione Numerica
di Problemi in Meccanica dei Continui
http://mox.polimi.it
MO X
MECH.1 - Modellistica e Approssimazione Numericadi Problemi in Meccanica dei Continui – p.
Contenuti della lezione
Il problema elastico modello: energia potenziale
Formulazione debole del problema elastico
Esistenza e unicità della soluzione debole
Formulazione agli spostamenti
Modelli bidimensionali per l’elasticità
Discretizzazione della formulazione agli spostamenti
Forma matriciale del problema discretizzato
Validazione numerica
Il caso di materiali incomprimibili
MECH.1 – p.
Riferimenti bibliografici
1. T.J.R. Hughes, The Finite Element Method - Linear static and dynamic
finite element analysis, (Cap. 2 e 4), Prentice-Hall (1987).
2. O.C. Zienkiewicz, The Finite Element Method, Mc-Graw Hill (1977).
3. S.C. Brenner, L.R. Scott, The Mathematical Theory of Finite Element
Methods, (Cap. 9), TAM 15, Springer-Verlag New York (1994).
4. D. Braess, Finite Elements: Theory, Fast Solvers, and Applications in
Solid Mechanics, (Cap. 6), Cambridge University Press, Cambridge
(1997).
5. R. Dautray, J.L. Lions, Mathematical Analysis and Numerical Methods
for Science and Technology, (Cap. 1), Vol. 1, Springer-Verlag Berlin
Heidelberg (1990,2000).
6. PdeTool Manual, The MathWorks Inc.
Inoltre, per ulteriori riferimenti al modello matematico del problema
dell’elasticità lineare nel caso di materiali comprimibili e incomprimibili, si può
fare riferimento alle Lezioni ELL1 e MECH2 di questo Corso.
MECH.1 – p.
Il problema elastico modello: energia potenziale
Consideriamo un corpo elastico lineare isotropo che occupa il dominio
, vincolato sul bordo
(con vincolo omogeneo), sottoposto
sulla parte di bordo
(
,
) e soggetto alla
al carico
forza di volume .
%
(-&
(1)
e
8
7
/
.
/
, con
;
:;* 9
<
: campo di spostamento tale che
7
.
dove
120 543
6 .
03
03 /
,+*
"
#
'(&
)
#$
!
Introduciamo il funzionale energia potenziale totale del sistema costituito dal
corpo e dalle forze e
: campo di deformazione associato a
.
8
tale che
?
30 !
=
>
?
03 !
7
!
e tale che
?
03
7
30 8
=
>
? : campo di sforzo associato a
;
MECH.1 – p.
! F
(3)
"A
)
EDC
!
BA)
A
@
Equazione costitutiva:
, : modulo elastico e modulo di Poisson del materiale;
@
A
dove
(2)
03 !
Equazione cinematica:
' 3 0
)
' 0? 3
Relazioni costitutive
8
tensore identità. Introducendo i parametri di Lamè
00 .
7
G
e
6
40
.
CD
F
,
)
GI
KL J
F
! M
è un tensore simmetrico e definito positivo per
AN
O
<
.
P
(4)
(5)
M
dove
o in forma compatta
!H
! F
)
DEC G
H!
H
del materiale si ha
MECH.1 – p.
Formulazione agli spostamenti
%
Sostituendo le (5) e (4) nella (1), otteniamo la seguente espressione
dell’energia potenziale
&? (-&
#
'(&
)
%
R+*
H! !
)
I G
KL J
Q
"
#
*+
"
#
)
&'
#
! )
GI
KLJ
!H
S
Il campo di spostamento corrispondente alla configurazione di equilibrio
assunta dal corpo elastico sotto l’azione dei carichi e dei vincoli esterni è la
soluzione del seguente problema
X TRUJ
VW
Y
(6)
? ]
su
:R* [ 6
\
7
[6
9
<
Z
nello spazio
MECH.1 – p.
S
La formulazione debole del problema (6) si ottiene annullando la variazione
prima di rispetto a e valutandola in
, e diventa:
^
?7
Z
(-&
(7)
R+*
"
#
#
('&
&'
!
tale che
! M Z
S 7
S
cercare
In Meccanica dei Continui, il problema (6) è noto come Principio di Minimo
dell’Energia Potenziale, mentre il problema (7) è il Principio dei Lavori Virtuali
(cf. lezione ELL1).
MECH.1 – p.
Esistenza e unicità della soluzione del pb. elastico
^
7
Z
)
GI
KLJ
SI
KL J
('&
H!
S ! ('&
! d
c
Z
! M
S è il funzionale lineare e continuo
? (-&
*+
"
#
#$
('&
`
# `
Z
c
d
mentre
(8)
è la forma bilineare
S
_
## _ baZ
dove
` S
_
S
7
Z
Il problema debole (7) è della forma astratta:
tale che
cercare
Z
## _ Per dimostrare l’esistenza e l’unicità della soluzione di (8) è necessario
sia coerciva sullo spazio . Ciò è
verificare che la forma bilineare
possibile grazie ad un’importante proprietà, detta disuguaglianza di Korn.
MECH.1 – p.
La disuguaglianza di Korn
<
f hg e
un dominio aperto con bordo regolare a tratti. Esiste allora una
tale che
^
?7
Z
99
!
('&
i
f 99
(9)
! Y Q
f
7
Sia
costante
Z
## _ è coerciva sullo
La (9) permette di concludere che la forma bilineare
spazio e pertanto, applicando il Lemma di Lax-Milgram, che il problema (8)
ammette soluzione unica.
j 99
mkl 99
kn
p?oo
N
Y
S 99
99
<
gj
S
Da ciò segue anche la stima a priori per la soluzione , ovvero che esiste una
costante
tale che
MECH.1 – p.
Forma forte del problema dell’elasticità
in
su
su
div
S S
# < S
)
<
tur
tsr q
S
S' Il problema forte corrispondente alla formulazione debole (8) è il seguente:
tale che
cercare
(10)
S )
GI
KLJ
SF
!H
S dove
v
|
zyx
w
${x
(cf. lezione ELL1).
MECH.1 – p. 1
Modelli bidimensionali per l’elasticità (I)
Sotto opportune ipotesi, è possibile semplificare le equazioni del modello
elastico 3D ottenendo due differenti modelli 2D, noti come modello plane
stress e modello plane strain.
Plane stress (stato di sforzo piano) - ipotesi:
il corpo ha una dimensione (spessore) trascurabile rispetto alle altre due
(es.: piastra);
il carico è applicato nel piano di simmetria dello spessore e non ha
componenti trasversali rispetto a tale piano.
z
L
y
L>>t
x
t
MECH.1 – p. 1
<
.
.
(
<
~
'}
Assumendo che il piano
coincida con il piano di simmetria dello spessore e
che l’asse sia orientato perpendicolarmente ad esso, si ha che
-
dove
! P
"A
<
<
tur sr q
- ! !
tr r
<
A
<
Q
A
@
tur
"A
sr q
.
tr r
e
.
è in generale non nulla e
? ! )
! "A
"
!
A
!
Osservazione: la componente di deformazione
può essere post-calcolata come
MECH.1 – p. 1
Modelli bidimensionali per l’elasticità (II)
Plane strain (stato di deformazione piana) - ipotesi:
si suppone che il corpo abbia una dimensione molto maggiore rispetto
alle altre due (es.: trave snella, cilindri, alberi,...);
non ci sono spostamenti lungo l’asse della trave, ma solo in direzione
perpendicolare all’asse stesso;
e hanno componente nulla
~
~
tutti i carichi applicati sono indipendenti da
rispetto all’asse .
L>>h
y
h
L
x
z
MECH.1 – p. 1
coincida con l’asse di simmetria si ha che
?<
-
-
!
<
~
Assumendo che l’asse
In questo caso la legge costitutiva (3) rimane quella valida nel caso generale,
previa opportuna eliminazione delle righe e delle colonne le cui componenti
sono nulle.
? è in generale non nulla e può
)
A
Osservazione: la componente di sforzo
essere post-calcolata come
MECH.1 – p. 1
Discretizzazione del problema elastico in 2D
tale che
(11)
^
7
? Z
&-
,+*
"
#
#
('&
('&
! M
S
!
S
7
cercare
Z
L’approssimazione Galerkin-elementi finiti del problema (7) diventa:
0 ' }
,
'} 6
40
0
0 S
0 ' }
6
40
}
S' ,
7
le componenti del campo di spostamento
Su ciascun elemento
possono essere espresse come
0
;
delle rispettive funzioni di base definite su
: restrizioni su
;
5
I
J
T dove
.
0
0 S
: gradi di libertà dell’approssimazione relativi all’elemento
MECH.1 – p. 1
<
b
n
¡ b ¡ b <
tr r
&
? ¢&
¡ b < ¡ b
b <
n
b b
n n
b b
<
utr sr q
!6 QQ ! 66 !
Q
tr r
tur sr q
b < b
Utilizzando la relazione cinematica (2) otteniamo la seguente espressione per
:
il tensore di deformazione
!
S
S Q
trrtrrrrtrr
utrrtrrrrtrr
Q
6
? &
<
<
<
tsrtrrtrrtrrr q
SQ 6 6 S
trtrrtrrtrrr
<
< Q
< 6
S
6
Limitiamo nel seguito la nostra attenzione al caso di elementi finiti .
il vettore degli spostamenti dei nodi
Indicando con
dell’elemento
(numerati secondo la usuale convenzione antioraria), si ha:
&
S6 6
Q S
Q
S
Forma matriciale (I)
MECH.1 – p. 1
Forma matriciale (II)
7
^
&
¢
M
¢
'(&
£
_
S £
¤
£
&
Nell’assemblaggio del sistema lineare associato a (11) è conveniente lavorare
a livello del singolo triangolo di griglia (cf. lezione IMPL2). Si ha
è la matrice di rigidezza dell’elemento
.
£
¤
dove
('&
£
¥
£
Il contributo al termine noto dovuto alla forza di volume è
? (-&
1 1
§£
+R*
¦
£
mentre quello dovuto al carico di superficie è
MECH.1 – p. 1
e
A
@
Consideriamo una piastra quadrata di lato unitario con
sottoposta a uno stato di trazione come indicato in figura.
?<
Patch test: stato di sforzo costante (I)
g N=1
g N=1
?<
.
S
@ } "A
@ "A
'
La soluzione esatta del problema è:
6¨
Il metodo numerico (che utilizza elementi finiti ) deve essere in grado di
risolvere esattamente (a meno della precisione macchina) il caso test in
oggetto.
MECH.1 – p. 1
Patch test: stato di sforzo costante (II)
Impostiamo con Pdetool il caso test proposto. Disegnamo un quadrato di
lato unitario con il lato sinistro disposto sull’asse delle ordinate e selezioniamo
i parametri nel modo seguente:
Application: Structural Mechanics Plane Stress.
.
,
A
@
PDE Coefficients:
?<
Boundary Conditions: sui lati inferiore e sinistro fissiamo
condizioni di Dirichlet imponendo il valore di u e v esatti, mentre sui lati
superiore e destro assegnamo lo sforzo normale, di valore unitario.
Usiamo una griglia molto rada (ad esempio 4 elementi) e risolviamo il
problema (cliccando sul bottone =). Quindi:
esportiamo i parametri della mesh (da Export Mesh);
esportiamo la soluzione (da Export Solution) che è rappresentata
da un vettore di lunghezza 2*nv, dove nv=size(p,2) è il numero dei
vertici della triangolazione;
estraiamo le componenti U=u(1:nv); V=u(nv+1:end) e facciamone
un disegno con pdesurf: il campo di spostamento lineare è
rappresentato in maniera esatta.
MECH.1 – p. 1
Patch test: stato di sforzo costante (III)
Per calcolare il tensore degli sforzi, è necessario:
}
ea
del campo di spostamento
S
'
calcolare le derivate parziali rispetto a
approssimato ;
calcolare le componenti del tensore di deformazione utilizzando la (2);
utilizzando la legge costitutiva dello stato di sforzo piano, calcolare lo
sforzo.
6¨
La funzione pdegrad calcola il gradiente dello spostamento a partire
(tale tensore è
dall’espressione del gradiente delle funzioni di base
costante su ).
La funzione pdesmech (che internamente utilizza pdegrad) calcola le
componenti del tensore di deformazione e di sforzo assegnando in ingresso i
parametri della mesh p,t, il tensore di rigidezza c (ottenuto da PDE, Export
PDE coefficients), le componenti della soluzione [U;V] e il nome della
quantità da calcolare.
sxx=pdesmech(p,t,c,[U;V],’tensor’,’sxx’);
syy=pdesmech(p,t,c,[U;V],’tensor’,’syy’);
txy=pdesmech(p,t,c,[U;V],’tensor’,’sxy’);
Anche il campo di sforzo è rappresentato in questo caso in maniera esatta,
essendo un tensore costante.
MECH.1 – p. 2
Caso test 2 (I)
?<
A
@
Consideriamo come secondo caso test una piastra quadrata di lato unitario
e
sottoposta a uno stato di trazione come indicato in figura
con
e vincolata allo spostamento orizzontale e verticale rispettivamente sui lati
sinistro e inferiore.
g N,x= 2
g N,y= 6x
g N,y= 2
g N,y= 2y
g N,x= 2
@ "A
? '}
@ } )A
Q
"
S
©
La soluzione esatta del problema è
Risolviamo il problema in modo analogo a prima. Poiché il metodo numerico
utilizza elementi finiti lineari in questo caso non ci aspettiamo di rappresentare
in modo esatto la soluzione.
MECH.1 – p. 2
Caso test 2 (II)
9
'
?
"
9
X T¬
9
'
'
"
S
X T¬
9
S' mª
«k
o
T¬
S
"
S
ª
Risolviamo il problema impostato (con la griglia proposta da Matlab) e, dopo
aver definito l’espressione della soluzione esatta e aver esportato la soluzione
del problema da PdeTool, calcoliamo la norma infinito dell’errore, definita
come
®<
L’errore risulta essere dell’ordine di
, e quindi, come previsto,
l’approssimazione con elementi finiti lineari non è in grado di rappresentare
esattamente una funzione quadratica.
Calcoliamo ora il campo di sforzo a partire dal campo di spostamento:
sxx=pdesmech(p,t,c,[uu;vv],’tensor’,’sxx’)
syy=pdesmech(p,t,c,[uu;vv],’tensor’,’syy’)
sxy=pdesmech(p,t,c,[uu;vv],’tensor’,’sxy’)
Chiaramente l’approssimazione ottenuta è costante a tratti, ovvero non vi è
continuità dello sforzo tra un triangolo e l’altro.
MECH.1 – p. 2
Post-calcolo del campo di sforzo
Come possiamo migliorare la rappresentazione del campo di sforzo?
S
il campo di spostamento
che abbiamo calcolato è una funzione
con gradi di libertà nei nodi della griglia;
lineare a tratti e continua su
S
vorremmo, pertanto, ottenere una rappresentazione del campo di sforzo
nel medesimo spazio di approssimazione utilizzato per .
Possiamo applicare la seguente procedura di ricostruzione del gradiente (da
cui poi calcolare le componenti della deformazione e quindi dello sforzo):
£
S¯
S ¯
calcoliamo
sulla triangolazione . Tale tensore risulta essere
costante su ogni elemento
e sarà indicato come
. Possiamo
attribuire tale valore al baricentro dell’elemento ;
di
0°
nel nodo
calcoliamo
X3
²,±
,³£
S¯
S ¯
per ricostruire il valore di
3 9
R²±
9
S ¯
0 °
(12)
.
0°
è l’insieme degli elementi che afferiscono al nodo
0´
dove
La routine pdeprtni effettua esattamente questa procedura. I risultati
prodotti sono quelli disegnati dal comando Plot di PdeTool.
MECH.1 – p. 2
Il caso di materiali quasi-incomprimibili (I)
Alcuni materiali di interesse industriale (ad esempio la gomma) sono
quasi-incomprimibili, ovvero è necessario spendere molta energia per
causarne un cambiamento di volume. Equivalentemente, un materiale quasi
incomprimibile sottoposto a uno stato di sforzo tende a deformarsi senza
modificare apprezzabilmente il proprio volume.
(13)
dove
Z
Z
µ
Z
EDC
I
LKJ
S
<
!
S
Z
µ
Matematicamente, la condizione di incomprimibilità si esprime come
rappresenta la deformazione volumetrica.
? H
)
"A
)A
A @
)A
H
H
)
G
¶
@
¶
È possibile (e interessante) formulare la (13) in funzione dei parametri di Lamè
del materiale introducendo il modulo di taglio e il modulo di elasticità
volumetrico
MECH.1 – p. 2
Il caso di materiali quasi-incomprimibili (II)
" A )A
¶
·
Definendo il rapporto
c ¸
G
c
? B¹)
A¸
c
B¹)
·
(13)
¸
segue che
(14)
»
½ V X URJ
Y½ ¼
ª
S
"
ª
Y
Y
S
"
S
N
ª
ª
º
G
Valori elevati del parametro rendono il problema elastico di ardua risoluzione
numerica. Dal lemma di Ceá si ha infatti che
## _ Hf
,
»
è la costante di coercività di
;
## _ è la costante di continuità della forma bilineare
H
)
º
G
dove
»
ºP
può diventare molto grande nel limite
e dunque segue che il rapporto
incomprimibile, ovvero l’errore commesso nell’approssimazione può diventare
molto grande (l’errore non è dovuto in questo caso all’uso di una griglia troppo
rada!).
MECH.1 – p. 2
Il caso di materiali quasi-incomprimibili (III)
Il problema del trattamento numerico di materiali nel limite incomprimibile si
manifesta attraverso l’insorgere del cosiddetto fenomeno di volume locking.
A
A titolo di esempio, mostriamo in figura il comportamento dell’errore di
discretizzazione (cf. la (14)) nel caso del problema test precedentemente
esaminato, al variare di .
−1
10
ν=0.49999
−2
10
ν=0.3−0.45
p=1
−3
10
−2
10
−1
10
0
10
Dal punto di vista numerico, il locking si traduce nel fatto che la matrice di
rigidezza della formulazione agli spostamenti diventa estremamente mal
condizionata nel limite incomprimibile. Ciò rende indispensabile adottare
formulazioni alternative.
MECH.1 – p. 2
Formulazione a più campi nel caso incomprimibile
KLJ
?S
"
¾
GI
Per ovviare alle difficoltà nel trattamento del caso incomprimibile
precedentemente discusse, è necessario modificare la relazione costitutiva (4)
introducendo la pressione idrostatica:
¾' ¾
div
in
in
su
su
turrr
? strrr
I
S S KLJ
# < S
<
S )
<
q
S
S' In questo modo la formulazione del problema elastico nel caso incomprimibile
diventa: cercare
,
tali che
(15)
"
?¿¾ F
!H S
S
dove
MECH.1 – p. 2
La (15) è nota come formulazione di Herrmann del problema dell’elasticità e
rappresenta un classico esempio di problema a più campi (in questo caso,
spostamento e pressione idrostatica).
È significativo osservare che la medesima formulazione (15) costituisce il
modello matematico del moto viscoso di un fluido viscoso (equazioni di
Stokes, saranno trattate nel Corso di Fluidodinamica Numerica).
La formulazione debole corrispondente alla (15) conduce alla risoluzione di un
problema di punto-sella, a differenza di quanto accade nel caso della (10) a cui
corrisponde un problema di minimo.
Ciò ha rilevanti conseguenze dal punto di vista numerico nella scelta di una
approssimazione stabile e convergente per la coppia spostamento-pressione
(velocità-pressione nel caso del problema di Stokes). Tale scelta è subordinata
alla condizione di compatibilità di Babuska-Brezzi (cf. Corso di Fluidodinamica
Numerica a cui rimandiamo per una discussione degli elementi finiti stabili
econvergenti per l’approssimazione di (15)).
MECH.1 – p. 2