Analisi ad elementi di contorno di lastre di Reissner

Analisi ad elementi di contorno
di lastre di Reissner
A. Leone
M. Aristodemo
Dipartimento di Strutture
Universitµ
a della Calabria
87030 Arcavacata di Rende (Cosenza)
Sommario
Questo report riguarda l'analisi ad elementi di contorno di lastre spesse,
basate sul modello continuo proposto da Reissner.
Il lavoro riporta la descrizione del modello continuo, la formulazione
integrale sul contorno e gli aspetti computazionali che coinvolgono il processo
di discretizzazione. Sono riportati i risultati dell'integrazione analitica, con
le espressioni complete degli integrali per la soluzione sul contorno e nel
dominio.
Le potenzialitµ
a del modello numerico sono state valutate attraverso una
serie di test numerici ed i risultati sono stati confrontati con altre soluzioni
ottenute sia con il metodo degli elementi di contorno che con il metodo
degli elementi ¯niti. Un particolare attenzione µe stata rivolta all'analisi dei
risultati relativi allo strato di bordo.
1
Introduzione
L'analisi di lastre in°esse caratterizzate da deformabilitµa tagliante interessa
numerose applicazioni della meccanica strutturale. L'interesse per questo argomento ha ricevuto un notevole impulso dallo sviluppo dei mezzi di calcolo
e dai problemi di modellazione posti dalle tecnologie piµ
u recenti.
Il modello meccanico piµ
u di®uso per la simulazione di lastre caricate
trasversalmente µ
e quello, sviluppato da Kirchho® [1], che descrive bene il
comportamento delle lastre sottili. La sua popolaritµa deriva dall'essenzialitµa
della formulazione matematica, che consente in casi semplici di trovare
soluzioni in forma analitica.
1
La necessitµ
a di de¯nire modelli che abbiano carattere generale e che permettano di analizzare la varietµa di dati che si presentano nelle applicazioni
tecniche, ha motivato le ricerche della prima metµa del Novecento sulle lastre
dotate di deformabilitµ
a tangenziale. I lavori piµu rilevanti su questo argomento sono quelli di Reissner [2][3] e Mindlin [4].
La strada piµ
u di®usa per la simulazione numerica delle lastre caricate
trasversalmente µe rappresentata dal metodo degli elementi ¯niti. Questo
metodo, basato su formulazioni variazionali e su rappresentazioni approssimate di tipo locale, µ
e stato oggetto di moltissime ricerche. Contributi signi¯cativi in questo ambito sono stati proposti da [5], [6] e [7].
Un modo alternativo per costruire una soluzione numerica del problema
delle lastre in°esse µe rappresentato dal metodo degli elementi al contorno. Il
metodo consente di formulare direttamente un problema con sole variabili di
contorno, usando la conoscenza della soluzione analitica in problemi de¯niti
su domini piµ
u semplici di quello in esame. Infatti, quando esiste la funzione
di Green del problema di®erenziale, il problema integrale sul dominio puµo
essere trasferito sul contorno geometrico della lastra senza introdurre approssimazioni.
Il metodo degli elementi al contorno di lastre in°esse sottili µe stato studiato da Hartmann [8] [9] ed altri autori nei primi anni 70. L'analisi ad
elementi di contorno di lastre spesse µe stata avviata dal 1982 con Van der
Weeen [10] che ha pubblicato un primo lavoro sull'argomento, fornendo
anche la soluzione fondamentale necessaria per la formulazione integrale
sul contorno. Piµ
u di recente sono stati pubblicati nuovi lavori [11], [12],
che riguardano il processo di discretizzazione delle lastre di Reissner con il
metodo degli elementi di contorno.
Il metodo degli elementi di contorno utilizza la forma discretizzata della
equazione integrale sul contorno per de¯nire un sistema algebrico nelle sole
variabili di contorno . L'analisi µe quindi capace di fornire soluzioni accurate
con poche risorse di calcolo ed µe caratterizzata da uno sviluppo progressivo del processo di soluzione che fornisce in una prima fase la soluzione
sul contorno e successivamente la soluzione nei punti di interesse del dominio. Quest'ultimo aspetto risulta estremamente vantaggioso, rispetto al
metodo degli elementi ¯niti, perchµe spesso non si ha interesse alla soluzione
sull'intero dominio ma all'interno di regioni limitate dove si veri¯cano le
condizioni che regolano il dimensionamento dell'elemento strutturale.
La formulazione integrale sul contorno delle lastre di Reissner utilizza
una soluzione fondamentale espressa attraverso le funzioni di Bessel modi¯cate. La presenza di queste funzioni richiede che i nuclei siano rappresentati
2
in forma approssimata. In particolare, la rappresentazione delle funzioni
di Bessel µe stata organizzata con di®erenti approssimazioni per piccoli e
grandi argomenti, usando un criterio automatico di selezione del valore di
transizione tra le due rappresentazioni.
L'integrazione dei coe±cienti del sistema sul contorno e delle quantitµa
necessarie alla valutazione della soluzione all'interno del dominio µe stata
sviluppata attraverso processi essenzialmente analitici. Viene cosi superata la problematica indotta dalla presenza di singolaritµa e quasi singolaritµ
a. La discretizzazione delle variabili di contorno µe e®ettuata attraverso
un'interpolazione quadratica di tipo B-spline che assicura la continuitµa della
funzione e della derivata con un solo parametro per elemento.
2
Il modello di lastra di Reissner
Il modello di lastra spessa di Reissner puµo essere derivato assumendo una
legge di variazione lineare lungo lo spessore della lastra delle componenti di
tensione
12
Mx z
¾x =
h3
12
M z
(1)
¾y =
h3 y
12
¿xy =
M z
h 3 xy
dove le componenti normali ¾x e ¾y quella tangenziale ¿xy sono espresse in
funzione dei momenti °ettenti Mx , M y e di quello torcente Mxy .
Le espressioni delle rimanenti componenti di tensione restano de¯nite
dalle equazioni di equilibrio nel dominio e sulle super¯ci di contorno § h2 .
Ã
µ
¿xz =
3Tx
2h
¿yz =
3Ty
2h
Ã
¾z =
3
q
4
2z 1
¡
h
3
Ã
1¡
1¡
µ
2z
h
2z
h
µ
¶2 !
¶2 !
2z
h
(2)
¶3 !
dove Tx e Ty sono le componenti di taglio.
La precedente descrizione del campo di tensione e le equazioni costitutive
dell'elasticitµ
a lineare implicano che le equazioni di compatibilitµa non possono
3
z, w
O
y, v
q
h/2
h/2
s
n
x, u
Figure 1: Variabili e sistema di riferimento
essere veri¯cate esattamente. Infatti, le componenti di tensione assunte
impongono che le componenti normali di deformazione ² x e ²y , e conseguentemente u, v e ° yz , siano funzioni cubiche di z. Questo µe in contrasto con
l'assunta variazione lineare della componente di tensione tangenziale ¿xy .
Questo errore nel soddisfacimento puntuale dell'equazione di compatibilitµa
puµ
o essere di®uso lungo lo spessore della lastra attraverso una descrizione
bi-dimensionale del campo di spostamento. Questa descrizione µe derivabile
dalla condizione di stazionarietµa dell'energia complementare totale su una
porzione di lastra - delimitata dal contorno ¡.
¦c =
1
2E
Z Z
-
+h
2
¡ h2
(¾x2 + ¾ 2y + ¾ 2z ¡ 2º(¾x ¾y + ¾y ¾ z + ¾x ¾ z ) +
2
2
2
2(1 + º)(¿xy
+ ¿xz
+ ¿yz
))dzd-
¡
Z Z
+h
2
¡ ¡ h2
(¾n u + ¿nt v + ¿nz w)dzd¡
dove n e t denotano il versore normale e tangenziale al contorno.
Ora, introducendo le assunte distribuzione delle tensioni, si trova
¦c
=
1
2E
Z
-
(
12
2
(Mx2 + M y2 ¡ 2ºM x My + 2(1 + º)Mxy
)+
h3
4
(3)
(4)
12(1 + º) 2
12º
(Tx + Ty2) ¡
q(Mx + M y ) +
5h
5h
¡
Z
¡
Z
+ h2
¡h
2
¾ 2z dz)d-
(M n Án + M nt Át + Tn w)
~ d¡
(5)
dove Án , Át , w
~ sono le variabili di spostamento bi{dimensionali sul contorno.
Per ricavare le espressioni Án , Át e w~ si considera l'integrale che esprime
il lavoro delle tensioni sul contorno:
Z Z
¡
+ h2
¡ 2h
(¾n u + ¿nt v + ¿nz w) dzd¡
(6)
Sostituendo le espressioni delle componenti di tensione, si ottiene
12
h3
Z Z
¡
+ h2
¡h
2
(M n u + v + Tn w) zdzd¡
(7)
e confrontando le espressioni dei due integrali sul contorno si ricavano le
espressioni delle variabili bidimensionali di spostamento sul contorno
Án =
12
h3
Át
12
h3
=
w
~ =
Z
¡ h2
Z
3
2h3
+ h2
+ h2
¡ h2
Z
uzdz
(8)
vzdz
(9)
+ h2
¡ h2
w(h 2 ¡ 4z 2 )dz
(10)
La condizione di stazionarietµa dell'energia complementare totale forniscono le equazioni di compatibilitµa da cui si ottengono le relazioni costitutive
Tx
=
Ty
=
M xy
=
Mx
=
My
=
5
Gh(w~;x + Áx )
6
5
Gh(w~;y + Áy )
6
1¡º
D(Áx;y + Áy;x )
2
¶
µ
3º
q
D Áx;x + ºÁy;y +
5Gh
µ
¶
3º
D Áy;y + ºÁx;x +
q
5Gh
5
(11)
(12)
(13)
(14)
(15)
dove
D=
Eh3
12(1 ¡ º 2 )
(16)
Le equazioni di equilibrio in termini di spostamento assumono la forma
0
D(1 ¡ º) B
B
@
2
r ¡ ¸2 +
³
1+º
1¡º
´
simm:
@2
@x2
r¡
³
´
1+º
@2
1¡º ³ @xy ´
@2
¸2 + 1+º
1¡º @y 2
@
¸2 @x
@
¸2 @y
¡¸2 r
1
0
1 0 1
0
C Áx
C @ Áy A = @ 0 A
A
w
q
dove
¸2 =
10
h2
(18)
rappresenta il parametro caratteristico della lastra.
La soluzione del sistema di®erenziale de¯nisce il campo di spostamenti
della lastra quando sono assegnate le condizioni al contorno che possono
riguardare le tre componenti cinematiche sul contorno, w,
~ Án e Át , e le tre
componenti statiche, Tn , Mn e M nt , con la possibilitµa di distinguere le due
situazioni
² vincolo soft: M nt = 0;
² vincolo hard: Át = 0.
Si noti che quando lo spessore della lastra tende a zero, la soluzione a
distanza su±ciente dal bordo tende a quella fornita dal modello senza scorrimenti di Kirchho®, mentre sul contorno rimane distinta per e®etto della
maggiore articolazione delle condizioni al contorno. Dovendo raccordare il
comportamento interno con le condizioni al contorno, la soluzione presenta
gradienti tanto piµ
u elevati quanto piµu piccolo µe il rapporto tra spessore e
lunghezza. Questo e®etto, che riguarda lo strato di bordo, µe stato oggetto di
studio da Arnold e Falk [13] e, con un approccio piµ
u ingegneristico, da Haggblad e Bathe [5]. In questi lavori si de¯nisce in forma asintotica la soluzione
vicino al bordo attraverso sviluppi in serie evidenzandone la dipendenza
dallo spessore.
6
(17)
Nel seguito si riportano sinteticamente i risultati conseguiti, per diverse
condizioni di vincolo, esplicitando la dipendenza delle caratteristiche di sollecitazione e degli spostamenti dallo spessore. Si osservi che il gradiente
della soluzione nello strato di bordo µe tanto piµ
u pronunciato quanto piµ
l'esponente di h µ
e algebricamente piccolo:
² appoggio soft:
T¿ (h ¡1 ); Tº (h 0 ); Mºº (h1 ); Mº¿ (h0 ); M ¿ ¿ (h 1 ); Á¿ (h 1 ); Áº (h2 )
² appoggio hard:
T¿ (h 0 ); Tº (h 1 ); M ºº (h 2); Mº¿ (h1 ); M ¿ ¿ (h 2 ); Á¿ (h2 ); Áº (h3 )
² incastro soft:
T¿ (h 1 ); Tº (h 2 ); M ºº (h 3); Mº¿ (h2 ); M ¿ ¿ (h 3 ); Á¿ (h3 ); Áº (h4 )
² incastro hard:
T¿ (h 0 ); Tº (h 1 ); M ºº (h 2); Mº¿ (h1 ); M ¿ ¿ (h 2 ); Á¿ (h2 ); Áº (h3 )
² lato libero:
T¿ (h ¡1 ); Tº (h 0 ); Mºº (h1 ); Mº¿ (h0 ); M ¿ ¿ (h 1 ); Á¿ (h 1 ); Áº (h2 )
Questi e®etti, che risultano molto pronunciati nella condizione di bordo
libero e appoggio soft, saranno analizzati nei test numerici.
3
Soluzione fondamentale
La formulazione integrale sul contorno richiede la conoscenza di una soluzione
analitica del problema di®erenziale che possa essere utilizzata, partendo da
una forma integrale sul dominio, per trasferire gli integrali sul contorno. In
particolare, si puµ
o utilizzare la soluzione relativa ad una lastra illimitata
soggetta ad un carico puntiforme. La forma di Navier del problema della
lastra di Reissner che de¯nisce la soluzione fondamentale, u ¤, in questo caso
diventa
¢ij u¤kj = ±(r)±ik
(19)
dove r µ
e la distanza tra il generico punto di campo e il punto sorgente
in cui sono applicate le due componenti del carico concentrato unitario e le
componenti dell'operatore di®erenziale ¢ij sono de¯nite dalle relazioni
7
Ã
2
2
¢®¯
= Dc1 (r ¡ ¸ )±®¯
¢®3
= ¢3® = ¡Dc1¸2
¢33
= Dc1 ¸2r 2
c2 @ 2
+
c1 @x ®x ¯
!
(20)
@
@x®
(21)
(22)
dove sono state introdotte le seguenti costanti
µ
¶
1¡º
2
µ
¶
1+º
=
2
c1
=
c2
(23)
(24)
Nelle precedenti espressioni gli indici latini variano da 1 a 3 e quelli greci
si riferiscono alle componenti bidimensionali e variano da 1 a 2.
La soluzione del sistema puµo essere espressa nella forma [14]
u¤j i = ¢co
ij e[r]
(25)
dove ¢co
ij sono i cofattori di posizione (i; j) dell'operatore ¢:
¢co
®¯
2
2
= ¸ D c1
Ã
4
r ±®¯ + (¡1)
(®+¯)
³
2 2 2
2
= ¢co
3® = ¡D ¸ c1 r ¡ ¸
¢co
33
= D2 c1 r2 ¡ c1 ¸2 (r 2 ¡ ¸2 )
´
@2
r ¡ c1(r ¡ ¸ )
@x® x ¯
2
2
´ @
2
¢co
®3
³
Ã
@x ®
2
!!
(26)
(27)
(28)
e la funzione e[r] deve soddisfare l'equazione di®erenziale
det(¢ ij )e([r] = ±[r]
(29)
In particolare, riformulando il problema in coordinate polari e mettendo
in conto la condizione di simmetria radiale, la funzione e(r) deve soddisfare
la relazione
±(r)
D3 c21 ¸2 r4 (r 2 ¡ ¸2 )e(r) = ¡
(30)
2¼r
8
la cui soluzione risulta
1
e[r] =
2¼D 3 ¸6
µ
2
1¡º
¶2 "
K0 [z] +
Ã
!
µ
z2
r
+ 1 ln
4
r0
¶
+1
#
(31)
dove K0 (z) µe la funzione di Bessel modi¯cata di ordine 0 e la costante r 0
rappresenta una distanza arbitraria. E' stata inoltre introdotta la distanza
adimensionale rappresentata dalla variabile
z = ¸r
(32)
Dalla soluzione e(r) si possono derivare le componenti u¤ij della soluzione
fondamentale
u ¤®¯
u¤®3
u¤33
= ¢co
®¯ e[r]
½
1
2
=
(A0 [z]±®¯ ¡ A1[z]r;® r;¯ )
2¼D 1 ¡ º
·
µ ¶
¾
¸
±
r
r r
¡ ®¯ 2 ln
+ 1 ¡ ;® ;¯
4
r0
2
¤
co
= ¡u3® = ¢®3e[r]
·
µ ¶
¸
rr;®
r
2 ln
=
+1
8¼D
r0
= ¢co
e[r]
33
1
= ¡
2¼D¸2
(
·
2
ln
1¡º
µ
r
r0
¶
¸
(33)
(34)
µ
z2
r
+1 ¡
ln
4
r0
¶)
(35)
dove A0 [z] e A1 [z] sono combinazioni delle funzioni di Bessel modi¯cate K 0
e K1
K1 [z]
1
¡ 2
z
z
K1 [z]
2
A1 [z] = K0 [z] + 2
¡ 2
z
z
A0 [z] = K0 [z] +
9
4
Forma integrale sul contorno
La formulazione integrale sul contorno µe dedotta attraverso il teorema di
reciprocitµ
a di Betti che stabilisce l'uguaglianza del lavoro tra due sistemi
equilibrati di forze. Considerando insieme al sistema di forze di volume qj
e di super¯cie pj , che producono gli spostamenti uj , il sistema costituito da
una azioni concentrate agenti nel punto sorgente », espressi dalle funzioni
delta di Dirac ±i[»; x], e dagli spostamenti prodotti nel punto di campo
x, indicati con u¤ij [»; x], il teorema di Betti fornisce l'equazione integrale
espressa in funzione delle sole variabili di contorno
u i[»] +
Z
¡
p¤ij [»; x]uj [x]d¡ =
Z
¡
u¤ij [»; x]pj [x]d¡ +
Z
-
u¤ij [»; x]q j [x]d-
(36)
dove si µe ottenuto conto della relazione
Z
-
±i[»; x]ui [x]d- = ui [»]
Nel caso della lastra di Reissner la forma integrale sul contorno risulta
ui [»] +
Z
Z¡
-
p¤ij [»; x]u j [x]d¡
q[x]
µ
u¤i3 [»; x]
=
¡
Z
¡
u ¤ij [»; x]pj [x]d¡ +
¶
º
u ¤ [»; x] d(1 ¡ º)¸2 ik;k
(37)
dove u ¤ij raccoglie le due componenti di rotazione (j = 1; 2) e quella di
traslazione (j = 3) indotti dalle due componenti di coppia (i = 1; 2) e
di forza (i = 3). I termini p¤ij rappresentano le componenti di momento
(j = 1; 2) e di taglio (j = 3) sul contorno.
Le componenti della soluzione fondamentale u ¤ij sono rappresentati dalle
relazioni
u ¤®¯
u¤®3
u¤33
½
1
2
=
(A0 [z]±®¯ ¡ A1[z]r;® r;¯ )
2¼D 1 ¡ º
·
µ ¶
¸
¾
±®¯
r
r;® r;¯
¡
2 ln
+1 ¡
4
r
2
·
µ 0¶
¸
rr;®
r
=
2 ln
+1
8¼D
r
( 0
µ ¶)
· µ ¶
¸
r
1
2
r
z2
= ¡
ln
+1 ¡
ln
2¼D¸2 1 ¡ º
r0
4
r0
10
(38)
(39)
(40)
dove µ
e stata indicata con r la distanza tra il punto sorgente ed il punto di
campo. Le componenti p¤ij , derivate da queste attraverso le equazioni di
compatibilitµ
a e di legame, assumono la forma
µ
¶
2º ¤
1¡º
u¤i®;¯ + u ¤i¯;® +
(41)
u ±®¯ º¯
2
1 ¡ º i°;°
´
1¡º 2³ ¤
(42)
p¤i3 = D
¸ ui® + u¤i3;® º®
2
dove le componenti º® rappresentano le componenti della normale nel punto
di campo.
Le espressioni esplicite dei termini p ¤ij sono
p¤i®
= D
p¤®¯
p¤®3
p¤3®
p¤33
1
[(4A2 ¡ 4A1 + c0 )(±®¯ r;º + r;® º¯ )
4¼r
+(4A1 + c1)r;® º¯ ¡ 2(2A2 + c0 )r;® r;¯ r;º ]
¸
=
[A0 º® ¡ A1 rº ]
2¼ ·µ
µ ¶¶
¸
1
r
= ¡
c2 + 2c1 ln
º® + 2c0 r® rº
8¼
r0
r
= ¡ ;º
2¼r
= ¡
(43)
(44)
(45)
(46)
avendo posto
c0 = 1 ¡ º
c1 = 1 + º
c2 = 1 + 3º
ed essendo stata introdotta la combinazione delle funzioni di Bessel
A2 (z) = K0 (z) +
µ
8
1
K 1(z) ¡
z
z
¶
+ zK1 (z)
L'integrale sul dominio del termine di carico puµo essere ricondotto, utilizzando il teorema della divergenza, ad un integrale sul contorno per varie
leggi di carico. Facendo riferimento a un carico uniforme e posto u¤i3 = vi;®®
(® = 1; 2) si ha
Z µ
¶
º
u ¤ [»; x] d(1 ¡ º )¸2 i®;®
¶
Z µ
º
¤
u
[»;
x]
º® d¡
= q
v¤i;®[»; x] ¡
(1 ¡ º )¸2 i®
¡
qi = q
u ¤i3[»; x] ¡
11
(47)
I termini della funzione ausiliaria vi necessaria a trasformare l'integrale
di dominio del carico in integrale di contorno sono
5
v®¤
=
v¤3
=
µ
µ
¶
¶
r
r 3 r®
4 ln
¡1
128¼D
r0
µ µ µ ¶
¶
µ ¶¶
r2
r
1
32
r
2
ln
z
¡
¡
ln
2
128¼D¸
r0
2
1¡º
r0
(48)
(49)
Rappresentazione integrale delle sollecitazioni
Partendo dalla forma integrale che esprime il campo di spostamenti
u i[»] =
Z
¤
u ij [»; x] pj [x]d¡ ¡
Z¡ µ
¡
¤
vi;®
[»; x]
Z
¡
pij [»; x]¤ u j [x]d¡ +
¶
º
u¤ [»; x] º® qd¡
¡
(1 ¡ º)¸2 i®
(50)
le componenti di sollecitazione sono derivate attraverso le relazioni puntuali
di legame
M ®¯
T®
µ
2º
1¡º
u ®;¯ + u ¯;® +
u°;° ±®¯
2
1¡º
1¡º 2
= D
¸ (u® + u3;® )
2
= D
¶
+
ºq
±
(1 ¡ º )¸2 ®¯
dopo aver introdotto la forma integrale delle derivate del campo di spostamento
M ®¯ [»] =
Z
+ q
T® [»] =
Z
+ q
¡
Z
u ¤®¯i [»; x]pi [x]ds ¡
¡
¡
Z
¤
w®¯
[»; x]d¡ +
¤
w3®
[»; x]d¡
¡
p¤®¯i [»; x]u i[x]d¡
º
q±
(1 ¡ º) ®¯
u ¤3i® [»; x]pi [x]ds ¡
¡
Z
Z
¡
(51)
p¤3i® [»; x]ui [x]d¡
(52)
I termini u¤ij h che intervengono nellle precedenti relazioni sono derivati
dalla soluzione fondamentale attraverso le relazioni
12
u¤3i®
½
2º
1¡º
u¤®i;¯ + u¤¯i;® +
±®¯ u ¤´i;´
2
1¡º
o
1¡º 2n ¤
= D
¸ u®i + u ¤3i;®
2
u ¤®¯i = D
¾
(53)
(54)
ed assumono la forma esplicita
1
[(4A2 ¡ 4A1 + c0 )(±®¯ r;® + ±®° r;¯ )
4¼r
+(4A1 + c1 )±®¯ r;° ¡ 2(2A2 + c0)r;® r;¯ r ;° ]
¸
=
[A0 ±®¯ ¡ A1 r® r¯ ]
2¼ ·µ
¸
µ ¶¶
r
1
c2 + 2c1 ln
±®¯ + 2c0 r® r¯
= ¡
8¼
r0
r;¯
=
2¼r
u¤®¯° = ¡
u¤3®¯
u¤®¯3
u¤3¯3
(55)
(56)
(57)
(58)
Analogamente, i termini p¤ijh sono de¯niti dalle relazioni
p¤®¯i
p¤3i®
½
1¡º
2º
p¤®i;¯ + p¤¯i;® +
± p¤
2
1 ¡ º ®¯ ´i;´
o
1¡º 2n ¤
= D
¸ p®i + p¤3i;®
2
= D
¾
(59)
(60)
ed assumono la forma:
p ¤®¯°
=
Dc0
f(2A2 ¡ 4A1 + c0 )(±°® º¯ + ±°¯ º® ) + (4A1 + c2 )±®¯ º°
4¼r2
¡(A3 ¡ 2A2 + 2c0 )((º® r;¯ + º¯ r ;®)r ;° + (±°® r;¯ + ±°¯ r;® )r;º )
¡2(2A2 + c1 )(±®¯ r;° r;º + º° r;® r;¯ ) + 4(A3 + 2c0 )r® r¯ r ;° r ;º g
p ¤®¯3
p ¤3®¯
p¤3®3
2
(61)
Dc0 ¸
f(A2 ¡ 2A1 )(r;¯ º® + r;® º¯ ) ¡ 2A2 r;® r;¯ r ;º + 2A1 ±®¯ r;º g
4¼r
(62)
2
Dc ¸
= ¡ 0 f(A2 ¡ 2A1 )(±®¯ r;º + r;¯ º® ) ¡ 2A2 r;® r;¯ r;º + 2A1 º¯ r;® g
4¼r
(63)
2
Dc0 ¸
=
f(z 2 A0 + 1)º® ¡ (z 2 A1 + 2)r;® r;º g
(64)
4¼r 2
=
13
dove µe stata introdotta la funzione A3 (z), combinazione delle funzioni di
Bessel
48
A3 (z) = 24K 0(z) +
z
µ
¶
1
K1 (z) ¡
+ 8zK 1 (z) + z 2 K0 (z)
z
(65)
I termini w¤ij necessari per la de¯nizione dell'integrale del carico, sono
de¯niti dalle relazioni
¤
w®¯
¤
w3®
½
1¡º
º
¤
v¤®;°¯ + v¯;°®
(u¤
+ u¤¯°;® )
¡
2
(1 ¡ º)¸2 ®°;¯
¾
·
¸
º
2º
¤
¤
(66)
v
¡
u
±®¯ º°
+
1 ¡ º ´;°´ (1 ¡ º)¸2 ´°;´
½
¾
1¡º 2 ¤
º
= D
¸ v®;° + v3;°® ¡
(u ¤ + u¤3°;® ) º° (67)
2
(1 ¡ º)¸2 ®°
= D
che coinvolgono le derivate della soluzione fondamentale u ¤i j e quelle delle
funzioni ausiliarie vi necessarie per il trasferimento dell'integrale di dominio
sul contorno. Le componenti w¤i j risultano
¤
w®¯
w¤3®
r
f(4 ln z ¡ 3)(c0 (r ;¯ º® + r;® º¯ ) + c2 ±®¯ r;º )
64¼
º ¤
+4(c0 r;® r;¯ + º±®¯ )r;º g ¡
u º
c0¸2 ®¯° °
1
º ¤
=
f(2 ln z ¡ 1)º¯ + 2r;¯ r;º g ¡
u º
8¼
c0 ¸2 3®¯ ¯
= ¡
14
(68)
(69)
6
Discretizzazione dell'equazione integrale
L'equazione integrale
u i[»] +
Z
¡
pij [»; x]¤ uj [x]d¡ =
Z
¡
q
uij [»; x]¤p j [x]d¡ +
Z µ
¡
v¤i;® [»; x] ¡
¶
º
u¤ [»; x] º®d¡
(1 ¡ º)¸2 i®
de¯nisce il campo di spostamento nel generico punto [»] del dominio in funzione degli spostamenti ui [x] e delle sollecitazioni pi[x] sul contorno ¡. Al
¯ne di utilizzare l'equazione integrale sul contorno per costruire un modello
discreto della lastra di Reissner µe necessario assumere una forma approssimata dei campi incogniti sul contorno.
Figure 2: Rappresentazione lineare a tratti del contorno geometrico.
La descrizione geometrica del contorno viene rappresentata attraverso
una linearizzazione a tratti che permette di valutare per via analitica i co15
e±cienti del sistema che forniscono la soluzione sul bordo e quelli necessari
al calcolo della soluzione all'interno del dominio.
I tratti di contorno senza singolaritµa, meccaniche o geometriche, individuano dei macroelementi, discretizzati in elementi all'interno dei quali si
adottano speci¯che interpolazioni delle variabili meccaniche.
La forma discreta delle equazioni integrali viene rappresentata attraverso sommatorie estese al numero di elementi in cui µe stato discretizzato il
contorno.
u i[»] +
n µZ
X
¡h
h=1
µ
Z
n
X
q
h=1
p¤ij [»; x]uj [x]d¡h
¡h
µ
¤
vi;®
[»; x]
¶
=
n µZ
X
¡
h=1
u¤ij [»; x]pj [x]d¡h
¶
º
¡
u ¤ [»; x] º® d¡h
(1 ¡ º )¸2 i®
¶
¶
+
(70)
Per la discretizzazione delle variabili meccaniche sul contorno sono state
utilizzate funzioni costanti a tratti e quadratiche. L'interpolazione costante
comporta notevoli sempli¯cazioni nella simulazione, sia per quanto riguarda
la preparazione del codice e sia nella sua esecuzione. La forma discreta
dell'equazione integrale in questo caso assume la forma
(»)
ui
+
n µZ
X
¡h
h=1
µZ
n
X
q
h=1
¶
p¤ij [»; x]d¡h u hj =
¡h
µ
v¤i;® [»; x] ¡
n µZ
X
h=1
¡
¶
u¤ij [»; x]d¡h phj +
¶
º
u¤ [»; x] º® d¡h
(1 ¡ º)¸2 i®
¶
(71)
L'interpolazione quadratica µe stata costruita attraverso una funzione
spline che riesce ad assicurare la continuitµa della funzione e della sua derivata
nei punti di interconnessione senza far intervenire i valori delle derivata come
paramentri dell'interpolazione [15]. I parametri necessari per ogni macroelemento risultano solo due in piµu di una interpolazione costante a tratti.
La rappresentazione approssimata della generica funzione f (´) µe del tipo
f[´] =
3
X
Ãi [´] fi
(72)
i=1
dove le funzioni di interpolazione, dedotte imponendo la continuitµa nei punti
di interconnessione, assumono le seguenti espressioni
16
f1
f3
f2
x
-a
a
Figure 3: Parametri di interpolazione
1 1
1
¡ ´ + ´2
8 4
8
3 1 2
¡ ´
4 4
1
1 1
+ ´ + ´2
8 4
8
Ã1 [´] =
Ã2 [´] =
Ã3 [´] =
dove ´ µe una ascissa adimensionale de¯nita nel dominio normalizzato [¡1; 1].
Questo tipo di interpolazione sarµa indicata nel seguito come H C (High
Continuity)[15].
Introducendo la rappresentazione approssimata dell'interpolazione delle
variabili meccaniche e del contorno, l'equazione integrale assume la forma
discreta
3 ³
X
Ãk [»] uki
k=1
+
n Z
X
p¤ij [»; x]Ãk [x]d¡ u hk
j
h=1 ¡h
3 X
n Z
X
k=1 h=1 ¡h
n Z
X
q º®
!
=
u¤ij [»; x]Ãk [x]d¡ phk
j +
h=1 ¡h
µ
¶
º
u ¤ [»; x] d¡ (73)
vi;® [»; x] ¡
(1 ¡ º )¸2 i®
Al ¯ne di costruire una relazione che contenga solo variabili di contorno
la precedente equazione integrale puµo essere specializzata adottando una
sorgente collocata sul contorno. Valutando la forma limite per la distanza
r, tra sorgente e punto di campo, che tende a zero, si ottiene
17
3 µ
X
1
k=1
2
Ãk [»] u ki
+
n Z
X
p¤ij [»; x]Ãk [x]d¡ uhk
j
h=1 ¡h
3 X
n Z
X
k=1 h=1 ¡h
n Z
X
q º®
!
=
u¤ij [»; x]Ãk [x]d¡ phk
j +
h=1 ¡h
µ
vi;® [»; x] ¡
¶
º
u¤ [»; x] d¡(74)
(1 ¡ º)¸2 i®
In forma compatta la precedente equazione si puµo scrivere
c[»]u + H [»]¹
u = G[»]¹
p
(75)
Questa forma dell'equazione dipende solo da parametri di spostamento
u
¹ e di sollecitazione p¹ sul contorno e puµo essere utilizzata per generare un
sistema di equazioni per collocazione. La strategia di localizzazione dei
punti di collocazione µ
e basata su una distribuzione uniforme dei punti sul
contorno. In ¯gura 4 µ
e indicata la distribuzione delle sorgenti. Si osservi
che negli elementi terminali del macroelemento sono collocate due sorgenti
nelle posizioni individuate dai parametri ® e ¯.
x
a
b
-a
a
b
a
Figure 4: Localizzazione delle sorgenti sul macroelemento.
Per ogni punto sorgente sono generate tre equazioni, due associate alle
sorgenti coppie ed una corrispondente alla sorgente forza. Si ottiene cosµi il
sistema di equazioni sul contorno
Hu
¹ = Gp¹
(76)
che, sulla base delle condizioni al contorno, viene riorganizzato nella forma
Ax = b
(77)
dove la matrice dei coe±cienti A si presenta tutta piena e non simmetrica
ed il vettore dei parametri incogniti x µe costituito da spostamenti o sforzi a
seconda del tipo di condizioni al contorno.
18
7
Rappresentazione delle funzioni di Bessel
La soluzione fondamentale della lastra di Reissner µe de¯nita attraverso le
funzioni di Bessel modi¯cate di ordine zero ed uno. La valutazione degli
integrali che contengono questi nuclei richiede quindi una rappresentazione
esplicita delle funzioni di Bessel che renda possibile l'integrazione in forma
analitica o numerica.
Le funzioni di Bessel modi¯cate Kn (z) sono de¯nite come soluzioni delle
equazioni di®erenziali
d2 y
dy
+z
¡ (z 2 + n2 )y = 0
(78)
dz 2
dz
Esse possono essere espresse in funzione di quelle di ordine piµu basso
attraverso la relazione
z2
2n
Kn (z) + K n¡1 (z)
(79)
z
mentre le derivate possono essere ottenute attraverso le seguenti regole
K n+1 (z) =
d
K 0 (z) = ¡K1 (z)
dz
d
1
K 1 (z) = ¡K0 (z) ¡ K1 (z)
dz
z
(80)
(81)
L'andamento delle funzioni di Bessel di ordine zero, K0 [z], e di ordine
uno, K 1[z], indicati nei gra¯ci seguenti (Fig. 4 e 5), mostrano un comportamento asintotico delle funzioni per z = 0 mentre per valori grandi
dell'argomento le funzioni tendono a zero. La loro rappresentazione in forma
esplicita µ
e basata usualmente su sviluppi in serie [16]. In particolare, quando
l'argomento z µe su±cientemente piccolo si utilizza lo sviluppo in serie di
potenze
1
K 0p [z] = (¡° ¡ ln( ) ¡ ln(z)) +
2
=
Np ³
X
´
ln( 12 ) + ln(z)
2 ¡ 2°
¡
8
4
C1i + C2i log(z) z 2i
i=0
K 1p [z] =
Ã
1
+
z
Ã
ln( 21 ) + ln(z)
¡ (1 ¡ 2 °)
+
4
2
19
!
z 2 + ::::
(82)
!
z + ::::
(83)
4
3
2
1
1
2
3
4
Figure 5: Gra¯co della funzione K0 (z)
10
8
6
4
2
1
2
3
Figure 6: Gra¯co della funzione K1 (z)
20
4
´
1 X³ i
+
C3 + C4i log(z) z 2i+1
z i=0
Np
=
(84)
dove ° µ
e la costante di Eulero (° = 0:577215::), np il numero di termini dello
(i)
sviluppo in serie ed i coe±cienti Cj sono de¯niti dalle relazioni
C1i
1
= ¡
(i!)2 22i
C2i = ¡
C3i
Ã
° ¡ log(2) ¡
1
(i!)2 22i
1
= ¡
(i + 1)!i!22i+1
C4i =
i
X
1
k=1
k
Ã
!
i
X
1
1
° ¡ log(2) ¡
¡
k
2(i + 1)
k=1
1
(i + 1)!i!22i+1
!
Per valori grandi di z la rappresentazione piµu conveniente µe lo sviluppo
in serie esponenziale che richiede un numero di termini tanto piµ
u limitato
quanto piµ
u sono grandi i valori di z
¡z
K0g [z] = e
=
µ
¼
1
9
1¡
+
+ ::::
2z
8 z 128 z 2
1
r
Ng Ã i !
X
¼ @
C5 A ¡z
1+
e
2z
¡z
K1g [z] = e
=
r
0
r
µ
Ã !1
r
Ng
X
¼ @
C6i A ¡z
e
1+
2z
zi
i=1
(85)
(86)
zi
i=1
3
¼
15
1+
¡
+ ::::
2z
8 z 128 z 2
0
¶
¶
(87)
(88)
dove ng µ
e il numero di termini dello sviluppo in serie esponenziale ed i
coe±cienti C5i e C6i hanno la forma
C5i
C6i
= ¡
= ¡
Qi
j=1 (2j ¡
8i i!
Qi
1)2
¡ 1)2 ¡ 4)
8ii!
j=1 ((2j
21
np nng
3
5
7
9
11
13
15
3
1.69
2.50
3.38
4.29
5.23
6.20
7.18
5
1.77
2.50
3.30
4.14
5.02
5.93
6.86
7
1.91
2.57
3.31
4.10
4.93
5.79
6.67
9
2.08
2.68
3.38
4.12
4.90
5.72
6.57
11
2.27
2.84
3.48
4.18
4.93
5.71
6.52
13
2.47
3.00
3.61
4.28
4.99
5.74
6.52
15
2.69
3.18
3.76
4.40
5.08
5.80
6.55
Table 1: Valori di zcr in funzione di NP e ng
Considerando la necessitµa di utilizzare entrambi gli sviluppi, µe necessario de¯nire il valore di transizione, zc , tra i due tipi di rappresentazione e
scegliere il numero di termini per entrambi gli sviluppi in serie.
Il valore z c da usare per un assegnato numero di termini degli sviluppi
puµ
o essere selezionato tenendo conto dell'errore che viene commesso nella
rappresentazione di entrambi gli sviluppi, individuandolo attraverso la condizione di stazionarietµ
a della quantitµa
e[z] =
°
° °
°
° K 0P + K0G ° ° K1P + K 1G °
°
°+ °
° = min
°
° °
°
K
K
0
(89)
1
dove K 0 e K1 indicano le valutazioni delle funzioni di Bessel di un codice di
manipolazione simbolica (ad esempio Mathematica) che le rappresenta con
un numero elevato di termini. I valori di zc ottenuti da questa condizione
sono forniti nella tabella 1 al variare del numero di termini inclusi nelle
rappresentazioni per piccoli argomenti np e grandi argomenti ng .
Per quanto riguarda la scelta del numero di termini da includere negli
sviluppi in serie si puµ
o considerare l'errore commesso nella valutazione dell'integrale della rappresentazione delle funzioni di Bessel. Nella rappresentazione per piccoli argomenti l'errore µe de¯nito come
e1 =
°
°R Z
R ZC R
°
° CR K (z)dz
K1P (z)dz
0P
° 0
0
° RZ
+ RZ
¡ 2°
°
CR
CR
°
°
K0 (z)dz
K 1 (z)dz
0
(90)
0
mentre nella rappresentazione per grandi argomenti assume la forma
e2 =
°R
°
R inf
° inf
°
° Z CR K 0G (z)dz
K
(z)dz
°
1G
ZC R
° R
+ R inf
¡ 2°
° inf
°
° ZC R K0 (z)dz
°
Z CR K 1(z)dz
22
(91)
np nng
3
5
7
9
11
13
15
3
1.1E-4
6.2E-6
5.8E-7
7.0E-8
1.0E-8
1.6E-9
2.9E-10
5
1.6E-4
6.2E-6
3.9E-7
3.4E-8
3.7E-9
4.6E-10
6.2E-11
7
2.9E-4
8.8E-6
4.2E-7
2.7E-8
2.3E-9
2.2E-10
2.2E-11
9
5.9E-4
1.5E-5
5.8E-7
3.1E-8
2.0E-9
1.6E-10
1.1E-11
11
1.2E-3
2.9E-5
9.7E-7
4.2E-8
2.3E-9
1.5E-10
7.6E-12
13
2.5E-3
6.0E-5
1.8E-6
6.7E-8
3.2E-9
1.8E-10
7.6E-11
15
5.0E-3
1.2E-4
3.4E-6
1.2E-7
4.9E-9
2.4E-11
9.8E-12
Table 2: Valori di e1 in funzione di nP e ng
np nng
3
5
7
9
11
13
15
3
2.9E-4
1.2E-4
5.3E-5
2.6E-5
1.4E-5
8.4E-6
5.2E-6
5
2.2E-5
1.9E-5
7.3E-6
2.8E-6
1.2E-6
5.5E-7
2.7E-7
7
9.0E-5
2.7E-7
1.6E-6
6.3E-7
2.3E-7
8.6E-8
3.5E-8
9
2.4E-4
1.2E-5
2.1E-7
1.5E-7
6.4E-8
2.2E-8
7.8E-9
11
5.4E-4
3.2E-5
1.8E-6
6.2E-8
1.5E-8
7.2E-9
2.5E-9
13
1.2E-3
7.1E-5
4.5E-6
2.9E-7
1.4E-8
1.3E-9
8.5E-10
15
2.4E-3
1.6E-5
9.9E-6
6.7E-7
4.7E-8
2.8E-9
8.6E-11
Table 3: Valori di e2 in funzione di NP e ng
Le tabelle precedenti evidenziano la convenienza di includere nelle rappresentazioni per piccoli e grandi argomenti un numero confrontabile di termini. Infatti si riscontra che, sia per i piccoli che per i grandi argomenti,
l'aggiunta di un termine nello sviluppo produce una riduzione di un ordine
di grandezza nell'errore. Sulla base di questi risultati, per la determinazione
della soluzione sul contorno sono stati adottati i valori np = 5, ng = 5,
zc = 2:50, mentre per l'analisi di dettaglio nei punti interni sono stati usati
i valori np = 7, ng = 7, zc = 3:31. L'inclusione di un numero piµu elevato di
termini nello sviluppo per piccoli argomenti puµo essere utile nelle situazioni
caratterizzate da gradienti elevati nello strato di bordo.
8
Valutazione degli integrali
Nel processo di discretizzazione ad elementi al contorno, la valutazione degli
integrali ha, insieme all'interpolazione utilizzata nella discretizzazione, un
ruolo importante sull'accuratezza dei risultati. In particolare, per e®etto
delle situazioni di singolaritµa e di quasi-singolaritµa non risultano e±caci gli
schemi di integrazione numerica standard. La valutazione degli integrali
23
viene preferibilmente condotta in modo analitico. L'integrazione analitica µe
agevolata dalla approssimazione linearizzata a tratti del contorno. Nel caso
delle lastre di Reissner i soli termini che richiedono l'integrazione numerica sono quelli che contengono le funzioni di Bessel nel dominio dei grandi
argomenti.
La rappresentazione approssimata delle variabili meccaniche comporta
che il termine tipico dell'equazione integrale si presenti nella forma
Z
¡h
u¤ t d¡ =
3 Z
X
a
u ¤Ãi [»] d» ti
¡a
i=1
(92)
dove a µ
e la semilunghezza dell'elemento rettilineo su cui si sta integrando e
Ãi [»] sono funzioni polinomiali in ». Rappresentando la funzione approssimante nella forma
Ãi [»] =
n
X
cij » j
(93)
j =0
dove n µ
e il grado del polinomio interpolante, µe possibile esprimere il termine
integrale nella forma
Z
¡h
¤
u t d¡ =
3 X
n
X
i=1 j =0
cij
µZ
a
¡a
¤ j
u » d»
¶
ti
(94)
E' da notare che il calcolo degli integrali deve esser diversi¯cato a seconda
che si debba utilizzare la rappresentazione delle funzioni di Bessel per piccoli
valori di z, dove si usa l'integrazione analitica, o per grandi valori di z, dove
si usa l'integrazione numerica. La selezione del tipo di rappresentazione µe
realizzata secondo lo schema indicato nella ¯gura 7.
Qualora entrambe le rappresentazioni devono essere usate nell'ambito
dello stesso elemento di integrazione, il generico termine viene calcolato come
somma degli integrali de¯niti sui sottodomini dell'elemento corrispondenti,
sulla base del valore z c giµ
a de¯nito, alle rappresentazioni per piccoli e grandi
argomenti. Utilizzando i simboli indicati in ¯gura 7, il generico integrale
viene rappresentato nella forma
j
U =
Z
a
¡a
¤ j
u » d» =
Z
a1
a0
¤ j
u » d» +
24
Z
a2
a1
¤ j
u » d» +
Z
a3
a2
u ¤» j d»
(95)
n
zc
a
x
t
analitica
numerica
a[0]
-a
numerica
a[2]
a[1]
a[3]
Figure 7: De¯nizione estremi di integrazione
Si noti che la suddivisione dell'integrale comporta un lieve aggravio computazionale, che in°uenza perµo poco i tempi di calcolo poichµe si veri¯ca
solo in un numero limitato di situazioni rispetto a quelle ricorrenti in cui
l'elemento µ
e interessato da un unico tipo di rappresentazione.
L'integrazione analitica µe stata organizzata attraverso l'uso di un sistema
di riferimento locale che consente di esprimere l'integrale in forma compatta.
La procedura di assemblaggio µe quindi preceduta da un'operazione di rotazione delle componenti di momento dal riferimento locale a quello globale.
Facendo riferimento alla ¯g.8, nel sistema di riferimento locale x; y, con
l'asse x coincidente con la direzione dell'elemento, le derivate della distanza
r tra il punto sorgente Q e il punto di campo P si sempli¯cano nelle
y¹
r
x¡x
¹
r
r ;y = ¡
r;x =
A±nchµ
e il generico termine possa essere valutato in una forma computazionalmente conveniente µe necessario introdurre le forme basi che ricorrono negli integrali da trattare analiticamente.
(k)
Ej
(k)
Fj
=
=
Z
x2
x1
x2
Z
x1
xk
dx
(z 2 )j
(z 2 )j xk dx
25
x2
t
n
x
Q
r
y
Gi
x1
n
P
t
W
Figure 8: Sistema di riferimento locale.
Gj(k) =
Z
x2
x1
log(z)(z 2 )j xk dx
dove z = ¸ r.
I risultati dell'integrazione possono essere organizzati in forma ricorrente
in modo da rappresentare un generico termine dell'interpolazione polinomiale e poter, quindi, includere un generico numero di termini della serie che
approssima le funzioni di Bessel.
(k)
La soluzione analitica degli integrali del tipo Ej puµo quindi essere
espressa nella forma
Ej(k)
=
( "
xk¡1
¡
(r 2 )j¡1
2¹
x(j ¡
#x 2
x1
(k¡1)
k)E j
+ (k ¡ 1)(¹
x 2 + y¹2 )E j(k¡2) +
o
1
2j ¡ k ¡ 1
In particolare, quando k = 2j ¡ 1, µe necessario utilizzare la seguente
rappresentazione
(k)
Ej
(k¡2)
= E j¡1
(k¡2)
¡ (x¹2 + y¹2)Ej
26
(k¡1)
+ 2¹
xE j
(k)
La rappresentazione ricorrente degli integrali del tipo Fj
espressa come
(k+1)
Fj
h
= ( x k (r 2 )j+1
ix2
x1
(k¡1)
¡ k(¹
x2 + y¹2 )Fj
(k)
2¹
x (j + k + 1)Fj )
(k)
Gj
=
( ·
µ
+
1
k + 2j + 2
(k)
mentre gli integrali appartenenti al tipo Gj
puµo essere
hanno la forma
¶¸
1 x2 k + 1 (k)
1
(r2 )j xk+1 2 log(r) ¡
Fj ¡
+
2
j x1
2
o 1
(k+2)
(k+1)
¹ Gj ¡1 )
¸2 j (Gj¡1 ¡ x
k+1
La valutazione degli integrali che contengono le funzioni di Bessel nella
rappresentazione per grandi argomenti µe stata sviluppata numericamente.
Questi integrali sono costituti da combinazioni della forma base
~ (k) =
H
j
Z
x2
x1
xj
p dx
ez z i z
I contributi non singolari sono stati valutati utilizzando la tecnica standard di Gauss. Nelle situazioni di quasi-singolaritµa l'integrazione di Gauss
fornisce risultati poco accurati e, talvolta, ina±dabili. In questi casi si ricorre a strategie che, attraverso opportune trasformazioni di variabile, rendono la funzione integranda piµu regolare, migliorando conseguentemente
l'accuratezza dell'integrale. L'e®etto di queste trasformazioni µe sostanzialmente quello di riposizionare i punti di integrazione de¯niti nello schema di
Gauss. A questa categoria appartengono i metodi di Telles [17] e quello di
Huang e Cruse [18]. Mentre il metodo di Telles si rivolge all'integrazione
di termini singolari e quasi singolari, la trasformazione di Huang e Cruse
µe adatta per trattare i contributi quasi singolari. Nel modello numerico
sviluppato µe stato usato questo secondo schema per trattare i contributi
quasi-singolari, mentre i termini singolari, che coinvolgono le funzioni di
Bessel per piccoli argomenti, sono stati trattati analiticamente.
27
Il metodo di Huang e Cruse µe basato sulla trasformazione
x=¡
2 ²(a + ²)
¡²
a+2²+ay
che collega l'ascissa ¯sica x con la posizione trasformata y, dove ² µe la distanza mimima della sorgente dall'elemento e a la semi-lunghezza del dominio
di integrazione.
La trasformazione µ
e valida per n > 1. Quando n = 1 la trasformazione
diventa
x=¡
1
¡²
c1 + c2 y
con
p
p
a+ ²+ m ²
p
p
2 ma+² m²
p
p
m
a+ ²¡ m ²
p
p
2 ma+² m²
m
c1
=
c2
=
dove il coe±ciente m µe de¯nito dalla condizione
·
ln ²
m = max 1; ¡
ln 2d
¸
in cui d µe una costante che deve soddisfare la disuguaglianza
1
p <d
2 ²
m
28
9
Organizzazione dell'algoritmo
La formulazione e la discretizzazione descritte nei precedenti capitoli sono
state utilizzate per sviluppare il programma di calcolo automatico TPLATE
(Thick / Thin PLATE) per l'analisi di lastre poligonali basate sul modello
di Reissner. Il programma µe divisibile nelle seguenti parti
² Input;
² PreProcessing;
² Analisi;
² PostProcessing;
Nella fase di Input vengono letti, attraverso un ¯le di testo, i dati relativi alla struttura da analizzare. Le informazioni necessarie per modellare
la struttura riguardano la sua geometria (posizioni dei vertici dei singoli
macroelementi e spessore) i dati per la discretizzazione (numero di elementi
per ogni macroelemento), le condizioni al contorno, i dati sul materiale
(E e º) ed i dati sul carico trasversale. In questa fase si stabiliscono alcune modalitµ
a dell'analisi, il tipo di interpolazione da usare e le opzioni
sull'output.
Nella fase di Preprocessing vengono svolte le operazioni di pre-elaborazione
dei dati di input. Viene, ad esempio, de¯nita la collocazione dei punti
sorgente in relazione al tipo di interpolazione scelta. Per l'interpolazione
costante a tratti sono collocate le sorgenti al centro di ogni elemento. Per
interpolazione di tipo quadratico sono collocati n + 2 sorgenti sugli n elementi di ogni macroelemento, seguendo il criterio di spaziatura descritto nel
paragrafo precedente.
Nella fase di Analisi viene costruito il sistema di equazioni che de¯nisce
la soluzione sul contorno. Gli aspetti piµu rilevanti di questa fase sono schematizzati nella ¯gura successiva.
29
Collocazione delle
sorgenti sul contorno
Contributo dell'integrale
del carico sul termine noto
Valutazione dei termini
della equazione integrale
i=1..Nl il=1..Nel
Valutazione dei nuclei nel
Gestione della rappresenta-
sistema di riferimento locale
zione delle funzioni di Bessel
Calcolo degli integrali base:
Calcolo degli integrali
E, F, G, K
base per z<zc
Valutazione degli integrali
U* e T*
Calcolo degli integrali
base per z>zc
Assemblaggio dei coefficienti nel sistema A x = b
j=1..Nl jl=1..Nel
Figure 9: Costruzione del sistema risolvente
La soluzione sul contorno µe fornita dal sistema
Ax = b
(96)
dove nella matrice A sono raccolti i coe±cienti corrispondenti ai parametri
incogniti raccolti nel vettore x e b µe il vettore dei termini noti. Alla ¯ne della
fase di Analisi µ
e nota completamente, in termini di parametri dell'interpolazione
scelta, la soluzione sul contorno nelle sue componenti di sforzo e di spostamento.
30
In¯ne, nella fase di PostProcessing viene valutata la soluzione nei punti
interni del dominio in funzione della soluzione sul bordo e del tipo di interpolazione adottata. In particolare quest'ultima fase µe articolata nel modo
indicato nella ¯gura seguente
Soluzione nei punti
interni al dominio
Contributo dell'integrale
del carico
Valutazione dei termini
dell'equazione integrale
Valutazione dei nuclei nel
sistema di riferimento locale
Gestione della rappresentazione delle funzioni di Bessel
Calcolo degli integrali base:
E, F, G, K
Calcolo degli integrali
Valutazione degli integrali
U* e T*
Calcolo degli integrali
base per z<zc
base per z>zc
j=1..Nl jl=1..Nel
Figure 10: Soluzione nei punti interni al dominio
Come si puµ
o osservare, la valutazione degli spostamenti ripercorre le
modalitµ
a seguite nella costruzione del sistema di equazioni. La de¯nizione
delle sollecitazioni interne risulta piµ
u impegnativa in quanto coinvolge le
derivate delle rotazioni e degli spostamenti. Risulta in questo caso ancora piµ
u
determinante la precisione con cui vengono calcolati gli integrali. Anche per
essi si utilizza la procedura generale di gestione basata sulla combinazione
degli integrali base.
31
10
Risultati numerici
I risultati numerici riportati nel seguito permettono di valutare le potenzialitµ
a del modello di discretizzazione sviluppato. I problemi test sono stati
scelti in modo da poter veri¯care le prestazioni del modello a confronto con
soluzioni disponibili e studiare le doti di accuratezza in situazioni piµu impegnative, quale il fenomeno del boundary layer, caratterizzato da elevati
gradienti nella distribuzione delle tensioni.
Nel calcolo degli integrali numerici, sono stati utilizzati in genere 8 punti
di Gauss, mentre per il calcolo della soluzione nelle zone di quasi singolaritµa
ed in presenza di rapporti tra spessore e lunghezze molto ridotti sono stati
usati 16 punti.
10.1
Lastra quadrata appoggiata
La lastra quadrata appoggiata sul contorno riportata in ¯gura 11 µe stata
analizzata sia come lastra sottile (tabella 4) e sia come lastra spessa (tabella
5), rappresentando entrambe le condizioni di appoggio con rotazione tangenziale libera (appoggio sof t) e nulla (appoggio hard). L'analisi, condotta per
diversi tipi di discretizzazione, mostra una buona convergenza alla soluzione
esatta e un'accuratezza soddisfacente anche impiegando un numero di variabili estremamente ridotto.
L
y
C
A
L
x
Figure 11: Lastra quadrata appoggiata.
32
Nel
2
4
8
16
32
[19]
103 W (C) D=(q L4)
hard
soft
0.41002 0.42514
0.40704 0.41113
0.40634 0.40776
0.40625 0.40690
0.40624 0.40667
0.406
My(C) 10=(q L2)
hard
soft
0.48211 0.49462
0.47955 0.48298
0.47895 0.48015
0.47887 0.47944
0.47887 0.47924
0.479
TA/(q L)
hard
soft
0.37484 0.49919
0.32112 0.39578
0.33339 0.41343
0.33730 0.42039
0.33762 0.42055
0.420
Table 4: Lastra sottile h=L = 1=1000
Nel
2
4
8
16
32
103 W (C) D=(q L4 )
hard
soft
0.42645 0.47733
0.42502 0.46388
0.42467 0.45947
0.42454 0.45833
0.42450 0.45817
(C)
My 10=(q L2 )
hard
soft
0.48239 0.52727
0.48113 0.51571
0.48083 0.51189
0.48072 0.51090
0.48068 0.51076
TA /(q L)
hard
soft
0.34453 0.46938
0.33907 0.42000
0.33860 0.42208
0.33851 0.42130
0.33848 0.42120
Table 5: Lastra spessa h=L = 1=10
10.2
Lastra di Razzaque
Il problema presentato in ¯gura 12 µe la lastra sghemba analizzata da Razzaque [20], libera su due lati e vincolata sugli altri con un appoggio di tipo
soft. I risultati nelle tabelle e nei gra¯ci che seguono, indicano la possibilitµ
a di raggiungere soluzioni accurate anche con pochi elementi. I gra¯ci
delle ¯gure successive, ottenuti impiegando 16 elementi per lato, mostrano
l'andamento del momento torcente e del taglio lungo la sezione a indicata in
¯gura 12. I gra¯ci evidenziano l'accuratezza del modello nel cogliere gli elevati gradienti del campo di tensione che si producono vicino al bordo quando
lo spessore della lastra diventa molto piccolo rispetto alla sua larghezza. In
questa situazione si attivano gradienti elevanti in corrispondenza dello strato
di bordo, descritti in dettaglio in [5], che rappresentano un test molto severo
per i modelli numerici.
I risultati delle tabelle 7 e 8 mettono a confronto le soluzioni ottenute
assumendo interpolazioni di tipo diverso per le variabili meccaniche. In
particolare, si puµ
o osservare la maggiore e±cienza dell'interpolazione HC
rispetto all'approssimazione costante a tratti.
33
L
y
C
a
a = 60°
x
3L/4
Figure 12: Lastra di Razzaque.
N el
2
4
8
16
32
[20]
103 w(c) D=(q L4 )
hard
sof t
0.80416 0.82858
0.79219 0.80100
0.79198 0.79487
0.79248 0.79345
0.79270 0.79303
0.7945
0.7945
(c)
my 103 =(q L2 )
hard
sof t
96.72379 98.25051
95.66029 96.05971
95.64432 95.78212
95.68113 95.72701
95.69849 95.71404
95.89
95.89
Table 6: Lastra di Razzaque (h = 0:1): spostamento e momento °ettente
nel punto centrale.
nel
2
4
8
16
[20]
errore
103
Interpolazione Costante
(c)
w (c) D=(q L 4 ) my 103=(q L 2 )
0.81032
90.0667
0.80634
96.7528
0.80733
97.1591
0.80390
96.9631
0.7945
95.89
1.18%
1.12%
103
Interpolazione HC
(c)
L4 ) my 103 =(q L2 )
0.82858
98.25051
0.80100
96.05971
0.79487
95.78212
0.79345
95.72701
0.7945
95.89
0.13%
0.17%
w(c) D=(q
Table 7: Lastra di Razzaque (h = 0:1): spostamento e momento °ettente
nel punto centrale.
34
0.005
0.000
-0.005
h/L:
0.100
-0.010
0.010
0.005
-0.015
0.001
-0.020
Mxy
-0.025
-0.030
-0.035
-0.040
-0.045
-0.050
0.00
0.01
0.02
0.03
0.04
0.05
x/L
Figure 13: Lastra di Razzaque: momento Mxy lungo la linea x = 3L=4.
0.040
h/L:
0.039
0.100
0.010
0.038
0.005
0.001
0.037
Ty
0.036
0.035
0.034
0.033
0.00
0.01
0.02
0.03
0.04
0.05
x/L
35
Figure 14: Lastra di Razzaque: taglio ty lungo la linea x = 3L=4.
h
0.1
1.0
10.0
Interpolazione Costante
(c)
10 w(c) D=(q L4 ) my 103 =(q L 2)
0.80390
96.9631
0.80621
96.1312
0.84577
97.5788
3
Interpolazione HC
(c)
10 w D=(q L 4 ) my 103 =(q L2 )
0.79345
95.72701
0.79935
96.18851
0.84702
98.10232
3
(c)
Table 8: Lastra di Razzaque (h = 0:1): spostamento e momento °ettente
centrali al variare dello spessore (Nel=16).
10.3
Mensola sghemba
Il problema della mensola obliqua di ¯gura 15 µe un test impegnativo per la
presenza di elevati gradienti del campo di tensione vicino agli angoli.
I risultati di alcuni valori dello spostamento trasversale, in corrispondenza degli spigoli liberi della lastra ed al variare dell'angolo di obliquitµa,
sono confrontati in tabella 9 con i valori forniti da modelli ad elementi ¯niti
ra±nati, adatti all'analisi di lastre sia sottili che spesse.
2
1
y
a
x
L
Figure 15: Mensola sghemba (E = 100, º = 0:3, h = 4, q = 1, L = 100).
36
N el
2
4
8
16
32
[21]
Q4BL(32x32)[6]
T3BL(64x64)[7]
® = 700
w1
w2
1.4168 1.0474
1.4313 1.0492
1.4338 1.0467
1.4329 1.0460
1.4320 1.0457
1.4297 1.0442
1.4309 1.0448
1.4309 1.0446
® = 500
w1
w2
1.2045 0.5673
1.2006 0.5551
1.1994 0.5518
1.1966 0.5509
1.1943 0.5505
1.1818 0.5468
1.1869 0.5481
1.1891 0.5485
® = 300
w1
w2
0.9001 0.1751
0.8898 0.1656
0.8801 0.1625
0.8743 0.1616
0.8706 0.1611
0.8502 0.1562
0.8518 0.1572
0.8590 0.1589
Table 9: Mensola sghemba: spostamenti agli spigoli del lato libero.
10.4
Lastra di Hartman
Il problema della lastra rettangolare parzialmente appoggiata lungo il contorno, mostrato in ¯gura 16 µe un caso interessante di condizioni al contorno
discontinue che provocano comportamenti singolari nella soluzione con forti
concentrazioni delle tensioni nell'intorno delle discontinuitµa.
y
4
x
4m
4m
Figure 16: lastra di Hartman con º = 0:3 D = 2700 KNm carico uniforme
p = 100 kN=m2 S = 0:004 m
Nelle tabelle 10 - 13 i valori dello spostamento trasversale e del momento °ettente lungo due sezioni signi¯cative della lastra sono confrontati
con le soluzioni, ottenute per lastre sottili in [8] e [22] con modelli ad ele37
menti di contorno ed in [23] con un modello ad elementi ¯niti ri¯nito, usato
nell'ambito di una strategia di soluzione di tipo multigrid con elevato numero
di variabili.
La soluzione, anche con una discretizzazione piuttosto rada, tende bene
ai valori di riferimento. Il vincolo di appoggio si µe modellato attraverso una
condizione di tipo hard.
Nel
x
-10/3
-8/3
-2.0
-4/3
-2/3
0.0
2/3
4/3
2.0
8/3
10/3
4.0
2
0.02860
0.05746
0.09375
0.15471
0.27129
0.47866
0.79469
1.20992
1.69953
2.23775
2.80384
3.38299
8
0.02712
0.05308
0.08228
0.12623
0.20610
0.35178
0.58837
0.91820
1.32280
1.77762
2.26152
2.75959
16
0.02701
0.05264
0.08089
0.12258
0.19764
0.33482
0.55948
0.87561
1.26602
1.70662
2.17636
2.66034
32
0.02698
0.05247
0.08029
0.12093
0.19375
0.32693
0.54593
0.85546
1.23897
1.67262
2.13543
2.61248
Riferimenti
[23]
[8]
0.027 0.03
0.054
0.082
0.120
0.191 0.20
0.323
0.528
0.832
1.210 1.27
1.619
2.076
2.549 2.65
Table 10: Problema di Hartman: spostamento trasversale w lungo y = 0
38
x
-10/3
-8/3
-2.0
-4/3
-2/3
0.0
2/3
4/3
2.0
8/3
10/3
4.0
2
0.02018
0.04160
0.06735
0.10811
0.19780
0.39135
0.71663
1.14789
1.65117
2.19956
2.77285
3.35829
N el
8
16
0.02015 0.01998
0.03900 0.03856
0.05967 0.05856
0.09031 0.08762
0.14855 0.14204
0.27236 0.25718
0.50732 0.47879
0.84988 0.80734
1.26734 1.21102
1.73015 1.66021
2.21710 2.13363
2.71534 2.61840
32
0.01993
0.03839
0.05809
0.08640
0.13908
0.25022
0.46550
0.78726
1.18422
1.62673
2.09350
2.57163
Riferimenti
[22]
[23]
0.020 0.020
0.038 0.038
0.058 0.057
0.085 0.084
0.136 0.134
0.243 0.243
0.451 0.444
0.766 0.763
1.155 1.158
1.591 1.579
2.050 2.044
2.521 2.521
Table 11: Problema di Hartman: spostamento trasversale w lungo y = 1
x
-10/3
-8/3
-2.0
-4/3
-2/3
0.0
2/3
4/3
2.0
8/3
10/3
2
14.18
-13.23
-96.27
-255.20
-440.93
-534.76
-480.71
-347.54
-211.31
-102.18
-26.66
N el
8
16
22.86
24.12
9.36
12.66
-43.85
-36.72
-148.51 -134.95
-302.65 -282.50
-435.65 -416.67
-445.24 -435.39
-345.63 -342.56
-217.55 -216.96
-108.52 -108.57
-32.38
-32.53
32
24.78
14.27
-33.36
-128.65
-272.99
-407.40
-430.30
-340.81
-216.50
-108.48
-32.51
Riferimenti
[23]
[8]
26.7
18.1
-25.2
-120.0 -133
-256.7
-388.8
-418.5 -433
-334.3
-214.5 -212
-104.3
-35.6
Table 12: Problema di Hartman: momento °ettente mx lungo y = 0
39
La distribuzione di alcune componenti signi¯cative della soluzione µe rappresentata nei gra¯ci a curve di livello delle ¯gure 17 e 18. Si noti la presenza
di forti concentrazione di tensioni nelle zone prossime al punto di cambiamento delle condizioni al contorno.
x
-10/3
-8/3
-2.0
-4/3
-2/3
0.0
2/3
4/3
2.0
8/3
10/3
2
4.989
-2.702
-48.788
-227.569
-532.012
-694.284
-534.208
-350.919
-205.479
-96.211
-25.885
N el
8
16
19.185
21.282
10.139
13.470
-27.960
-21.777
-114.952 -101.482
-309.034 -280.274
-583.489 -550.658
-545.936 -537.398
-366.718 -360.806
-212.293 -205.419
-100.250 -93.699
-27.875
-26.079
32
21.610
14.492
-19.485
-96.380
-268.177
-536.586
-535.384
-360.941
-205.536
-93.763
-26.136
Riferimenti
[22]
[23]
22.2
22.1
16.9
16.5
-16.7
-13.8
-90.6
-88.7
-256.3 -248.3
-523.4 -516.9
-523.5 -533.7
-360.3 -361.6
-205.1 -207.8
-94.4
-91.0
-30.1
-25.3
Table 13: Problema di Hartman: momento °ettente mx lungo y = 1
40
W:
1
0
-1
-3
-2
-1
0
1
2
3
-3
-2
-1
0
1
2
3
-3
-2
-1
0
1
2
3
Mxx:
1
0
-1
Myy:
1
0
-1
Figure 17: Problema di Hartman: spostamento trasversale e momenti M xx,
Myy .
41
Mxy
1
0
-1
-3
-2
-1
0
1
2
3
-3
-2
-1
0
1
2
3
-3
-2
-1
0
1
2
3
Tx
1
0
-1
Ty
1
0
-1
Figure 18: Problema di Hartman: momento M xy e sollecitazione di taglio
Tx e Ty .
42
10.5
Un solaio a piastra
Il problema presentato in ¯gura 10.5, che schematizza un solaio vincolato
su pilastri e soggetto a carico uniforme, mostra l'e±cacia del modello nella
simulazione, con un numero relativamente basso di variabili ( 900 variabili
utilizzando una discretizzatione del contorno esterno in 16 elementi e 4 per
ciascun lato dei pilastri ), di un problema strutturale articolato. Lo spessore
della lastra µ
e pari a 0:20 m ed il coe±ciente di Poisson µe stato assunto pari
a 0:3.
Figure 19: Solaio su pilastri
43
W
12
10
8
6
4
2
2
4
6
8
10
12
2
4
6
8
10
12
Mxx
12
10
8
6
4
2
Figure 20: Solaio su pilastri: spostamento e momento M xx
44
Tx
12
10
8
6
4
2
2
4
6
8
10
12
Mxy
12
10
8
6
4
2
2
4
6
8
10
12
Figure 21: Solaio su pilastri: taglio Tx e momento torcente Mxy
45
11
Appendice
Le espressioni degli integrali delle funzioni fondamentali che sono necessari
per la costruzione del sistema di equazioni per la soluzione sul contorno,
utilizzando gli integrali base sopra de¯niti, per quanto riguarda il nucleo
Uij¤ , sono:
(k)
Uºº
=
1
(k)
(k)
(4(A00
¡ y¹2¸2 A(k)
¹2 ¸2 E 1(k) +
11 ) ¡ c0 (G0 + y
4¼Dc0
(k)
(0:5 ¡ log(z0))E0 ))
y¹¸2
(k+1)
(k)
(k+1)
(k)
(4(A11
¡x
¹ A11 ) + E 1
¡x
¹ E1 )
4¼D
y¹
(k)
(k)
(G0 + (0:5 ¡ log(z0))E 0 )
= ¡
4¼D
(k)
= Uº¿
1
(k)
(k)
=
(2(A00
¡ A(k)
¹2 ¸2 A(k)
¹2 ¸2 E1(k) +
10 + y
11 ) ¡ c0(G0 ¡ y
4¼Dc0
(k)
Uº¿
=
(k)
Uºw
U¿(k)
º
U¿(k)
¿
(k)
U¿(k)
=
w
(k)
Uwº
=
(k)
Uw¿ =
(k)
U ww
=
(1:5 ¡ log(z0))E0 ))
1
(k+1)
(G(k+1) ¡ x
¹G(k)
¡x
¹E 0(k) ))
0 + (0:5 ¡ log(z0))(E0
4¼D 0
(k)
¡U ºw
(k)
¡U ¿w
1
(k)
(k)
(k)
(k)
¡
(8(G 0 + (1 ¡ log(z0))E0 ) ¡ c0 (G1 ¡ log(z0)F1 ))
8¼D¸2 c0
dove sono stati indicati con
c0
c1
= 1¡º
= 1+º
c2
= 3¡º
Il nucleo delle funzioni Pij¤ diventa:
(k)
Pºº
=
¸2 y¹
(k)
(k)
(k)
(4A(k)
y2 ¸2 (2A(k)
21 ¡ 4A 11 + c2 E1 ¡ 2¹
22 + c0 E 2 ))
4¼
46
¸2
(k+1)
(k)
(k+1)
(k)
(h+1)
(2(A21
¡x
¹A21 ¡ 2A11
+ 2¹
xA11 ) + c0 (E1
¡
4¼
(k+1)
x
¹E 1(k) ) ¡ y¹2 ¸2(4(A(k+1)
¡x
¹A(k)
¡x
¹E 2(k) )))
22
22 ) + 2c0 (E2
¸2 (k)
(A ¡ y¹2 ¸2 A(k)
11 )
2¼ 00
¸2
(k+1)
(k)
(k+1)
(k)
¡x
¹A11 ) + c1 (E 1
¡x
¹E1 )
¡ (4(A11
4¼
(k+1)
(k)
(k+1)
(k)
¡¹
y 2¸2 (4(A22
¡x
¹A22 ) + 2c0 (E 2
¡x
¹E2 )))
¸2y¹
(k)
(k)
(k)
(k)
(k)
(2A21 + 4A11 + c0 E1 ¡ 2¹
y 2¸2 (2A22 + c0 E 2 ))
¡
4¼
y¹¸4 (k+1)
(A11
¡x
¹ A(k)
11 )
2¼
1
((3c0 ¡ 4 + 2c1 log(z0))F0(k) ¡ 2c0 ¸2 y¹2 E 1(k) ¡ 2c1 G(k)
0 )
8¼
c0 y¹¸2 (k+1)
(k)
(E1
¡x
¹E 1 )
4¼
(k)
Pº¿
= ¡
(k)
Pºw
=
P¿(k)
=
º
P¿(k)
=
¿
P¿(k)
=
w
(k)
Pwº
=
(k)
Pw¿
=
(k)
(k)
Pww
=
y¹¸2 E1
2¼
dove le quantitµ
a A(h)
ij sono espresse come:
(k)
A00
(k)
A1h
(k)
A2h
(k)
(k)
(k)
= K 00 + K10 ¡ E 1
(k)
(k)
(k)
= K 0h + 2 K1h ¡ 2 E h+1
(k)
(k)
(k)
(k)
= 4K 0h + K1h¡1 + 8 K1h ¡ 8Eh+1
(k)
(k)
avendo indicato con K 0h
e K 1h
gli integrali che coinvolgono le funzioni
di Bessel modi¯cate
(k)
K0h
(k)
K1h
=
=
Z
x2
x1
x2
Z
x1
K0 [z]xk
dx
z 2h
K1 [z]xk
dx
z 2h+1
che possono essere ricondotti alla somma di due parti integrali di cui la
prima, estesa alla zona di validitµa delle rappresentazioni delle funzioni di
47
Bessel per piccoli argomenti, µe espressa come combinazione degli integrali
analitici prima presentati:
(k)
K 0h
=
h
X
(k)
C1i Eh¡i
i=0
(k)
K 1h
=
(k)
Eh+1
+
Ntssp ³
X
(k)
(k)
C1i Fi¡h + C2i G i¡h
i=h
+
h
X
(k)
C3i E h¡i
i=0
+
Ntssp
X ³
(k)
´
(k)
C3i Fi¡h + C4i Gi¡h
i=h
´
e da una seconda parte, in cui intervengono le rappresentazioni in serie esponenziale, i cui termini possono essere integrati numericamente adottando
accorgimenti atti ad evitare situazioni di quasi singolarit¶a in corrispondenza
del passaggio tra piccoli e grandi argomenti.
Gli integrali dei termini v® necessari per de¯nire l'integrale del termine
di carico, sono:
Vºº
V¿ º
Vwº
1
(0)
(4G(0)
1 ¡ (1 + 4 log(z0))F 1 +
128¼D¸2
(0)
(0)
¸2 y¹2 (8G0 + (2 ¡ 8 log(z0))F0 ))
y¹¸2
(1)
(0)
= ¡
(8(G 0 ¡ x
¹G0 ) +
2
128¼D¸
(1)
(0)
(2 ¡ 8 log(z0))(F0 ¡ x
¹F0 ))
y¹
(0)
(0)
=
(64(G 0 + (0:5 ¡ log(z0))F0 )=c0 ¡
128¼D¸2
(0)
(0)
(4G1 ¡ (1 + 4 log(z0))F1 ));
=
Le espressioni degli integrali delle ulteriori funzioni fondamentali necessari per la costruzione della soluzione all'interno del dominio, utilizzando
ancora gli integrali base prima de¯niti, sono:
(k)
Uººº
(k)
Uº¿
º
(k)
U wºº
1 2
(k)
(k)
(k)
(k)
(k)
¸ y¹(4A21 ¡ 4A11 + c2 E1 ¡ 2y¹2 ¸2 (2A22 + c0 E2 ))
4¼
1 2
(k+1)
=
¸ (2(A(k+1)
¡x
¹ A(k)
+ 2¹
xA(k)
21
21 ¡ 2A 11
11 ) +
4¼
(k+1)
(k)
(k+1)
(k)
(k+1)
(k)
c0 (E1
¡x
¹E 1 ) ¡ y¹2 ¸2 (4(A22
¡ xA
¹ 22 ) + 2c0 (E 2
¡x
¹ E2 )))
¸2 (k)
(k)
=
(A ¡ y¹2 ¸2 A11 )
2¼ 00
= ¡
48
(k)
Uºº¿
=
(k)
Uº¿
¿
=
(k)
U wº¿
=
(k)
U ººw
=
(k)
U º¿w
=
¸2
(k+1)
(k)
(k+1)
(k)
(4(A11
¡x
¹A11 ) + c1 (E 1
¡x
¹E1 ) ¡
4¼
(k+1)
y¹2¸2 (4(A(k+1)
¡x
¹A(k)
¡x
¹E2(k) )))
22
22 ) + 2c0 (E 2
¸2
(k)
(k)
(k)
y¹(2A(k)
y2 ¸2 (2A(k)
21 + 4A 11 + c0E 1 ¡ 2¹
22 + c0 E2 ))
4¼
y¹ 4
¸ (a11(k+1) ¡ x
¹ A(k)
11 )
2¼
1
((3c0 ¡ 4 + 2c1)F0(k) ¡ 2c0 ¸2 y¹2 E 1(k) ¡ 2c1 G (k)
0 )
8¼
c0 y¹¸2 (k+1)
(k)
(E1
¡x
¹ E1 )
4¼
(k)
(k)
Uwºw
U¿(k)
¿º
y¹¸2 E 1
= ¡
2¼
(k)
(k)
(k)
2
= ¡¹
y ¸ (4A11 + (c1 ¡ 2c0 )E1 ¡ 4A21 +
(k)
U¿(k)
¿¿
=
U ¿(k)
¿w
=
(k)
U w¿
º
=
(k)
U w¿
¿
=
(k)
Uw¿w
=
(k)
Pººº
=
(k)
Pºº¿
(k)
Pº¿
º
(k)
Pº¿
¿
(k)
¸2 y¹2 (4A22 + 2c0 E 2 ))c2=2
¸2
(k)
(k+1)
(k+1)
(k)
(¹
x 4A11 ¡ 4A11
+ c1 (E 1
¡x
¹E1 )
4¼
(k+1)
+¸2 y¹2 (4(A(k+1)
¡x
¹A(k)
¡x
¹E 2(k) )))
22
22 ) + 2c0(E2
1
(k)
(k)
(k)
(¡c0 E 0 ¡ 2c1 G0 + 2c0 ¸2 y¹2 E1 )
4¼
(k)
Uwº¿
¸2 (k)
(k)
(A ¡ A10
+ ¸2y¹2 A(k)
11 )
2¼ 00
¸2 (k+1)
(k)
(E
¡ xE
¹ 1 )
2¼ 1
¸2 Dc0
(k)
(k)
(k)
(4(A21 ¡ A11 ) + c3 E1 ¡
4¼
(k)
(k)
(k)
(k)
4¸2 y¹2 (A32 + c2 E 2 ) + 4¸4 y¹4 (A33 + 2c0 E3 ))
y¹¸4 Dc0
(k+1)
(k)
(k+1)
(k)
= ¡
(¡2(A32
¡x
¹ A32 + c2 E2
¡x
¹c2 E 2 ) +
4¼
(k+1)
4¸2 y¹2 (A(k+1)
¡x
¹A(k)
¡ x2c
¹ 0 E3(k) ))
33
33 + 2c0 E 3
(k)
= Pºº¿
¸2 Dc0
(k)
(k)
(k)
(k)
(k)
=
(4(A11 ¡ A21 ) ¡ c0 E1 + A10 ¡ 2A00 +
4¼
49
(k)
(k)
(k)
(k)
4¸2 y¹2 (A32 + 2c0 E 2 ) ¡ 4¸4 y¹4 (A33 + 2c0 E3 ))
¸2 Dc0
(k)
(k)
(k)
(k)
(k)
(4(A11 ¡ A21 ) ¡ c0 E1 + 4¸2y¹2 (A32 + 2c0 E 2 ) ¡
4¼
(k)
(k)
4¸4 y¹4 (A33 + 2c0 E 3 ))
y¹¸4 Dc0
(k+1)
(k)
(k+1)
(k)
¡
((A32
¡x
¹A32 + (4c0 ¡ c2 )(E2
¡x
¹ E2 )) ¡
2¼
(k+1)
2¸2 y¹2 (A(k+1)
¡x
¹A(k)
¡ x2c
¹ 0 E3(k) ))
33
33 + 2c0 E 3
y¹¸4 Dc0
(k)
2 2 (k)
¡
(2(A21
¡ A(k)
¹ A22 )
11 ) ¡ 2¸ y
4¼
¸4 Dc0 (k+1)
(k)
(k+1)
(k)
(k+1)
(k)
(A21
¡x
¹ A21 ¡ 2A11
+ 2¹
xA11 ¡ 2¸2 y¹2(A22
¡x
¹A22 ))
4¼
(k)
¡Pººw
P¿(k)
¿º
=
P¿(k)
¿¿
=
(k)
Pººw
=
(k)
Pº¿w
=
P¿(k)
¿w
=
(k)
Pwºº
(k)
Pwº¿
(k)
= ¡Pººw
(k)
Pw¿
º
(k)
Pw¿
¿
(k)
Pwºw
(k)
Pw¿w
(k)
= ¡Pº¿w
¸4 Dc0
2 2
=
(¡2A(k+1)
+ 2¹
xA(k)
¹ (A(k+1)
¡x
¹A(k)
11
11 + 2¸ y
22
22 ))
4¼
y¹¸4 Dc0
(k)
(k)
(k)
= ¡
(2A11 + A21 ¡ 2¸2 y¹2A22 )
4¼
¸4 Dc0 (k)
2 2
(E1 + A(k)
¹ (2E2(k) + A(k)
=
00 ¡ ¸ y
11 ))
4¼
y¹¸6 Dc0
(k+1)
(k)
(k+1)
(k)
=
(2(E2
¡x
¹E 2 ) + A11
¡x
¹A11 )
4¼
avendo indicato con
c0
c1
= 1¡º
= 1+º
c2
c3
= 3¡º
= 3+º
Gli ulteriori integrali dei termini necessari per de¯nire l'integrale del
termine di carico, sono:
y¹
Uººº
(0)
(0)
((8 + c0 )E00 ¡ 4c3 G0 ¡ 4c0 ¸2 y¹2 E 1 ) ¡ º
64¼
(c0 ¸2 )
c0
U º¿º
(3(E 0(1) ¡ x
¹E 0(0) ) ¡ 4¸2 y¹2 (E 1(1) ¡ x
¹E 1(0) ) ¡ 4(G(1)
¹ G(0)
0 ¡x
0 )) ¡ º
64¼
c0 ¸2
Wºº
= ¡
Wº¿
=
50
W¿ ¿
Wwn
Wwt
y¹
U¿ ¿ º
(0)
(0)
(0)
((8 ¡ 9c0 )E0 ¡ 4(4 ¡ 3c0 )G0 + 4c0 ¸2y¹2 E1 ) ¡ º
64¼
c0 ¸2
1
Uwºº
(0)
(0)
=
(2G 00 ¡ E0 + 2¸2 y¹2 E 1 ) ¡ º
8¼
c0 ¸2
Uw¿ º
y¹¸2 (1)
(0)
(E 1 ¡ x
¹E 1 ) ¡ º
= ¡
4¼
c0 ¸2
= ¡
(h)
in cui i termini A3j sono ulterioni combinazioni delle funzioni di Bessel
e ugualmente a quanto si veri¯ca con i termini di ordine piµ
u basso, sono
rappresentabili in maniera numerica per la rappresentazione per grandi argomenti mentre per quella a piccoli argomenti si rappresenta attraverso una
espressione del tipo:
(k)
A3h
(k)
(k)
(k)
(k)
= K 0(k)
h+1 + 24K 0h + 8K1 h¡1 + 48K 1h ¡ 48Eh+1
References
[1] G.R. Kirchho®, Ueber das Gleichgewicht und die Bewegung einer elastischen
Scheibe, J.reine angew.Math.(Crelle), 1850.
[2] E. Reissner,On the Theory of Bending of Elastic Plates, J. of Math,. Phys.,
23, 1944.
[3] E. Reissner,The E®ect of Transverse Shear deformation on the Bending of
Elastic Plates, J.of Applied Mech.. Phys., 12, 1945.
[4] R.D. Mindlin,In°uence of Rotatory Inertia and Shear on Flexural Motions of
Isotropic Elastic Plates, J. of Applied Mechanics,1951.
[5] B. Haggblad and K.J.Bathe, Speci¯cation of boundary conditions for ReissnerMindlin plate bending ¯nite elements, Int. j. num. method. eng.,30,981-1011
(1990).
[6] O.C. Zienkiewicz, Z. Xu, L.F. Zeng, A. Samuelsson and N.E. Wiberg,Linked
interpolation for Reissner-Mindlin plate elements: part I - A Simple Quadrilateral,I nt. j. num. methods eng., 36, 3043-3056(1993).
[7] R.L. Taylor and F. Auricchio,Linked interpolation for Reissner-Mindlin plate
elements: part II - a simple triangle,Int. j. num. methods eng., 36, 30573066(1993).
51
[8] F. Hartmann and R. Zotemantel, The direct boundary element method in plate
bending, Int. j. numer. methods eng., 23, 2049{2069(1986).
[9] F. Hartmann, Introduction to Boundary Elements. Theory and Applications,
Springer-Verlag, Berlin, 1989.
[10] F. Vander Weeen,Application of the boundary integral equation methods to
Reissner's plate model, Int. j. numer. methods eng.,18, 1{10,1982.
[11] J.T. Katsikadelis A.J. Yotis, A new boundary element solution of tkick plates
modelled by Reissner's theory, Engeneering Analysis whith Boundary Element,
Vol. 12, pp. 65{74, 1993.
[12] V.J. Karaman and J.C.F. Telles, On boundary element for Reissner's plate
theory, Engineering Analysis, Vol. 5, No. 1, pp. 21{27, 1988.
[13] D.N. Arnold and R.S. Falk, The boundary layer for Reissner-Mindlin model ,
SIAM j: Math. Anal.,21,281-312 (1990).
[14] H. HÄormander,Linear Partial Di®erential Operators., Springer, Berlin, 1963.
[15] M. Aristodemo, A high-continuity ¯nite element model for two dimensional
elastic problems, Computer.& Structures, 21, 987-993(1985).
[16] M. Abramowitz and I. A. Stegun,Eds., Handbook of Mathematical Function,
Dover, New York (1972).
[17] J.C.F. Telles,A self-adaptive co-ordinate trasformation for e±cient numerical
evaluation of general boundary element integrals, Int. j. numer. methods eng.,
24, 959{973(1987).
[18] Q. Huang and T.A. Cruse, Some Notes On Singular Integral Techiniques In
Boundary Element Analysis, Int. J. Numer. Method. Eng.,36, 2643 { 2659
(1993).
[19] S.P. Timoshenko, S. Woinowsky-Krieger, Theory of Plates and Shells, McGraw Hill, 1959.
[20] A. Razzaque,Program for triangular bending elements with derivative smoothing, Int. j. numer. methods eng., 6, 333{345(1973).
[21] F.G. Yuan and R.E. Miller,A cubic triangular ¯nite element for °at plates
with shear, Int. j. numer. eng., 28, 109-126(1989).
[22] M. Aristodemo D. Archinµ
a, Boundary elements analisys of Kirchho® plate,
Report 95/1 D.A.S.Te.C. Universitµ
a di Reggio Calabria.
52
[23] S. Lopez and S. Fortino, Soluzione adattativa multigrid di problemi bidimensionali elastici, Department of Structures, Technical Report n.164, University
of Calabria, Arcavacata di Rende, Italy.
[24] V.J. Karaman and J.C.F. Telles,On boundary element for Reissner's plate
theory, Engineering Analysis,Vol. 5,No. 1,pp. 21{27,1988.
53

Download Report