Analisi di telai piani con tecniche FEM:
Maple® e il package <beam2D_FEMpack>
Roberto Lot ©2000
Dipartimento di Ingegneria Meccanica dell'Università degli Studi di Padova
<beam2D_FEMpack> è una raccolta di procedure scritte in Maple che permette di studiare i telai piani con il metodo degli elementi finiti (FEM).
Illustriamo il funzionamento del software attraverso un esempio, analizzando il telaio mostrato nella figura
Innanzitutto carichiamo il pacchetto
> restart: read "beam2D_FEMpack.m";
Risolveremo il nostro problema affrontando una ad una le seguenti fasi:
La prima
operazione da effettuare è la definizione dei nodi della struttura, specificando
la posizione di ciascuno, ossia le coordinate X e Z.
In Maple la definizione dei nodi viene fatta con una table, ciascun
elemento è un nodo così costituito:
nodenumber = [ X coordinate , Y coordinate ]
>
node := table ( [1= [0.,1.] , 2= [+1.6,+0.4] , ,
seq( seq( 10*r+c = [ 0.8*(c-1) , 0.15*(r-1)] , r=1..2) , c=1..3)
]);
![node := TABLE([1 = [0., 1.], 2 = [1.6, .4], 11 = [0...](images/esempio_master2.gif){kind=small}
2. Definizione degli elementi della struttura
Cominciamo definendo le caratteristiche delle sezioni impiegate, utilizzando unità SI
>
E := 210000e6:
modulo di elasticità dell'acciaio
>
rho := 7800: densità
elemento nero, tubo rettangolare
>
bz0 := 0.100: by0:=0.100: s0 := 0.003:
lati del tubo e suo spessore
>
mu0 :=rho*2*(bz0+by0)*s0: massa
lineare
>
EA0 := E*2*(bz0+by0)*s0:
modulo di elasticità assiale
>
EJ0 := (1/12*bz0^3*s0 + 1/4*bz0^2*by0*s0 )*E:
modulo di elasticità flessionale
elementi orizzontali (blu), tubo rettangolare
>
bz1 := 0.050: by1:=0.020: s1 := 0.003:
lati del tubo e suo spessore
>
mu1 :=rho*2*(bz1+by1)*s1:
>
EA1 := E*2*(bz1+by1)*s1:
> EJ1 := (1/12*bz1^3*s1 + 1/4*bz1^2*by1*s1 )*E:
elementi verticali e inclinati (rossi), tubo
rettangolare
>
bz2 := 0.020: by2:=0.020: s2 := 0.003:
lati del tubo e suo spessore
>
mu2 :=rho*2*(bz2+by2)*s2:
>
EA2 := E*2*(bz2+by2)*s2:
>
EJ2 := (1/12*bz2^3*s2 + 1/4*bz2^2*by2*s2 )*E:
Passiamo ora alla definizione degli elementi che
costituiscono il telaio.
Anche in questo caso, ogni
elemento
è descritto tramite una lista contenente:
beam = [ i,j , mu, EA, EJ, l, theta ];
con il seguente significato dei simboli:
i,j = nodi di estremità
mu = massa lineare
EA = modulo di elasticità assiale
EJ = modulo di elasticità flessionale
theta = inclinazione elemento
l = lunghezza
In base alla figura, defiamo l'insime dei suoi elementi
>
beam := table([ [ 1, 2, mu0,EA0,EJ0],
[ 2, 23, mu0,EA0,EJ0], elementi
neri
seq( seq( [ 10*r+c, 10*r+c+1, mu1,EA0,EJ1], c=1..2),r=1..2),
elementi orizzontali
seq( [ 10+c, 21+c, mu2,EA2,EJ2], c=1..2),
elementi diagonali
seq( [ 10+c, 20+c, mu2,EA2,EJ2], c=2..3)
]);
elementi verticali
La definizione dei vincoli del sistema è semplicissima: è sufficiente raccogliere in una lista le coordinate dei nodi solidali al telaio
(In reltà i vincoli possono essere definiti in maniera molto più generale, attraverso una lista di relazioni lineari implicite, non necessariamente indipendenti)
>
impl_constraint :=
[
x1 ,z1,
cerniera nel nodo 1
x21,z21, phi21,
incastro nel nodo 21
x11,z11 ,phi11
];
incastro nel nodo 11
4. Assemblaggio della struttura
E' sufficiente eseguire una chiamata alla la procedura built_structure, la quale automaticamente verifica i nodi, gli elementi e le condizioni di vincolo, separa le coordinate dipendenti da quelle indipendenti e determina il numero di gradi di libertà del sistema.
> built_structure(node,beam,impl_constraint):
Warning, global variables 'NN','nix','coords' were been defined
Sono stati definiti 8 nodi: n = [ X [m] Z [m] ] -------------------------- 1 = [ +0.000 +1.000 ] 2 = [ +1.600 +0.400 ] 11 = [ +0.000 +0.000 ] 12 = [ +0.800 +0.000 ] 13 = [ +1.600 +0.000 ] 21 = [ +0.000 +0.150 ] 22 = [ +0.800 +0.150 ] 23 = [ +1.600 +0.150 ] --------------------------
Warning, global variables 'NE','bix' were been defined
Sono stati definiti 10 elementi di tipo trave: m , ( i, j) = mu[kg/m] EA[N] EJ[Nm2] l[m] theta[rad] -------------------------------------------------------------------- 1 , ( 1 , 2) = 9.4 2.52e+08 2.10e+05 1.709 -0.359 2 , ( 2 ,23) = 9.4 2.52e+08 2.10e+05 0.250 -1.571 3 , (11 ,12) = 3.3 2.52e+08 1.44e+04 0.800 +0.000 4 , (12 ,13) = 3.3 2.52e+08 1.44e+04 0.800 +0.000 5 , (21 ,22) = 3.3 2.52e+08 1.44e+04 0.800 +0.000 6 , (22 ,23) = 3.3 2.52e+08 1.44e+04 0.800 +0.000 7 , (11 ,22) = 1.9 5.04e+07 1.68e+03 0.814 +0.185 8 , (12 ,23) = 1.9 5.04e+07 1.68e+03 0.814 +0.185 9 , (12 ,22) = 1.9 5.04e+07 1.68e+03 0.150 +1.571 10 , (13 ,23) = 1.9 5.04e+07 1.68e+03 0.150 +1.571 --------------------------------------------------------------------
Warning, global variables 'constraint','NC','dip','indip','NDOF' were been defined
La struttura possiede 16 gradi di libertà.
E' possibile disegnare la struttura assemblata utilizzando le opportune subroutine.
Innanzitutto definiamo i colori:
> blue :=COLOUR(RGB,0,0,1): red:=COLOUR(RGB,1,0,0): black:=COLOUR(RGB,0,0,0):
creiamo le etichette dei nodi (i parametri tra parentesi (xoff,zoff) corrispondono all'offset di visualizzazione)
> nlab := node_label(-0.05,-0.05):
creiamo il disegno della struttura con una chiamata alla procedura draw_beams (beams,disp,sc,np) dove:
beams: INPUT lista degli elementi
disp: INPUT spostamenti dei nodi
sc: INPUT scala di amplificazione degli spostamenti
np: INPUT numero di punti con cui rappresentare ciascuna trave
> P_undeformed := draw_beams(beam,0,0,2,red):
e infine plottiamo la struttura
> PLOT(P_undeformed,nlab,SCALING(CONSTRAINED),AXESSTYLE(NONE),THICKNESS(2));
![[Maple Plot]](images/esempio_master10.gif)
Cominciamo con la definizione delle forze esterne applicate alla struttura, specificando per ciascun nodo la componente orizzontale e verticale della forza nonché il momento:
node = [ Fx=horizontal , Fz=vertical, My=torque ];
raccogliamo tutte le forze in una table
> force := table([ 12 = [ 10e3 , 10e3 , 200.0 ] , 13 = [ 10e3 , 10e3 , 0.0 ] ]);
![force := TABLE([12 = [.10e5, .10e5, 200.0], 13 = [....](images/esempio_master11.gif){kind=small}
a questo punto è sufficiente effettuare la chiamata alla procedura:
> static_analysis(beam, force, static_disp):;
dove:
beam: INPUT descrizione degli elementi
force: INPUT forze esterne
static_disp: OUTPUT spostamenti dei nodi
TABLE([1 = [0, 0, .822436218037365594e-2], 2 = [-.187080498859635864e-2, .537291752250592272e-2, -.646329522016223478e-2], 11 = [0, 0, 0], 12 = [.459845758293222998e-4, .487773564099060372e-2, .436063095571934818e-2], 13 = [.877094633005099592e-4, .541307776872617785e-2, .720702090808683963e-3], 21 = [0, 0, 0], 22 = [.112079660452303926e-3, .484974542212747266e-2, .519457235796320048e-2], 23 = [.604855443368294710e-4, .538180807255102604e-2, -.810007060701085045e-2]])
Il risultato è una table in cui per ciascun nodo sono indicate le tre componenti di deformazione
node = [ Dx=horizontal displacement, Dz=vertical displacement, y=rotation ];
E' possibile vedere la configurazione finale utilizzando le apposite procedure grafiche
>
scd:=10:
fattore di amplificazione delle
deformazioni
>
npt:=10:
numero di punti di rappresentazione per
ciascun elemento
>
P_stat_def := draw_beams(beam,static_disp,scd,8,blue): grafico
della struttura deformata
>
scf:=0.00003:
scala delle forze
>
scm:=0.0005:
scala dei momenti
>
PP_force := draw_forces(force,scf,scm,static_disp,scd,red,THICKNESS(1)):
diagramma delle forze
finalmente
>
PLOT(P_stat_def,PP_force,node_label(0.05,0.05),SCALING(CONSTRAINED),AXESSTYLE(NONE),THICKNESS(1));
![[Maple Plot]](images/esempio_master12.gif)
Effettuare il calcolo dei modi di vibrare è molto semplice, è sufficiente chiamare la procedura relativa nella seguente maniera:
>
modal_analysis(beam, myomega, mymode):
dove:
beam:
INPUT descrizione degli elementi
myomega: OUTPUT
pulsazioni dei modi
mymode:
OUTPUT forme modali, ossia spostamenti di ciascun
nodo
Warning, pre-load conditions are not take into account in modal analysis
frequenze proprie
omega( 1) = 4.7213e+04 rad/s = 7514.1812 Hz omega( 2) = 3.0699e+04 rad/s = 4885.9649 Hz omega( 3) = 2.7180e+04 rad/s = 4325.8458 Hz omega( 4) = 2.3147e+04 rad/s = 3683.9492 Hz omega( 5) = 1.8510e+04 rad/s = 2945.9845 Hz omega( 6) = 1.6652e+04 rad/s = 2650.1939 Hz omega( 7) = 9930.3616 rad/s = 1580.4661 Hz omega( 8) = 7559.3493 rad/s = 1203.1078 Hz omega( 9) = 5642.4704 rad/s = 898.0271 Hz omega(10) = 4235.3831 rad/s = 674.0822 Hz omega(11) = 3641.8951 rad/s = 579.6256 Hz omega(12) = 2250.4807 rad/s = 358.1751 Hz omega(13) = 1800.066 rad/s = 286.4894 Hz omega(14) = 323.1182 rad/s = 51.4259 Hz omega(15) = 900.1761 rad/s = 143.2675 Hz omega(16) = 1370.4958 rad/s = 218.1212 Hz
Anche per i modi è facile avere una rappresentazione grafica
>
P_undeformed := draw_beams(beam,0,0,2,red,LINESTYLE(4)):
struttura indeformata
>
scd := 0.5:
fattore di scala per gli spostamenti
>
nn := 14; seleziono il modo
>
myomega[nn]; stampo
la pulaszione
>
P_mode := draw_beams(beam,mymode[nn],scd,6,blue):
struttura deformata
>
PLOT(P_undeformed,P_mode,node_label(0,0),SCALING(CONSTRAINED),AXESSTYLE(NONE),THICKNESS(1));
Ancora più efficace è l'animazione:
>
nfr := 12: numero
dei frame
>
frame:='frame':
>
for i from 1 to nfr do frame[i]:= draw_beams(beam,mymode[nn],evalf(scd*cos(2*Pi*(i-1)/nfr)),6,blue):
od:
> PLOT(ANIMATE(seq([frame[i],P_undeformed],i=1..nfr)),AXESSTYLE(NONE),SCALING(CONSTRAINED));
![[Maple Plot]](images/esempio_master15.gif)
Se volete utilizzare <beam2D_FEMpack>, contattatemi!
ulteriori modi di vibrare della struttura
|
|
|
|
|
|
Altri esempi:
![[Maple Plot]](images/esempio_master18.gif)
![[Maple Plot]](images/esempio_master21.gif)
![[Maple Plot]](images/esempio_master24.gif)
![[Maple Plot]](images/esempio_master27.gif)
