Esercitazioni in Ansys PDF [PDF]

  • 0 0 0
  • Gefällt Ihnen dieses papier und der download? Sie können Ihre eigene PDF-Datei in wenigen Minuten kostenlos online veröffentlichen! Anmelden
Datei wird geladen, bitte warten...
Zitiervorschau

Part I

Esercitazioni in ANSYS Mechanical APDL 1

Introduzione

Tipicamente in ANSYS per input file si intende il file che contiene tutte le istruzioni testuali necessarie all’implementazione e risoluzione del problema in esame. Esso è normalmente composto da 3 parti, corrispondenti a 3 distinti ambienti: • preprocessing, dove viene definito il modello geometrico con tecnica bottomup o top-down (o mista), viene definita la tipologia di elementi da impiegare con le relative proprietà del materiale, viene effettuata la discretizzazione del modello geometrico in elementi finiti e, infine, sono introdotti i vincoli e i carichi. Quest’ultima operazione può anche essere effettuata nella successiva fase di solution. Quando la geometria del modello è particolarmente semplice può essere conveniente costruire il modello agli elementi finiti definendo direttamente i nodi e gli elementi (senza, quindi, passare dalla preventiva modellazione geometrica); • solution: dove viene scelto il solutore da utilizzare, definiti eventualmente i controlli di soluzione e risolto infine il problema; • postprocessing: dove vengono elaborati e visualizzati i risultati ottenuti. Nel seguito, le caratteristiche principale del software verrano esposte attraverso alcuni esempi dimostrativi.

1

2

Esercitazione 1: modellazione geometrica di una giunzione tubolare

Si risolva in ANSYS il seguente problema di modellazione geometrica 3D. Si supponga di voler realizzare tramite il modellatore del programma la struttura in Figura 1:

Figure 1: Cilindri accoppiati. Si definiscono: • d 1 diametro esterno del tubo 1; • d 2 diametro esterno del tubo 2; • t 1 spessore del tubo 1; 2

• t 2 spessore del tubo 2; • h altezza del tubo verticale, che presenta alla quota

h 2

il centro della giunzione

con il tubo 2; • l lunghezza del tubo 2. Come primo passo, è consigliato eliminare dalla memoria qualsiasi comando residuo utilizzando le righe: –––––––––– finish /clear –––––––––– Segue la fase preliminare di definizione dei parametri geometrici. –––––––––– d1 = 0.1 t1= 0.005 d2 = 0.07 t2 = 0.004 h = 0.2 l = 0.15 –––––––––– In ANSYS non è obbligatorio specificare le unità di misura dei dati numerici. Sta all’utente mantenere la coerenza dei dati inseriti, in maniera da ritrovare la stessa coerenza nel risultato dell’analisi. Si può entrare ora in ambiente di pre-processing. –––––––––– /prep7 cyl4,0,0,d1/2-t1„d1/2„h –––––––––– 3

Il comando cyl4, permette di modellare cilindri o aree circolari con diverse caratteristiche. La sintassi è la seguente: CYL4,XCENTER,YCENTER,RAD1,THETA1,RAD2,THETA2,DEPTH quindi occorre nell’ordine inserire le coordinata x e y del centro della base, il raggio esterno del cilindro, l’angolo di partenza della rivoluzione (non necessariamente di 360◦ ), il raggio interno se esistente, l’angolo di arrivo della rivoluzione, e l’altezza del cilindro. Sono stati lasciati vuoti i campi relativi ai due angoli, in quanto è richiesta una rivoluzione completa. Nell’ambito della modellazione solida in ANSYS, il sistema di riferimento utilizzato è quello relativo al cosiddetto working plane. Il working plane rappresenta un diverso sistema di riferimento rispetto a quello globale, e può essere spostato nello spazio virtuale a seconda delle esigenze di modellazione. Gli assi del working plane sono wx , wy e wz . Di default tale sistema di riferimento coincide con quello globale cartesiano. Si deve tener presente che per molte funzioni in ANSYS, le coordinate richieste dalla sintassi delle funzioni sono quelle del working plane, infatti nel caso del comando cyl4 utilizzato precedentemente, le coordinate da immettere sono relative al working plane. Ciò significa che bisogna porre particolare attenzione al posizionamento del working plane nello spazio per far sì che la funzione utilizzata produca il risultato desiderato. A questo punto, generato il cilindro 1, per generare il secondo tramite cyl4 bisogna traslare e ruotare il working plane. –––––––––– wprota„90

!rotazione working plane di 90◦ attorno all’asse x

wpoffs„h/2

!traslazione working plane di

h 2

–––––––––– Lo spostamento del working plane è stato diviso in due step successivi, dati rispettivamente una rotazione prodotta dal comando wprota e da una traslazione data da wpoffs.

4

La sintassi del primo comando è: WPROTA,THXY,THYZ,THZX che richiede nell’ordine l’angolo di rotazione desiderate del working plane attorno agli assi z,x,y. La sintassi del secondo è: WPOFFS,XOFF,YOFF,ZOFF che richiede nell’ordine gli offsets di spostamento lungo gli assi x, y e z. Alternativamente, invece di traslare il working plane, è possibile assegnare nuove coordinate per l’origine della geometria. Ciò comunque comporta che venga effettuata la rotazione del sistema di riferimento. Si procede quindi con la generazione del secondo cilindro: –––––––––– cyl4,0,0,d2/2-t2„d2/2„l+d1/2 –––––––––– Infine, si conclude la modellazione eliminando le parti in eccesso e rifinendo la struttura. La prima operazione da effettuare è rendere solidali i cilindri, tramite il comando vovlap. Esso richiede di indicare i numeri identificativi dei volumi da rendere solidali. Si ricordi a tal proposito che la numerazione delle entità numeriche è automatica e progressiva, per cui il comando sarà in questo caso: –––––––––– vovlap,1,2 –––––––––– Il comando vlist permette in ogni caso di mostrare la numerazione dei volumi. L’unione tra le due entità ha prodotto nuovi volumi. A questo punto si procede all’eliminazione dei volumi da cancellare tramite il comando vdele: VDELE,NV1,NV2,NINC,KSWP

5

La sintassi richiede che vengano indicati i volumi da cancellare passando dal volume nv1 al volume nv2 con passo ninc. Il campo kswp permette di specificare se si vogliono o meno eliminare, oltre ai volumi, anche le strutture geometriche gerarchicamente inferiori (ma non condivise con altri volumi): 0 per non cancellarle, 1 per cancellarle. –––––––––– vdele,3,5,2,1 –––––––––– La geometria è realizzata.

Figure 2: Geometria dei cilindri accoppiati. File input completo. –––––––––– finish /clear !PARAMETRIZZAZIONE 6

!DEFINIZIONE PARAMETRI GEOMETRICI d1 = 0.1 t1= 0.005 d2 = 0.07 t2 = 0.004 h = 0.2 l = 0.15 !PRE-PROCESSING /prep7 !GENERAZIONE TUBO PRINCIPALE cyl4,0,0,d1/2-t1„d1/2„h !SPOSTAMENTO WORKING PLANE wprota„90 wpoffs„h/2 !GENERAZIONE TUBO SECONDARIO cyl4,0,0,d2/2-t2„d2/2„l+d1/2 !UNIONE VOLUMI vovlap,1,2 !CANCELLAZIONE VOLUMI ECCEDENTI vdele,3,5,2,1 –––––––––– NB: Qualsiasi riga di codice preceduta da un punto esclamativo è trascurata dal compilatore, ciò permette di inserire commenti nello script.

7

3

Esercitazione 2: modellazione geometrica di una tazza

Si risolva in ANSYS il seguente problema di modellazione geometrica 3D. Si supponga di voler realizzare tramite il modellatore del programma la struttura in Figura 3:

Figure 3: Tazza Si definiscono: - d 1 diametro esterno della parte cilindrica; - d 2 diametro interno della parte cilindrica; - h 1 altezza della parte cilindrica; 8

- b 1 diametro esterno all’estremità superiore della parte tronco-conica; - b 2 diametro interno all’estremità superiore della parte tronco-conica; - h 1 altezza della parte tronco-conica; - t diametro del manico; - r ingombro del manico; - h 1 /2 altezza del centro del manico rispetto al piano d’appoggio . Come visto in precedenza, il primo passo da effettuare è la pulizia della memoria: –––––––––– finish /clear –––––––––– Si passa ora alla definizione dei parametri geometrici: –––––––––– d1 = 0.08 d2 = 0.06 b1 = 0.06 b2 = 0.04 h1 = 0.12 h2 = 0.03 t = 0.015 r = 0.05 h = 0.005 –––––––––– Dopo essere entranti nell’ambiente di preprocessing, occorre, tenuto conto del posizionamento di default del working plane, generare la parte cilindrica della tazza: –––––––––– /prep7 cyl4,0,0,d2/2„d1/2„h1

9

–––––––––– per la quale si ricorre al comando cyl4 utilizzato già in precedenza. A questo punto si può procedere con la generazione della parte conica con il comando cone, generando due volumi tronco-conici e successivamente tramite l’operazione booleana di sottrazione, eliminare il volume del tronco più piccolo. La sintassi del comando è: CONE,RBOT,RTOP,Z1,Z2,THETA1,THETA2 che richiede nell’ordine l’indicazione del raggio della base inferiore, del raggio della base superiore, le quote delle due basi e gli eventuali angoli di inizio e fine della rivoluzione. Quindi, per il tronco di cono più grande: –––––––––– cone,d1/2,b1/2,h1,h1+h2 –––––––––– mentre per la parte tronco-conica interna: –––––––––– cone,d2/2,b2/2,h1,h1+h2 –––––––––– A questo punto, si procede con l’operazione di sottrazione. Si ricordi che con il comando vlist si ha l’opportunità di visualizzare la numerazione dei volumi. Il comando per l’operazione di sottrazione è vsbv: VSBV,NV1,NV2,SEPO,KEEP1,KEEP2 che richiede nell’ordine il volume da cui sottrarre, il volume che deve essere sottratto, la generazione di un’area condivisa o di aree coincidenti all’interfaccia tra i due volumi, il mantenimento o l’eliminazione dei due volumi. Il comando sarà in questo caso: –––––––––– vsbv,2,3„,delete

10

–––––––––– Bisogna ora "saldare" il volume cilindrico ed il volume tronco-conico, operazione che garantisce la continuità strutturale per l’analisi agli elementi finiti. Il comando da utilizzare è vglue: –––––––––– vglue,all –––––––––– Il comando vovlap in generale non può essere utilizzato in questo caso. Si lascia come esercizio il tentativo. Si genera ora il fondo: –––––––––– cyl4,0,0,d1/2„„,h vovlap,all –––––––––– Per modellare il manico si utilizzerà il comando che permette di generare un toro. La prima operazione è muovere il working plane: –––––––––– wprota„90 wpoffs,d1/2,h1/2 –––––––––– A questo punto, si genera il volume col comando torus: –––––––––– torus,0,t/2,r-t/2 –––––––––– la cui sintassi è: TORUS,RAD1,RAD2,RAD3,THETA1.THETA2 che richiede il raggio interno (il toro può infatti essere realizzato anche cavo), il raggio esterno, il raggio di rivoluzione, gli angoli di inizio e fine rivoluzione.

11

L’operazione ha però generato un toro completo, che va sezionato eliminando le parti superflue. Tramite l’input vsel,r,p si visualizza la numerazione dei volumi. Tramite le due righe di codice seguenti, l’operazione è portata a termine: –––––––––– vovlap,1,6 vdele,8„,1 ––––––––––

Figure 4: Tazza realizzata con Ansys Mechanical APDL. File input completo. –––––––––– finish /clear !DEFINIZIONE PARAMETRI GEOMETRICI d1 = 0.08 d2 = 0.06 12

b1 = 0.06 b2 = 0.04 h1 = 0.12 h2 = 0.03 t = 0.015 r = 0.05 h = 0.005 /prep7 !GENERAZIONE PARTE CILINDRICA cyl4,0,0,d2/2„d1/2„h1 !GENERAZIONE PARTE TRONCO-CONICA ESTERNA cone,d1/2,b1/2,h1,h1+h2 !GENERAZIONE PARTE TRONCO-CONICA INTERNA cone,d2/2,b2/2,h1,h1+h2 !SOTTRAZIONE VOLUMI vsbv,2,3„,delete !FUSIONE VOLUMI vglue,all !GENERAZIONE FONDO cyl4,0,0,d1/2„„h vovlap,all !MOVIMENTAZIONE WORKING PLANE wprota„90 wpoffs,d1/2,h1/2 !GENERAZIONE MANICO torus,0,t/2,r-t/2 vovlap,1,6 vdele,8„,1

13

––––––––––

14

4

Esercitazione 3: elemento asta

Si risolva il problema di analisi strutturale mostrato in Figura 5.

Figure 5: Problema strutturale. La sezione trasversale di ciascun elemento è rappresentata in Figura 6

Figure 6: Sezione trasversale dell’asta. Si chiede di determinare, date le forze applicate e i vincoli, la tensione risultante negli elementi. Inizialmente, si procede con la definizione dei parametri geometrici e di carico. –––––––––– 15

l=1.7 h=1.2 d=0.03 t=0.01 E=210e+9 F1=5000 F2=3000 F3=2000 F4=1000 –––––––––– L’unico parametro del materiale di interesse è il modulo di Young, poichè trattandosi di una travatura reticolare, gli elementi possono essere modellati come aste soggette a sforzo normale. Nell’ambiente di pre-processing come prima cosa definiamo la tipologia dell’elemento, il materiale e le proprietà della sezione. –––––––––– /prep7 et,1,link180 mp,ex,1,E sectype,1,link secdata,3.1416*((d/2)**2-(d/2-t)**2) –––––––––– Il comando et ha la seguente sintassi: ET, ITYPE, Ename, KOP1, KOP2, KOP3, KOP4, KOP5, KOP6, INOPR dove nel campo ITYPE occorre inserire un numero identificativo dell’elemento, mentre nel campo Ename viene specificata la tipologia, in questo caso link180. Il comando può anche reggere ulteriori key-options, in questo caso non utilizzate. Il comando mp viene usato per definre le proprietà del materiale, con la seguente

16

sintassi: MP, Lab, MAT, C0, C1, C2, C3, C4 in cui ex è il modulo di Young, 1 è l’identification number del materiale, E è il valore della proprietà. Il comando Sectype: SECTYPE, SECID, Type, Subtype, Name, REFINEKEY richiede di inserire un identification number, seguito dal tipo dell’elemento di cui si sta definendo la sezione. Invece Secdata permette di assegnare il valore dell’area, la sua sintassi sarà: SECDATA, VAL1, VAL2, VAL3, VAL4, VAL5, VAL6, VAL7, VAL8, VAL9, VAL10, VAL11, VAL12 Si prosegue con la generazione dei nodi. Il comando utilizzato è n N,NODE,X,Y,Z che richiede il numero identificativo del nodo e le sue coordinate spaziali. –––––––––– n,1 n,2,l/2,-h n,3,l n,4,3*l/2,-h n,5,2*l n,6,5*l/2,-h n,7,3*l n,8,7*l/2,-h n,9,4*l –––––––––– Per la generazione degli elementi asta, si indicano i nodi di estremità col comando e. Prima di generare gli elementi, è necessario specificarne il tipo, il materiale e il

17

numero di riferimento della sezione, coi comandi type, mat, secnum, che richiedono identification numbers precedentemente definiti. In questo caso tale passaggio può essere omesso in quanto non si è definita più di una variante per ciascuna delle caratteristiche sopraelencate. –––––––––– type,1 mat,1 secnum,1 e,1,2 $ e,1,3 $ e,2,3 $ e,2,4 $ e,3,4 e,3,5 $ e,4,5 $ e,4,6 $ e,5,6 $ e,5,7 e,6,7 $ e,6,8 $ e,7,8 $ e,7,9 $ e,8,9 –––––––––– NB: Il simbolo $ serve per inserire più comandi sulla stessa riga di codice.

Figure 7: Travatura reticolare. A questo punto si possono specificare i vincoli e i carichi. Il comando d permette di indicare il nodo su cui applicare i vincoli, seguiti dai gradi di libertà vincolati (in questo caso, all). Per quanto riguarda i carichi, si usa il comando f, che richiede di specificare il nodo 18

di applicazione del carico, la direzione e l’entità del carico. Per i nodi su cui sono applicati uguali carichi si procede dapprima selezionando i nodi di interesse tramite il comando nsel, successivamente si applica il carico e poi si annulla la selezione effettuate con nall. –––––––––– d,1,all d,9,all f,5,fy,-F1 nsel,r,node„3,7,4 f,all,fy,-F2 nall nsel,r,node„2,8,6 f,all,fy,-F3 nall nsel,r,node„4,6,2 f,all,fy,-F4 nall –––––––––– A questo punto si può risolvere il problema e visualizzare i risultati. Il comando /dscale è utile per specificare la scale di visualizzazione. –––––––––– /solu antype,static solve /post1 /dscale,1,1 –––––––––– Nel caso in cui si volesse visualizzare la deformata rispetto all’indeformata si può

19

usare il comando pldisp,1 ; mentre plvect,u permette di visualizzare gli spostamenti nodali come vettori. –––––––––– pldisp,1 plvect,u ––––––––––

Figure 8: Spostamenti nodali. Il comando Etable permette di creare una tabella, di nome sxx1, in cui caricare i valori della tensione (ls) sul nodo i-esimo. Nella stessa maniera si può creare una tabella, di nome sxx2, con i valori della tensione sul nodo j-esimo. Plls plotta a video tali valori. –––––––––– etable,sxx1,ls,1 etable,sxx2,ls,2 plls,sxx1,sxx2 –––––––––– Di seguito il comando etable è riutilizzato per creare una tabella per lo sforzo normale. Pletab,N visualizza a video lo sforzo normale, mentre pretab,N stampa a video i risultati ottenuti. 20

Figure 9: Diagramma delle tensioni. –––––––––– etable,N,smisc,1 pletab,N pretab,N ––––––––––

Figure 10: Valori dello sforzo normale. File input completo. 21

–––––––––– finish /clear !PARAMETRIZZAZIONE !DEFINIZIONE PARAMETRI GEOMETRICI l=1.7 h=1.2 d=0.03 t=0.01 !DEFINIZIONE PARAMETRI DI CARICO E=210e+9 F1=5000 F2=3000 F3=2000 F4=1000 !PRE-PROCESSING /prep7 et,1,link180 mp,ex,1,E sectype,1,link secdata,3.1416*((d/2)**2-(d/2-t)**2) n,1 n,2,l/2,-h n,3,l n,4,3*l/2,-h n,5,2*l n,6,5*l/2,-h n,7,3*l

22

n,8,7*l/2,-h n,9,4*l !GENERAZIONE ELEMENTI type,1 mat,1 secnum,1 e,1,2 $ e,1,3 $ e,2,3 $ e,2,4 e,3,4 $ e,3,5 $ e,4,5 $ e,4,6 e,5,6 $ e,5,7 $ e,6,7 $ e,6,8 e,7,8 $ e,7,9 $ e,8,9 !DEFINIZIONE VINCOLI d,1,all d,9,all !DEFINIZIONE CARICHI f,5,fy,-F1 nsel,r,node„3,7,4 f,all,fy,-F2 nall nsel,r,node„2,8,6 f,all,fy,-F3 nall nsel,r,node„4,6,2 f,all,fy,-F4 nall !SOLVING /solu antype,static solve

23

!POST-PROCESSING /post1 /dscale,1,1 etable,sxx1,ls,1 etable,sxx2,ls,2 plls,sxx1,sxx2 etable,N,smisc,1 pletab,N pretab,N ––––––––––

24

5

Esercitazione 4: elemento trave

In ANSYS vi è la possibilità di utilizzare due tipologie di elemento trave: una basata sulla teoria di Eulero-Bernoulli (beam4 e beam44 ), l’altra sulla teoria di Timoshenko (beam188 ). In particolare, con il beam44 e il beam188 è possibile definire la sezione sfruttando una apposita libreria. Si risolva il problema di analisi strutturale mostrato in Figura 11

Figure 11: Schema della struttura. Gli elementi (1) e (2) sono travi a sezione rettangolare, l’elemento (3) è una trave a sezione quadrata, mentre l’elemento (4) è un asta. Le travi (1) e (2) sono disconnesse a momento. Analoga considerazione vale per l’elemento (4), tuttavia essendo schematizzato come un’asta il grado di libertà di rotazione è assente. Il file input inizia con la definizione dei parametri geometrici della struttura, del materiale e di carico: –––––––––– l=15

!lunghezza

h=20

!altezza pilone 25

b=7

!larghezza sezione carreggiata

a=1.2

!altezza sezione carreggiata

s=2

!lato sezione pilone

d=0.4

!diametro sezione strallo

Ec= 25e+9

!modulo Young cemento armato

Ea=210e+9

!modulo Young acciaio

ni=0.3

!coefficiente Poisson acciaio

Q=120e+5

! carico distribuito su campate

–––––––––– Come si può notare è stato necessario definire il coefficiente di Poisson poichè si è deciso di schematizzare il problema con elementi trave che utilizzano la teoria di Timoshenko. Succesivamente, si definiscono i keypoints 1, 2, 3, 4 e 5 e le linee. I keypoints 2 e 3 sono coincidenti, in questo modo è possibile accoppiarli lasciando disaccoppiati i gradi di libertà di rotazione. –––––––––– /prep7 k,1 k,2,l/2 k,3,l/2 k,4,l/2+L k,5,3*l/2,h l,1,2 l,3,4 l,4,5 l,3,5 –––––––––– Per l’elemento beam188 è possibile definire l’ordine delle funzioni di forma attra-

26

verso il keyopt(3): • Keyopt(3)=0, l’elemento beam188 è basato su funzioni di forma lineari; • Keyopt(3)=2, l’elemento beam 188 ha un solo nodo interno per lo schema di interpolazione ed è basato su funzioni di forma quadratica; • Keyopt(3)=3, l’elemento beam 188 ha due nodi interni ed è basato su funzioni di forma cubica. –––––––––– et,1,beam188„,3 et,2,link180 mp,ex,1,Ec mp,ex,2,Ea mp,nuxy,2,ni sectype,1,beam,rect secdata,b,a sectype,2,beam,rect secdata,s,s sectype,3,link secdata,3.1416*(d/2)**2 –––––––––– Segue la generazione della mesh. –––––––––– esize,l/20 Type,1 $ mat,1 $ secnum,1 Lmesh,1,2 Type,1 $ mat,1 $ secnum,2 Lmesh,3 27

–––––––––– Le aste non conviene discretizzarle in più elementi, ma utilizzarne uno solo. –––––––––– esize„1 Type,2 $ mat,2 $ secnum,3 Lmesh,4 –––––––––– A questo punto si inseriscono i vincoli direttamente sul modello FEM. Si selezionano i keypoints e i nodi ad essi associati e si assegnano gli spostamenti con il comando d, la cui sintassi è: D, Node, Lab, VALUE, VALUE2, NEND, NINC, Lab2, Lab3, Lab4, Lab5, Lab6 indicando "all" si assegnano gli spostamenti ai nodi precedentemente selezionati, in "label" si indica il grado di libertà. –––––––––– ksel,r,kp„1 nslk,r d,all,ux„„,rotz,rotx,roty,uz allsel,all ksel,r,kp„4 nslk,r d,all,all allsel,all ksel,r,kp„5 nslk,r d,all,ux„„,rotz,rotx,roty,uz allsel,all ksel,r,kp„2,3 nslk,r

28

cp,next,ux,all $ cp,next,uy,all $ cp,next,uz,all $ cp,next,rotx,all $ cp,next,roty,all allsel,all –––––––––– Successivamente si provede con l’accoppiamento dei keypoint 2 e 3 attraverso il comando: CP, NSET, Lab, NODE1, NODE2, NODE3, ... , NODE17 lasciando disaccoppiato il grado di libertà di rotazione. –––––––––– ksel,r,kp„2,3 nslk,r cp,next,ux,all cp,next,uy,all allsel,all –––––––––– La struttura considerata, pur essendo piana e non avendo carichi fuori dal piano stesso, è definita nello spazio, quindi per evitare la labilità nella direzione uscente dal piano (direzione z ) è opportuno vincolare la struttura in tale direzione. Ciò è stato già fatto vincolando tutti i gradi di libertà per il keypoint 4. Di seguito si assegnano i carichi con il comando sfbeam SFBEAM, Elem, LKEY, Lab, VALI, VALJ, VAL2I, VAL2J, IOFFST, JOFFST, LENRAT Per stabilire su quale faccia agisce il carico si deve assegnare un valore al campo LKEY. Osservando lo schema dell’elemento beam188 riportato in figura, è possibile intuire la faccia su cui agisce il carico distribuito, infatti: • i loadkeys 4 e 5, che giacciono sull’asse x, identificano le estremità della trave; • il loadkey 1 identifica la faccia di normale z ; • il loadkey 2 identifica la faccia di normale y; 29

• al loadkey 3 corrisponde un carico distribuito lungo l’asse della trave.

Figure 12: Schema dell’elemento beam188. In questo caso, essendo il sistema di riferimento locale orientato come il sistema di riferimento globale, si evince che il carico deve agire sulla faccia di normale y per cui il valore del loadkey è 2. –––––––––– lsel,r,line„1,2 esll,r nsle,r,1 sfbeam,all,2,pres,q allsel,all –––––––––– A questo punto si può procedere con la soluzione. –––––––––– /solu Solve –––––––––– La struttura modellata con il software ANSYS è la seguente: Con i seguenti comandi è possibile entrare in ambiente post-processing al fine di visualizzare la deformata della struttura. 30

Figure 13: Struttura modellata con ANSYS. –––––––––– /post1 /dscale,1,1 pldisp,1 –––––––––– File input completo. –––––––––– finish /clear !DEFINIZIONE PARAMETRI GEOMETRICI l=15 h=20 b=7 a=1.2 s=2 d=0.4

31

Figure 14: Struttura deformata. !DEFINIZIONE PARAMETRI MATERIALE Ec= 25e+9 Ea=210e+9 ni=0.3 !DEFINIZIONE PARAMETRI DI CARICO Q=120e+5 /prep7 !GENERAZIONE KEYPOINTS k,1 k,2,l/2 k,3,l/2 k,4,l/2+l k,5,3*l/2,h !GENERAZIONE LINEE CHE COLLEGANO I KEYPOINT l,1,2 32

l,3,4 l,4,5 l,3,5 !DEFINIZIONE ELEMENTO, MATERIALE E SEZIONI et,1,beam188„,3 et,2,link180 mp,ex,1,Ec mp,ex,2,Ea mp,nuxy,2,ni sectype,1,beam,rect secdata,b,a sectype,2,beam,rect secdata,s,s sectype,3,link secdata,3.1416*(d/2)**2 !CREAZIONE MESH esize,l/20 Type,1 $ mat,1 $ secnum,1 Lmesh,1,2 Type,1 $ mat,1 $ secnum,2 Lmesh,3 esize„1 Type,2 $ mat,2 $ secnum,3 Lmesh,4 !INSERIMENTO DEI VINCOLI ksel,r,kp„1 nslk,r d,all,ux„„,rotz,rotx,roty,uz

33

allsel,all ksel,r,kp„4 nslk,r d,all,all allsel,all ksel,r,kp„5 nslk,r d,all,ux„„,rotz,rotx,roty,uz allsel,all ksel,r,kp„2,3 nslk,r cp,next,ux,all cp,next,uy,all cp,next,uz,all cp,next,rotx,all cp,next,roty,all allsel,all !APPLICAZIONE CARICHI lsel,r,line„1,2 esll,r nsle,r,1 sfbeam,all,2,pres,q allsel,all /solu Solve /post1 /dscale,1,1 pldisp,1

34

–––––––––– Si vuole visualizzare a video i diagrammi delle caratteristiche delle sollecitazioni. Per visualizzare il diagramma del momento flettente Mx per l’elemento beam188, si scrive: –––––––––– Esel,r,type„1

!si selezionano tutti gli elementi con identification number pari a

1 Etable,Mz1,smisc,3

!si scrive 3 per ottenere il valore sul nodo i

Etable,Mz2,smisc,16 !si scrive 16 per ottenere il valore sul nodo j Plls,Mz1,Mz2 ––––––––––

Figure 15: Diagramma del momento flettente. Dalla Figura 15 si nota come i valori negativi del momento flettente siano riportanti in basso, mentre i valori positivi verso l’alto. La convenzione utilizzata dal 35

software ANSYS può essere invertita. Per estrapolare lo sforzo normale si scrive: –––––––––– Etable,N1,smisc,1 Etable,N2,smisc,14 Plls,N1,N2 ––––––––––

Figure 16: Diagramma dello sforzo normale per la trave. Per estrapolare i valori di sforzo normale dell’asta: –––––––––– Allsel,all Esel,r,type„2 Etable,Na1,smisc,1 Pletable,Na1 ––––––––––

36

6

Esercitazione 5: esempio di capannone

Si consideri il problema strutturale mostrato in Figura 17.

Figure 17: Problema strutturale. Definito il sistema di riferimento globale come indicato le sezioni dei ritti e dei correnti superiori, con la relativa orientazione, sono evidenziate nella Figura 18. I ritti sono quindi dei profili HEA100, mentre i correnti superiori sono degli UPN100. Definizione dei parametri geometrici: –––––––––– l=8 b=5 37

Figure 18: Profili dei ritti: HEA100 e UPN100 h1=4 h2=2 w1h=0.1 w2h=0.1 w3h=0.1 t1h=0.01 t2h=0.01 t3h=0.006 w1u=0.05 w2u=0.05 w3u=0.1 t1u=0.0085 t2u=0.0085 t3u=0.006 –––––––––– Per ciò che riguarda le sezioni, le indicazioni sulle dimensioni delle aree fanno riferimento a tabelle normate. In particolare, per il profilo HEA100, w1h e w2h costituiscono le larghezze delle due ali, mentre w3h è l’altezza della sezione; t1h, t2h e t3h sono rispettivamente gli spessori delle due ali e lo spessore dell’anima. Per 38

quanto riguarda l’UPN100 invece, w1u e w2u costituiscono le larghezze delle due ali, w3u è l’altezza della sezione, e t1u, t2u e t3u sono ancora rispettivamente gli spessori delle due ali e lo spessore dell’anima. Si può procedere con la definizione dei parametri del materiale e di carico. –––––––––– E=210e+9 ni=0.3 rho=7850 fx1=5000 fx2=10000 fy1=800 fy2=400 –––––––––– Modello geometrico. Il primo passaggio è la defininizione dei keypoint, seguita da quella delle linee, secondo la numerazione riportata nella Figura 19: –––––––––– /prep7 k,1 k,2,b k,3,b,h1 k,4,b/2,h1+h2 k,5„h1 kgen,2,1,5,1„,-l –––––––––– L’utilizzo del comando kgen, la cui sintassi è: KGEN, ITIME, NP1, NP2, NINC, DX, DY, DZ, KINC, NOELEM, IMOVE ha permesso di copiare i nodi da 1 a 5 della facciata principale ad una distanza -l

39

Figure 19: Numerazione dei keypoints e delle linee. sull’asse Z. –––––––––– l,1,5 $ l,2,3 $ l,5,3 $ l,3,4 $ l,5,4 l,3,8 $ l,5,10 $ l,4,9 l,6,10 $ l,7,8 $ l,10,8 $ l,8,9 $ l,10,9 –––––––––– Si deve ora trasformare il modello geometrico in un modello di calcolo agli elementi finiti. Si prosegue quindi con la definizione di elementi, materiali e sezioni. –––––––––– et,1,beam188„,3 mp,ex,1,E mp,nuxy,1,ni mp,dens,1,rho

40

sectype,1,beam,i secdata,w1h,w2h,w3h,t1h,t2h,t3h sectype,2,beam,chan secdata,w1u,w2u,w3u,t1u,t2u,t3u –––––––––– Si ricordi che l’elemento beam188, per default, utilizza shape functions lineari. Per utilizzare shape functions di tipo cubico, bisogna assegnare al keyopt(3) il valore 3. Per la generazione della mesh si parte dai ritti, che secondo la numerazione data sono le linee l1,l2,l9,l10. Si richiede che venga specificato il numero di elementi desiderato per ogni ritto. Ciò viene effettuato con il comando esize. specificando la dimensione dell’elemento, in modo che esso comporti automaticamente il ricorso ad un certo numero di elementi per ritto. –––––––––– esize,0.2 –––––––––– A questo punto si deve assegnare l’elemento. Per assegnare in questo caso, la tipologia di elemento, il materiale, la tipologia di sezione, ed eventualmente il terzo nodo di orientazione, si ricorre al comando latt, la cui sintassi è: LATT, MAT, REAL, TYPE, —, KB, KE, SECNUM. Occorre assegnare nell’ordine: materiale, eventuali costanti reali, tipo di elemento, keypoint di orientamento (KB) per la sezione del nodo i e keypoint di orientamento per la sezione del nodo j (in modo che sia possibile conferire due orientamenti diversi alle sezioni estreme dell’elemento), numero identificativo della sezione. Mesh ritti: –––––––––– latt,1„1„„1 lmesh,all

41

allsel,all –––––––––– Si supponga quindi inizialmente di non assegnare nessun keypoint di orientazione. Come si può vedere però, il risultato dell’operazione non coincide con il risultato desiderato, a causa di una diversa orientazione delle sezioni rispetto a quella richiesta dalla traccia. Per come sono state costruite le linee di riferimento degli elementi, è necessario definire il keypoint di orientazione.

Figure 20: Nodi di riferimento per l’elemento "beam188". Dal confronto tra il sistema di riferimento locale e quello globale riportati nella Figura 20, si deduce l’orientazione corretta della sezione, espressa dalle righe di codice sostitutive seguenti: –––––––––– k,100,2*b,h1/2 latt,1„1„100„1 lmesh,1,2 k,101,2*b,h1/2,-l latt,1„1„101„1 lmesh,9,10 –––––––––– 42

Ripeto il procedimento senza specificazione dei keypoints anche per i longheroni, ovvero i correnti superiori longitudinali. –––––––––– latt,1„1„„2 lmesh,6,8,1 –––––––––– Anche in questo caso risulta necessaria la definizione dei keypoints di orientamento. I comandi corretti sono: –––––––––– k,102,-b,h1,-l/2 latt,1„1„102„2 lmesh,6,7 k,103,-b,h1+h2,-l/2 latt,1„1„103„2 lmesh,8 –––––––––– Si ripete infine il procedimento senza specificazione dei keypoints anche per i traversi. –––––––––– latt,1„1„„2 lmesh,3,5,1 –––––––––– Anche in questo caso risulta necessario correggere l’orientazione delle sezioni: –––––––––– latt,1„1„„2 lmesh,4 k,104,b/2,h1,-l/2 latt,1„1„104„2

43

lmesh,3 k,105,b/4,h1+h2/2,-l latt,1„1„105„2 lmesh,5 latt,1„1„„2 lmesh,12 k,106,b/2,h1,-2*l latt,1„1„106„2 lmesh,11 k,107,b/4,h1+h2/2,-2*l latt,1„1„107„2 lmesh,13 –––––––––– A questo punto si procede con la definizione di vincoli e carichi. –––––––––– nsel,r,loc,y,0 d,all,all nall lsel,r,line„6,7$esll,r$nsle,r,all sfbeam,all,2,pres,-fy2 allsel,all lsel,r,line„8$esll,r$nsle,r,all sfbeam,all,2,pres,-fy1 allsel,all lsel,r,line„1,9,8$esll,r$nsle,r,all sfbeam,all,1,pres,-fx1 allsel,all lsel,r,line„2,10,8$esll,r$nsle,r,all

44

sfbeam,all,1,pres,-fx2 allsel,all –––––––––– Ora si procede con il calcolo della soluzione e con la visualizzazione dei risultati. –––––––––– /solu solve /post1 /dscale,1,1 pldisp,1 ––––––––––

Figure 21: Struttura deformata. Con le stringhe di codice: –––––––––– eshape,1,1 plnsol,s,eqv –––––––––– 45

vengono stampati a video i valori delle tensioni equivalenti di Von Mises all’interno degli elementi.

Figure 22: Struttura con i rispettivi valori delle tensioni. Come al solito sarà possibile visualizzare le caratteristiche della sollecitazione. Per i momenti ad esempio: –––––––––– etable,my1,smisc,2 etable,my2,smisc,15 plls,my1,my2 etable,mz1,smisc,3 etable,mz2,smisc,16 plls,mz1,mz2 –––––––––– File input completo. –––––––––– finish

46

Figure 23: Caratteristiche delle sollecitazioni. /clear !DEFINIZIONE PARAMETRI GEOMETRICI (in m) !DEFINIZIONE LUNGHEZZE ELEMENTI TRAVE l=8 b=5 h1=4 h2=2 !DEFINIZIONE PROFILO DI SEZIONE HEA100 w1h=0.1 w2h=0.1 w3h=0.1 t1h=0.01 t2h=0.01 t3h=0.006 47

!DEFIZIONE PROFILO DI SEZIONE UPN100 w1u=0.05 w2u=0.05 w3u=0.1 t1u=0.0085 t2u=0.0085 t3u=0.006 !DEFINIZIONE PARAMETRI DEL MATERIALE E=210e+9 ni=0.3 rho=7850 !DEFINIZIONE PARAMENTRI DI CARICO (in N/m) fx1=5000 fx2=10000 fy1=800 fy2=400 /prep7 !DEFINIZIONE KEYPOINT k,1 k,2,b k,3,b,h1 k,4,b/2,h1+h2 k,5„h1 kgen,2,1,5,1„,-l !DEFINIZIONE LINEE l,1,5 $ l,2,3 $ l,5,3 $ l,3,4 $ l,5,4 l,3,8 $ l,5,10 $ l,4,9 l,6,10 $ l,7,8 $ l,10,8 $ l,8,9 $ l,10,9

48

!DEFINIZIONE ELEMENTI, MATERIALI E SEZIONI et,1,beam188„,3 mp,ex,1,E mp,nuxy,1,ni mp,dens,1,rho sectype,1,beam,i secdata,w1h,w2h,w3h,t1h,t2h,t3h sectype,2,beam,chan secdata,w1u,w2u,w3u,t1u,t2u,t3u esize,0.2 !MESH RITTI k,100,2*b,h1/2 latt,1„1„100„1 lmesh,1,2 k,101,2*b,h1/2,-l latt,1„1„101„1 lmesh,9,10 !MESH LONGHERONI k,102,-b,h1,-l/2 latt,1„1„102„2 lmesh,6,7 k,103,-b,h1+h2,-l/2 latt,1„1„103„2 lmesh,8 !MESH TRAVERSI latt,1„1„„2 lmesh,4 k,104,b/2,h1,-l/2

49

latt,1„1„104„2 lmesh,3 k,105,b/4,h1+h2/2,-l latt,1„1„105„2 lmesh,5 latt,1„1„„2 lmesh,12 k,106,b/2,h1,-2*l latt,1„1„106„2 lmesh,11 k,107,b/4,h1+h2/2,-2*l latt,1„1„107„2 lmesh,13 !VINCOLI nsel,r,loc,y,0 d,all,all nall !CARICHI lsel,r,line„6,7 $ esll,r $ nsle,r,all sfbeam,all,2,pres,-fy2 allsel,all lsel,r,line„8 $ esll,r $ nsle,r,all sfbeam,all,2,pres,-fy1 allsel,all lsel,r,line„1,9,8 $ esll,r $ nsle,r,all sfbeam,all,1,pres,-fx1 allsel,all lsel,r,line„2,10,8 $ esll,r $ nsle,r,all

50

sfbeam,all,1,pres,-fx2 allsel,all /solu solve /post1 /dscale,1,1 pldisp,1 ––––––––––

51

7

Esercitazione 6: elementi piani

Si supponga di dover risolvere il problema piano schizzato in Figura 24.

Figure 24: Piastra rettangolare forata. Trattasi di una piastra forata, sollecitata da una tensione costante σ 0 . Si voglia in particolare determinare le tensioni nel punto A e nel punto B. Data l’evidente doppia simmetria, è possibile condurre lo studio solo su un quarto della struttura, imponendo gli opportuni i vincoli sugli assi di simmetria.

Figure 25: Quarta parte della piastra forata. Si definiscono i parametri del problema, i keypoint e le linee che collegano quest’ultimi: –––––––––– 52

b=1 w=0.6 r=0.1 E=210e+9 ni=0.3 s0=100e+6 /prep7 k,1,r $ k,2,2*r $ k,3,b/2 $ k,4„r $ k,5„2*r $ k,6,2*r,2*r k,7,b/2,2*r $ k,8„w/2 $ k,9,2*r,w/2 $ k,10,b/2,w/2 l,1,2 $ l,2,3 $ l,4,5 l,2,6 $ l,3,7 $ l,5,6 l,6,7 $ l,5,8 $ l,6,9 l,7,10 $ l,8,9 $ l,9,10 –––––––––– Per tracciare l’arco di circonferenza, si ricorre alla funzione larc, la cui sintassi è: LARC, P1, P2, PC, RAD che richiede nell’ordine i punti iniziale (P1) e finale dell’arco (P2), il punto di riferimento della curvatura (PC) e il raggio (RAD). Si definisce il keypoint aggiuntivo 11 in corrispondenza dell’origine. –––––––––– k,11 larc,1,4,11,r ldiv,13,0.5 l,12,6 –––––––––– Il comando ldiv è servito per suddividere l’arco in due parti uguali. La sintassi è: LDIV, NL1, RATIO, PDIV, NDIV, KEEP che richiede l’inserimento del numero della linea da suddividere e del rapporto di

53

suddivisione, oltre che eventualmente del numero desiderato del keypoint generato dalla divisione. Non avendo specificato alcun numero preferenziale per il nuovo keypoint, il software assume il primo numero di keypoint disponibile (il 12). I keypoints 12 e 6 vengono collegati da una linea.

Figure 26: Geometria della piastra con i relativi keypoints. Scrivendo nel command prompt i due comandi riportati di seguito, è possibile visualizzare a video la numerazione delle linee: –––––––––– /pnum,line,1 lplot –––––––––– A questo punto si procede con la realizzazione delle aree: –––––––––– 54

Figure 27: Visualizzazione a video della numerazione delle linee.

55

al,3,6,15,14 al,4,1,13,15 al,4,2,5,7 al,8,6,9,11 al,9,7,10,12 –––––––––– Successivamente, si definiscono le proprietà dell’elemento e del materiale: –––––––––– et,1,plane42„1 Mp,ex,1,E Mp,nuxy,1,ni –––––––––– Per la risoluzione del problema si ricorre all’elemento plane42. In corrispondenza del keyopt(2) si sceglie il tipo di formulazione, inserendo i seguenti valori nel corrispettivo campo: 1. keyopt(2)=0, formulazione con gli elementi incompatibili (correzione del fenomeno dello shear-locking); 2. keyopt(2)=1, formulazione base (nessuna correzione dello shear-locking). Inoltre, attraverso il keyopt(3) si stabilisce il tipo di comportamento: • keyopt(3)=0, plane stress; • keyopt(3)=1, per l’elemento assialsimmetrico; • keyopt(3)=2, plane strain. Generazione mesh: –––––––––– Lsel,r,line„4,6,2 $ lsel,a,line„13,14 $ lsel,a,line„5,11,6 56

Lesize,all„,6 Lsel,all Lsel,r,line„1,3,2 $ lsel,a,line„15 Lesize,all„,6 Lsel,all Lsel,r,line„8,10 Lesize,all„,4 Lsel,all Lsel,r,line„2,12,5 Lesize,all„,8 Lsel,all Type,1 $ mat,1 Amesh,all –––––––––– Il comando lsel è servito per selezionare le linee interessate. La sintassi è: LSEL, Type, Item, Comp, VMIN, VMAX, VINC, KSWP. In seguito, si sono suddivise le linee in un numero opportuno, in modo tale da avere una mesh regolare anche in corrispondenza dell’area in cui è presente il foro. La sintassi che permette di realizzare la suddivisione delle linee è: LESIZE, NL1, SIZE, ANGSIZ, NDIV, SPACE, KFORC, LAYER1, LAYER2, KYNDIV, che richiede in ordine, la selezione delle linee (NL1), la dimensione di suddivisione (SIZE ) o l’angolo di suddivisione di una linea curva (ANGSIZE) o il numero di suddivisione delle linee (NDIV). Per i restanti campi del comando LSEL si rimanda all’help del software. La sintassi del comando Amesh, che permette di generare la mesh, è AMESH, NA1, NA2, NINC Gli elementi ottenuti sono tutti elementi quadrilateri.

57

Figure 28: Mesh dell’elemento piano bidimensionale. –––––––––– nsel,r,loc,x,0 Dsym,symm,x Nall nsel,r,loc,y,0 Dsym,symm,y Nall –––––––––– Si sono introdotti i vincoli di simmetria. Con la sintassi: NSEL, Type, Item, Comp, VMIN, VMAX, VINC, KABS, si sono selezionati tutti i nodi che giacciono sulla retta di equazione x=0. Infine, utilizzando il comando Dsym il software imposta in automatico il vincolo opportuno, la sintassi è:

58

DSYM, Lab, Normal, KCN. Il campo Lab consente di selezionare la tipologia di vincolo di simmetria: 1. SYMM, per i vincoli di simmetria; 2. ASYM, per i vincoli di antisimmetria; invece, il campo Normal indica la superficie normale alla direzione scelta. In questo caso, si bloccano tutti gli spostamenti normali alla direzione y. L’ultimo campo KCN specifica il sistema di riferimento. Per applicare il carico: –––––––––– nsel,r,loc,x,b/2 Sf,all,pres,-s0 Nall –––––––––– dove ,con il comando: NSEL, Type, Item, Comp, VMIN, VMAX, VINC, KABS si sono selezionati tutti i nodi che giacciono sulla retta di equazione x=b/2 e con sf si è applicato il carico distribuito. Il comando SF, Nlist, Lab, VALUE, VALUE2, richiede che sino indicati i nodi su cui è applicato il carico (Nlist) e la tipologia di carico (Lab). In riferimento al problema si scrive pres per indicare un carico di pressione. –––––––––– /solu solve –––––––––– /post1 /dscale,1,1 59

plnsol,s,x ––––––––––

Figure 29: Soluzione finale dell’elemento piastra. Come si evince dalla soluzione finale, rappresentata in figura, il punto ’A’ è in trazione mentre il punto ’B’ è in compressione. Il codice sorgente può essere scritto anche in forma parametrica, defininendo i parametri di mesh. In questo modo, è possibile cambiare la dimensione dell’elemento costituente la mesh in modo tale da ottenere una soluzione del problema più accurata. File input completo. –––––––––– finish /clear !DEFINIZIONE DEI PARAMETRI GEOMETRICI b=1 60

w=0.6 r=0.1 !DEFINIZIONE DEI PARAMETRI DEL MATERIALE E=210e+9 ni=0.3 !DEFINIZIONE DEI PARAMENTRI DI CALCOLO s0=100e+6 !DEFINIZIONE PARAMETRI DI MESH inf=1 Mesh1=3*inf Mesh2=2*inf Mesh3=4*inf /prep7 !MODELLAZIONE GEOMETRIA k,1,r $ k,2,2*r $ k,3,b/2 $ k,4„r $ k,5„2*r $ k,6,2*r,2*r k,7,b/2,2*r $ k,8„w/2 $ k,9,2*r,w/2 $ k,10,b/2,w/2 !GENERAZIONE DELLE LINEE l,1,2 $ l,2,3 $ l,4,5 l,2,6 $ l,3,7 $ l,5,6 l,6,7 $ l,5,8 $ l,6,9 l,7,10 $ l,8,9 $ l,9,10 !GENERAZIONE DELL’ARCO k,11 larc,1,4,11,r ldiv,13,0.5 l,12,6 !GENERAZIONE DELLE AREE al,3,6,15,14

61

al,4,1,13,15 al,4,2,5,7 al,8,6,9,11 al,9,7,10,12 !ELEMENTO, MATERIALE et,1,plane42„1 Mp,ex,1,E Mp,nuxy,1,ni !MESH Lsel,r,line„4,6,2 $ lsel,a,line„13,14 $ lsel,a,line„5,11,6 Lesize,all„,Mesh1 Lsel,all Lsel,r,line„1,3,2 $ lsel,a,line„15 Lesize,all„,Mesh1 Lsel,all Lsel,r,line„8,10 Lesize,all„,Mesh2 Lsel,all Lsel,r,line„2,12,5 Lesize,all„,Mesh3 Lsel,all Type,1 $ mat,1 Amesh,all !VINCOLI DI SIMMETRIA nsel,r,loc,x,0 Dsym,symm,x Nall nsel,r,loc,y,0

62

Dsym,symm,y Nall !CARICHI nsel,r,loc,x,b/2 Sf,all,pres,-s0 Nall /solu solve /post1 /dscale,1,1 plnsol,s,x –––––––––– Imponendo differenti valori dei parametri di mesh, attraverso la sintassi: –––––––––– inf=N Mesh1=3*inf Mesh2=2*inf Mesh3=4*inf –––––––––– con N=1,2,3,4, si ottengono diversi valori delle soluzioni, poichè la mesh si infittisce. Per visualizzare i valori delle tensioni: –––––––––– prnsol,s,prin –––––––––– Effettuando la risoluzione del problema per N=1, il coefficiente di concentrazione pari al rapporto del valore ottenuto e la tensione costante σ 0 è: K1 = 3.33.

63

Infittimento della mesh (N=2): –––––––––– inf=2 Mesh1=3*inf Mesh2=2*inf Mesh3=4*inf –––––––––– il valore ottenuto del coefficiente di concentrazione è pari a K2 = 3.55. In maniera analoga per N=3: –––––––––– inf=3 Mesh1=3*inf Mesh2=2*inf Mesh3=4*inf –––––––––– il valore ottenuto è K3 = 3.57. Infine, per N=4: –––––––––– inf=4 Mesh1=3*inf Mesh2=2*inf Mesh3=4*inf –––––––––– si ottiene il valore K4 = 3.57. Attraverso tale analisi di convergenza, è possibile stabilire il livello di accuratezza della mesh necessario per avere risultati affidabili. Lo stesso procedimento può essere ripetuto utilizzando l’elemento Plane182 in condizioni plane stress (keyopt(3)=0). ¯ method (keyopt(1)=0) (integrazione completa sui 4 punti Nella formulazione del B

64

di Gauss e correzione della matrice di compatibilità) –––––––––– Et,1,plane182,0„0 –––––––––– i valori dei coefficienti di concentrazione che si ottengono per gli stessi livelli di mesh prima definiti, sono K1 = 3.33, K2 = 3.55, K3 = 3.57, e K4 = 3.57, che risultanoidentici ai valori ottenuti con l’elemento Plane42 in formulazione base. L’elemento Plane182 da anche la possibilità di utilizzare una formulazione semplificata con integrazione ridotta e controllo dell’hourglassing (keyopt(1)=1). La relativa sintassi è: –––––––––– Et,1,plane182,1„0 –––––––––– I valori che in questo caso si ottengono per il K sono: K1 = 2.50, K2 = 2.93, K3 = 3.10, e K4 = 3.20, leggermente inferiori rispetto a quelli ottenuti precedentemente (plane42 con formulazione con nodi incompatibili). Questo si verifica perchè si riduce l’ordine di integrazione e, quindi, l’accuratezza di calcolo. Nella formulazione Enhanced strain formulation per l’elemento plane182 (keyopt(1)=2): –––––––––– Et,1,plane182,2„0 –––––––––– i valori ottenuti di K sono: K1 = 3.27, K2 = 3.49, K3 = 3.54, e K4 = 3.54. Quest’ultimi sono più vicini alla formulazione base dell’elemento Plane42. La stessa analisi può essere effettuata con l’utilizzo dell’elemento Plane183, il quale è l’elemento piano quadratico a 8 nodi con funzioni di forma quadratiche. L’elemento Plane183 consente di utilizzare sia l’elemento quadrilatero che triangolare impostando rispettivamente un valore del keyopt(1)=0 oppure keyopt(1)=1. Inoltre è dotato della

65

forma degenere cioè l’elemento triangolare quadratico, poichè il software potrebbe non riuscire a utilizzare tutti elementi quadrilateri a causa di alcune irregolarità della geometria. Di seguito, si riportano in tabella i risultati ottenuti attraverso l’utilizzo dei due elementi, Plane42 e Plane182 con le diverse formulazioni. K

42 K(1)=1 42 K(1)=0 182 K(1)=0 182 K(1)=1 182 K(1)=2

K1 (inf = 1) 3.33

3.27

3.33

2.50

3.27

K2 (inf = 2) 3.55

3.49

3.55

2.93

3.49

K3 (inf = 3) 3.57

3.54

3.57

3.10

3.54

K4 (inf = 4) 3.56

3.54

3.57

3.20

3.54

Figure 30: Analisi di convergenza sul valore di K.

66

8

Esercitazione 7: meccanica della frattura linear elastica

Si supponga di dover calcolare il coefficiente di intensificazione degli sforzi KI nella piastra criccata soggetta al carico di trazione σ 0 mostrata in Figura 31

Figure 31: Piastra con cricca. Trattandosi di problema simmetrico, ci si può limitare a modellarne la metà (Figura 32). Ovviamente, in corrispondenza della sezione di simmetria andranno imposti i dovuti vincoli di simmetria (spostamento normale alla superficie pari a zero). Si scelga il sistema di riferimento con origine coincidente con l’apice della cricca. Si può procedere con la definizione dell’area della semipiastra tramite i seguenti comandi noti: –––––––––– finish 67

Figure 32: Sistema simmetrico. /clear /prep7 k,1 $ k,2,-2 $ k,3,14 $ k,4,14,6 $ k,5,-2,6 a,2,1,3,4,5 –––––––––– A questo punto si specificano la tipologia di elementi e le proprietà del materiale: –––––––––– et,1,plane183,1„2 mp,ex,1,206000 $ mp,nuxy,1,0.3 –––––––––– Poiché si tratta di un problema in cui è presente una cricca, andranno utilizzati, in corrispondenza dell’apice della stessa, elementi singolari che si ottengono a partire dagli elementi iso-parametrici quadratici spostando il nodo di mezzeria dei lati convergenti verso l’apice della cricca di una quantità pari a L/4 (dove L è la dimensione del lato). In questo modo gli elementi riproducono correttamente la variazione degli spostamenti e delle tensioni in prossimità della cricca, così come previsto dalle formule di Irwin.

68

Per tale ragione, si fa ricorso all’elemento plane183, che è l’elemento piano quadratico. Il primo keyopt con il valore 1 serve per specificare la forma triangolare dell’elemento. Il terzo keyopt a valore 2 indica che il problema è di tipo plain-strain. Per poter trasformare l’elemento scelto in elemento singolare in ambiente ANSYS si ricorre al comando kscon. –––––––––– kscon,1,1/24,1,12,0.5 –––––––––– kscon richiede prima di tutto di indicare il keypoint corrispondente all’apice della cricca (in questo caso, il numero 1); il secondo campo definisce il raggio della zona attorno all’apice della cricca all’interno della quale ricorrere agli elementi singolari; il valore 1 del terzo campo permette di trasformare l’elemento quadratico in elemento singolare spostando il nodo di mezzeria del lato che converge verso l’apice della cricca ad 1/4 della dimensione del lato stesso dall’apice; il campo 4 serve per indicare il numero di elementi da generare in direzione circonferenziale attorno all’apice della cricca (in questo caso, 12); infine, il campo 5 indica il rapporto tra la dimensione della seconda fila di elementi (ovviamente non singolari) e la prima, in modo da effettuare un "passaggio graduale" dalla zona singolare alla zona non singolare. Si definisce quindi la dimensione di riferimento dell’elemento e si genera la mesh. –––––––––– esize,1 type,1 $ mat,1 $ amesh,1 –––––––––– I comandi che seguono servono ad imporre i vincoli di simmetria in direzione y sulla linea 2, che unisce i keypoints 1 e 3, e il carico di pressione sulla linea 4, che unisce i keypoints 4 e 5. Si noti che la pressione è stata espressa col segno "-" in quanto il programma di default considera positivi i carichi entranti (in questo caso, ciò avrebbe prodotto una compressione). Il comando /psf è utile per visualizzare a

69

video il carico di pressione applicato, specie per un check sul verso dello stesso. –––––––––– lsel,r,line„2 $ nsll,r,1 $ dsym,symm,y $ allsel,all lsel,r,line„4 $ nsll,r,1 $ sf,all,pres,-100 $ allsel,all /psf,pres,norm,2 –––––––––– A questo punto si può entrare in /solution e risolvere, e subito dopo entrare in /post1. –––––––––– /solu solve /post1 –––––––––– Per calcolare il KI è necessario definire un sistema di riferimento locale con origine sull’apice della cricca, assi x ed y parallelo e perpendicolare alla faccia della cricca, rispettivamente. –––––––––– cskp,11,0,1,3,5 $ rsys,11 –––––––––– Con il comando cskp si crea un nuovo sistema di riferimento locale, tramite la definizione di 3 nuovi keypoints. Il primo campo indica il numero che verrà utilizzato da qui in poi per richiamare il nuovo riferimento (11); il secondo campo definisce un numero utile a specificare la tipologia di sistema di riferimento generato: 0 corrisponde ad un sistema di coordinate cartesiane; dopodichè vengono specificati i 3 keypoints utili alla definizione del riferimento, secondo quest’ordine: la prima posizione è riservata al keypoint indicante l’origine del riferimento (1), la seconda serve per il keypoint che individua l’asse x (3) e la terza individua il piano xy (posso indicare uno qualunque dei keypoints tra 4 e 5). Con il comando rsys,11 si sta "attivando"

70

il sistema di riferimento in modo che i calcoli da qui in poi vengano effettuati in tale sistema. Per poter calcolare il KI bisogna definire un path con 3 o 5 punti, a seconda che si stia modellando metà cricca o una cricca intera. In questo caso sono richiesti 3 punti. Il percorso deve passare per un punto sull’apice della cricca e per altri due punti sulla faccia della cricca. Non è importante esattamente a quale distanza debbano trovarsi gli ultimi due punti. Il percorso verrà definito con il comando path: –––––––––– path,ki,3 ppath,1,2 ppath,2,16 ppath,3,10 kcalc,0„1 –––––––––– Il nome del path è "ki", ed è dotato di 3 punti come detto. Con il comando ppath si definiscono tali punti: il primo è il nodo 2, il secondo è il nodo 16 e il terzo è il nodo 10. kcalc permette di calcolare il KI : il primo campo serve per indicare se il problema è plain strain o assialsimmetrico (valore 0) o plain stress (valore 1); la seconda posizione non è generalmente importante; il terzo campo serve per specificare se la cricca è stata modellata a metà o per intero: in questo caso, essendo stata fatta modellazione di metà cricca per un problema simmetrico il valore del campo 3 è 1. Si voglia ora plottare l’andamento della tensione in prossimità della cricca per controllare se effettivamente in quella zona si riscontra una singolarità delle tensioni. –––––––––– pdef,sy,s,y pdef,sy,s,y

71

pdef,sy,s,y plpath,sy,sx,sz –––––––––– Con il comando pdef si sono definite tre variabili, sy, sx, sz, nelle quali vengono memorizzati i valori rispettivamente delle tensioni σ x , σ y , σ z . Con il comando plpath vengono plottate tali variabili lungo il percorso prima definito (Figura 33).

Figure 33: Andamento delle tensioni. File input completo. –––––––––– finish /clear /prep7 k,1 $ k,2,-2 $ k,3,14 $ k,4,14,6 $ k,5,-2,6 a,2,1,3,4,5 et,1,plane183,1„2 mp,ex,1,206000 $ mp,nuxy,1,0.3 72

kscon,1,1/24,1,12,0.5 esize,1 type,1 $ mat,1 $ amesh,1 lsel,r,line„2 $ nsll,r,1 $ dsym,symm,y $ allsel,all lsel,r,line„4 $ nsll,r,1 $ sf,all,pres,-100 $ allsel,all /psf,pres,norm,2 /solu solve /post1 cskp,11,0,1,3,5 $ rsys,11 path,ki,3 ppath,1,2 ppath,2,16 ppath,3,10 kcalc,0„1 pdef,sy,s,y pdef,sy,s,y pdef,sy,s,y plpath,sy,sx,sz ––––––––––

73

9

Esercitazione 8: analisi dinamica strutturale

9.1

Linee guida per la scelta del passo di integrazione (time stepping)

L’accuratezza della soluzione di un problema dinamico dipende dal passo temporale di integrazione. Appare ovvio che la scelta del passo temporale di integrazione è conseguenza di un trade-off tra accuratezza e tempo di calcolo. In questo capitolo, ci si limiterà a parlare del solutore che utilizza il metodo di Newmark. In prima approssimazione, la stima della dimensione del passo temporale può essere ricavata attraverso la seguente relazione:

IT S =

1 20f

(1)

in cui f è la frequenza più di cui si voglia cogliere gli effetti nella risposta del sistema. Come visto in altri capitoli, la risoluzione di un problema avviene seguendo diversi step: 1. applicazione del metodo dell’analisi modale; 2. calcolo delle frequenze naturali del sistema; 3. identificazione della frequenza massima di interesse; Si consideri un sistema con un solo grado di libertà costituito da massa e molla. È possibile plottare il grafico che descrive il periodo di elongazione in funzione del numero di steps per ciclo. Dalla Figura 34 si evince che diminuendo il passo temporale (aumentando il numero di steps temporali per ciclo), i risultati per ogni ciclo per un periodo di elongazione risultano inferiori all’1%. Si prenda in considerazione un sistema sul quale sia applicato un carico che segua una legge temporale ben precisa. In questo caso, il time step deve essere opportu-

74

Figure 34: Periodo di elongazione in funzione del numero di steps temporali per ciclo. namente piccolo da poter seguire in maniera accurata la storia temporale del carico applicato. Nella Figura 35, si nota come la risposta del sistema è in ritardo.

Figure 35: Ritardo della risposta rispetto al carico applicato. Nel caso in cui si verificasse un cambiamento di step di carico conviene utilizzare dei valori più piccoli del passo temporale, il quale è possibile calcolarlo attraverso la seguente relazione:

IT S =

1 . 180f

75

(2)

In questo modo si riduce l’effetto di ritardo della risposta.

9.2

Automatic Time Stepping

La scelta dell’opzione Automatic Time Stepping permette di ottimizzare il passo temporale. In problemi non lineari, inoltre, permette di incrementare i carichi in modo appropriato sulla base delle soluzioni di convergenza ottenute in precedenza.

9.3

Definizione delle condizioni iniziali

Per default in ANSYS Mechanical APDL le condizioni iniziali sono assunte di quiescienza. In generale, per assegnare le condizioni iniziali (spostamento non nullo e/o velocità non nulla), ci si può avvalere del comando IC, la cui sintassi è: IC, NODE, Lab, VALUE, VALUE2, NEND, NINC In ordine occorre indicare il nodo (o i nodi selezionati) su cui imporre la condizione iniziale (NODE), il grado di libertà interessato dalla condizione iniziale (LABEL), il valore iniziale di spostamento (VALUE) e velocità (VALUE2). Il comando IC è consigliabile utilizzarlo quando le condizioni iniziali sono le stesse per l’intera struttura. Quando, invece, si hanno condizioni iniziali diverse in varie parti del sistema è preferibile utilizzare altri metodi per assegnare le condizioni iniziali (onde evitare situazioni di incompatibilità). 9.3.1

Spostamento iniziale zero e velocità iniziale non nulla

La velocità non nulla è applicata imponendo piccoli spostamenti su un piccolo intervallo di tempo sulla parte di struttura dove la velocità deve essere specificata. Ad esempio se u˙ 0 = 0.25, si può applicare uno spostamento di 0.001 su un intervallo di 0.004. ––––––––––

76

TIMINT,OFF nsel,.....

! Selezionare i nodi dove imporre la velocità

d,all,ux,0.001 time,0.004 lswrite,1 ddele,all,ux TIMINT,ON .. . –––––––––– Il comando TIMINT,OFF disattivagli effetti inerziali. 9.3.2

Spostamento e velocità iniziali non nulli

Se, ad esempio, occorre imporre la condizione u0 = 1 e u˙ 0 = 2.5, si può applicare lo spostamento 1 nell’intervallo 0.4 –––––––––– TIMINT,OFF nsel,.....

! Selezionare i nodi dove imporre la velocità

d,all,ux,1.0 time,0.4 lswrite,1 ddele,all,ux TIMINT,ON .. . –––––––––– 9.3.3

Spostamento iniziale non nullo e velocità iniziale nulla

In questo caso, è richiesto l’uso di due substeps (NSUBST,2) con un cambio nello step di applicazione dello spostamento (KBC,1). Senza il campbiamento nel tipo di step (o 77

con solo un substep), lo spostamento imposto varierebbe con il tempo determinando una velocità iniziale non nulla. Quindi, nel caso, ad esempio, di u0 = 1 e u˙ 0 = 0: –––––––––– TIMINT,OFF nsel,.....

! Selezionare i nodi dove imporre la velocità

d,all,ux,1.0 time,.001 nsubst,2 kbc,1

! Stepped loads

lswrite,1 ! Transient solution TIMINT,ON time,...

! Intervallo di tempo realistico

ddelemall,ux kbc,0 .. .

! Ramped loads (se appropriato)

–––––––––– 9.3.4

Accelerazioni iniziale non nulla

E’ possibile applicare un’accelerazione iniziale non nulla, specificando l’accelerazione richiesta (con il comando ACEL) in un piccolo intervallo di tempo. Supponendo di voler applicare un’accelerazione di 9.81 si ha –––––––––– acel„9.81 time,.001 nsubst,2 kbc,1

78

ddele,...

! Perché la struttura deve essere non vincolata nel primo

load step, altrimenti l’accelerazione specificata non avrebbe effetti lswrite,1 .. . ––––––––––

9.4

Esempio di analisi al transitorio

Si vogliano studiare le oscillazioni di un trampolino in seguito al tuffo dell’atleta (Figura 36).

Figure 36: Trampolino. All’estremità del trampolino verrà applicato un carico il cui andamento è riportato nel diagramma di Figura 37. Si procede partendo dalla definizione dei parametri geometrici d’interesse. Si creano prima di tutto le dimensioni del trampolino: –––––––––– /prep7 /TITLE,Transient analysis of a diving board 79

Figure 37: Andamento del carico applicato. !DEFINIZIONE PARAMETRI GEOMETRICI l=4.5 $ w=0.75 $ t=0.2 –––––––––– Per creare la geometria, si definiscono due keypoints: –––––––––– !DEFINIZIONE KEYPOINTS k,1 k,2,0,0,w KGEN,2,1,2,1,l –––––––––– Con il comando kgen così formulato i due keypoint appena definiti vengono copiati a distanza l lungo la direzione x. Il primo 2 indica il numero di volte per il quale l’operazione va ripetuta, ed è seguito dal "pattern" di keypoints da copiare, ovvero i keypoints dal numero 1 al 2 con passo 1; infine, è indicata la distanza a cui copiarli, ovvero in direzione x della quantità l. A questo punto viene generata l’area di base del trampolino che ha per vertici i 4 keypoints appena definiti: –––––––––– 80

!GENERAZIONE AREA A,1,2,4,3 –––––––––– Trattandosi di una struttura in cui una dimensione è trascurabile rispetto alle altre due, con carichi agenti trasversalmente al piano medio, con buona approssimazione la struttura potrà essere trattata con la teoria della piastra, e quindi potranno essere utilizzati elementi shell. –––––––––– !ELEMENT TYPE AND SIZE ET,1,shell63 ESIZE,w/4 /ESHAPE,1,2 !PROPRIETA’ DEL MATERIALE mp,ex,1,30e6$mp,nuxy,1,0.3$mp,dens,1,2000 !SELEZIONE SPESSORE ELEMENTO R,1,t –––––––––– In ANSYS esistono diverse tipolgie di elementi shell. Si ricordi che la matrice di rigidezza dell’elemento shell è ottenuta combinando la matrice di rigidezza dell’elemento piastra (teoria di Reissner-Mindlin) per tener conto degli effetti flessionali con la matrice di rigidezza dell’elemento piano (in condizion plane stress) per tener conto degli effetti membranali. Così come nel caso piano si hanno elementi triangolari e quadrilateri con funzioni di forma lineari o quadratiche, anche per gli elementi shell, in base al tipo di elemento piano scelto, si avranno elementi triangolari, quadrilateri di tipo lineare o quadratico. Si è inoltre visto che per gli elementi piani lineari in ANSYS è possibile utilizzare il plane42 o il plane182, che differiscono per la tipologi di formulazione (con particolare riferimento alla modalità di correzione dello shear-locking). Analogamente, per

81

l’elemento shell è possibile operare una scelta tra l’elemento shell63, che si basa sul plane42 per tener conto degli effetti membranali, o lo shell18 1, in cui i termini di rigidezza relativi ai gradi di libertà membranali sono calcolati sulla base della formulazione del plane182. Con il comando esize si è poi definita la dimensione della mesh. Il comando /eshape serve per una semplice visualizzazione a video dell’elemento: con la sintassi riportata, vengono visualizzati a video gli elementi con il loro spessore. Con il comando che segue viene infatti assegnato lo spessore del trampolino, tramite la definizione di una costante reale. Con il comando amesh si effettua la mesh dell’area 1. –––––––––– ! MESH AREAS AMESH,1 –––––––––– Bisogna ora imporre i vincoli. –––––––––– ! SELEZIONE DEI VINCOLI FISSI SUL BORDO DEL TRAMPOLINO nsel,r,loc,x,0 $ d,all,all $ nsel,all –––––––––– Si sono selezionati tutti i nodi che hanno coordinata x pari a 0 e su tutti i gradi di libertà è stato imposto spostamento nullo. Infine, sono stati riselezionati tutti i nodi. Si inserisce inoltre un ulteriore appoggio, ad una distanza dal bordo fisso pari al 25% della lunghezza del trampolino. –––––––––– !CREAZIONE VINCOLO SUPPORTO AGGIUNTIVO nsel,r,loc,x,0.25*l $ d,all,uy $ nsel,all –––––––––– Si sono selezionati i nodi a distanza 0.25 × l dal bordo vincolato, e si è imposto 82

su questi lo spostamento uy pari a zero. Il modello FEM a questo punto è ultimato. Bisogna applicare i carichi e risolvere il problema. Poiché il problema prevede l’applicazione di un carico variabile nel tempo, si procederà con una analisi al transitorio. La prima grandezza da fissare è l’intervallo temporale ∆t d’integrazione. Si ricordi che lo step temporale deve essere scelto in modo tale da poter permettere di cogliere la natura fisica del fenomeno: non dovrà essere nè troppo grande, in quanto si correrebbe il rischio di non cogliere alcune frequenze di oscillazione, nè troppo piccolo perché si rischierebbe di appesantire il costo computazionale dell’analisi. Si ricordi, inoltre, che la scelta dell’intervallo temporale può avvenire a valle di una analisi a convergenza, ovvero, partendo da un valore dello step non troppo piccolo, effettuare una prima analisi e progressivamente ridurre il ∆t, fermandosi al valore al di sotto del quale si osserva che la soluzione inizia a non variare più apprezzabilmente. In ogni caso, una prima indicazione sull’entità dello step temporale si può ottenere effettuando una preliminare analisi modale. L’analisi modale permette di calcolare le frequenze naturali della struttura e i modi di vibrare ad esse associate. Una volta determinate le frequenze proprie, la scelta del ∆t risulterà facilitata. –––––––––– !ANALISI MODALE PRELIMINARE /solution antype,modal modopt,subsp,5 solve fini –––––––––– Si è quindi entrati in ambiente /solution, richiamando il solutore modal per effettuare l’analisi modale. Il comando modopt serve per specificare l’algoritmo numerico

83

scelto da utilizzare per risolvere il problema autovalori-autovettori: in questo caso è stato scelto l’algoritmo subspace. Infine, è stato specificato il numero di frequenze proprie da estrarre: le prime 5. A questo punto, si esce e si rientra in ambiente /solution (nelle versioni più recenti di ANSYS questo passaggio non è necessario), per immettere i seguenti comandi: –––––––––– /solution expass,on mxpans,5 solve fini –––––––––– Le stringhe di codice appena specificate sono utilizzate per estrarre i modi di vibrare associati alle frequenze naturali individuate poc’anzi, ovvero per calcolare gli autovettori del problema, associati ai 5 autovalori (le frequenze) prima trovati. A questo punto si può entrare in ambiente /postprocessing. –––––––––– /post1 set,list –––––––––– Il comando set,list viene visualizzata una tabella che riporta le 5 frequenze naturali estratte, in ordine crescente. Per poter visualizzare i modi di vibrare, verranno utilizzate le stringhe seguenti: –––––––––– set,1,1 pldisp,1 set,1,2 pldisp,1

84

set,1,3 pldisp,1 set,1,4 pldisp,1 set,1,5 pldisp,1 –––––––––– le quali, nell’ordine, hanno permesso di visualizzare ciascuno dei 5 modi di vibrare. Scegliendo un intervallo temporale pari all’inverso della più piccola frequenza individuata, probabilmente si riusciranno già a cogliere gli aspetti più rilevanti della risposta della struttura; per agire con un grado maggiore di sicurezza, si potrà utilizzare tale valore come punto di partenza per un’analisi a convergenza. A questo punto si può procedere con l’analisi al transitorio. –––––––––– /solution antype,trans trnopt,full !APPLICAZIONE DEL CARICO SECONDO UN APPROCCIO STEPPED kbc,1 –––––––––– antype,trans è utilizzato per richiamare il solutore al transitorio; trnopt,full serve invece per specificare il solutore scelto per l’analisi al transitorio (in questo caso, full, che è basato sull’utilizzo del metodo di Newmark). Per spiegare l’utilità del comando kbc bisogna effettuare delle premesse. La storia temporale imposta permette di individuare fondamentalmente 3 loadsteps, ciascuno corrispondente ai 3 intervalli di tempo in cui è apprezzabile una variazione del carico: - il primo loadstep, da t0 a t1;

85

- il secondo loadstep, da t1 a t2; - il terzo loadstep, da t2 a t3. A questo punto è necessario specificare la modalità con cui il carico deve variare all’interno di ogni loadstep. Col comando kbc è possibile specificare se si desidera applicare il comando come stepped o come ramped. L’approccio ramped è sintetizzato nella Figura 38.

Figure 38: Schematizzazione Ramped. Il carico applicato al’interno del secondo loadstep come ramped sale da 0, valore assunto in t1, a F0 assunto in t2, ovvero come una rampa applicata nell’intero loadstep. Ripetendo l’approccio ramped per il terzo loadstep, si avrebbe una rampa discendente da F0 assunto in t2 a 0 in t3. E’ chiaro, dunque, che l’utilizzo di uno schema ramped richiederebbe una ridefinizione differente dei loadstep. Si ricordi che in generale l’approccio ramped dà meno problemi in un’analisi a convergenza. Altrimenti, per applicare una storia di carico più simile a quella specificata inizialmente, specie se ∆t è sufficientemente piccolo, bisognerà ricorrere ad un approccio stepped. Assegnando il carico come stepped, esso viene portato al valore assegnato nel 86

Figure 39: Schematizzazione Stepped. primo substep del secondo loadstep, e mantenuto costante in tutto il loadstep. Nel terzo loadstep, il carico viene riportato a zero nel primo substep, e mantenuto costante in tutto il loadstep. Dal punto di vista sintattico, per definire un carico ramped occorre emettere il comando kbc,0, per definrlo stepped, kbc,1. Si fissa lo step temporale come –––––––––– f5=4.07 dt=1/f5 –––––––––– Dopo aver definito l’intervallo, esso viene assegnato con il comando deltim. –––––––––– deltim,dt $ autots,on alphad,0.5 $ betad,0.002 outres„all $ outpr,all –––––––––– 87

autots,on serve per fare in modo che il software monitori automaticamente il valore dello step temporale scelto, correggendone l’entità sia che risulti troppo piccolo che troppo grande. alphad e betad servono per definire i coefficienti utili alla determinazione della matrice di dissipazione viscosa secondo il metodo di Rayleigh. outres„all e outpr,all servono per specificare che in output e nel file dei risultati vanno salvati i calcoli. A questo punto si procede con la creazione del carico variabile nei loadsteps. –––––––––– nsel,r,loc,x,l f,all,fy,0 $ nsel,all time,0.01 lswrite,1 –––––––––– Si sono selezionati i nodi all’estremità del trampolino ed è stato loro assegnato il carico fy pari a 0. Si è assegnato il tempo t1=0.01 s coincidente con l’inizializzazione della forza. Il loadstep è definito tramite il comando lswrite,1. –––––––––– nsel,r,loc,x,l fdel,all,all $ f,all,fy,-100 $ nsel,all time,1.5 lswrite,2 –––––––––– Si sono selezionati i nodi all’estremità del trampolino ed è stato loro assegnato il carico fy pari a -100 N, concorde con il riferimento. Si è poi assegnato il tempo t2=1.5 s e definito il secondo loadstep con il comando lswrite,2. –––––––––– nsel,r,loc,x,l fdel,all,all $ f,all,fy,0 $ nsel,all

88

time,12 lswrite,3 lssolve,1,3 fini –––––––––– Si sono selezionati i nodi all’estremità del trampolino ed è stato loro assegnato il carico fy pari a 0. Si è assegnato il tempo t3=12 s e definito il terzo loadstep (lswrite,3). Con il comando lssolve,1,3 vengono risolti i loadstep dal numero 1 al numero 3. A questo punto si può entrare in post-processing per la visualizzazione dei risultati, per esempio la variazione dello spostamento uy di un nodo specifico, il nodo 32 (Figura 40). –––––––––– /post26 nsol,2,32,u,y plvar,2 ––––––––––

Figure 40: Andamento dello spostamento lungo y. 89

Il comando nsol permette di creare una variabile che può inglobare una grandezza a scelta: in questo caso, assegnando il nome 2 alla variabile, e specificando il nodo 32, rimangono due campi per la specificazione della variabile, ovvero lo spostamento u lungo y. Col comando plvar,2 viene plottata la variabile 2. –––––––––– File input completo. /prep7 /TITLE,Transient analysis of a diving board l=4.5 $ w=0.75 $ t=0.2 k,1 k,2,0,0,w KGEN,2,1,2,1,l A,1,2,4,3 ET,1,shell63 ESIZE,w/4 /ESHAPE,1,2 mp,ex,1,30e6 $ mp,nuxy,1,0.3 $ mp,dens,1,2000 R,1,t AMESH,1 nsel,r,loc,x,0 $ d,all,all $ nsel,all nsel,r,loc,x,0.25*l $ d,all,uy $ nsel,all /solution antype,modal modopt,subsp,5 solve fini /solution expass,on

90

mxpans,5 solve fini /post1 set,list set,1,1 pldisp,1 set,1,2 pldisp,1 set,1,3 pldisp,1 set,1,4 pldisp,1 set,1,5 pldisp,1 /solution antype,trans trnopt,full kbc,1 f5=4.07 dt=1/f5 deltim,dt $ autots,on alphad,0.5 $ betad,0.002 outres„all $ outpr,all nsel,r,loc,x,l f,all,fy,0 $ nsel,all time,0.01 lswrite,1

91

nsel,r,loc,x,l fdel,all,all $ f,all,fy,-100 $ nsel,all time,1.5 lswrite,2 nsel,r,loc,x,l fdel,all,all $ f,all,fy,0 $ nsel,all time,12 lswrite,3 lssolve,1,3 fini /post26 nsol,2,32,u,y plvar,2 ––––––––––

92

10

Esercitazione 9: submodeling

Si consideri il problema in Figura mostrato in Figura 41.

Figure 41: Schema di carico e dimensioni di un albero rotante. Dati: dp =30 mm da =40 mm bc =10 mm Lp =30 mm La =30 mm r=2 mm F=4000N La tecnica del submodelling consente di eseguire analisi di dettaglio (con mesh fitte) senza la necessità di avere grosse risorse di calcolo. Infatti, è possibile fare analisi accurate, ad esempio, dello stato di tensione in una parte di un componente effettuando dapprima un calcolo su un modello ‘grossolano’ con mesh non fitta e, successivamente, effettuare un’analisi solo sulla parte del componente di interesse, il ‘sottomodello’, utilizzando una mesh fitta e imponendo come carichi gli spostamenti interpolati dall’analisi precedente sulle sezioni di ‘taglio’ del sottomodell. Il problema è affrontato in 4 steps: 1. Realizzazione della geometria, priva di raccordo (modello grossolano) e con mesh 93

rada; 2. Realizzazione di un sottomodello contenente il raccordo, e con mesh fitta; 3. Interpolazione degli spostamenti sulle superfici di contorno del sottomodello; 4. Soluzione del sottomodello. Con questa procedura, il "modello grossolano" è utilizzzato per determinare gli spostamenti in corrispondenza delle sezioni di ‘estrazione’ del "modello dettagliato". Di seguito gli steps per la risoluzione. • Fase-1: modello coarse (Figura 42)

Figure 42: Modello coarse.

• Fase-2: sottomodello (Figura 43).

• Fase-3: interpolazione spostamenti (Figura 44). 94

Figure 43: Sottomodello.

Figure 44: Interpolazione spostamenti.

95

–––––––––– finish /clear resume,coarse,db /post1 –––––––––– Per effettuare l’interpolazionde sui nodi salvati nel file: –––––––––– cbdof,sottomodelo,node„sottomodello,cbdo finish –––––––––– • Fase-4: soluzione sottomodello. –––––––––– finish /clear resume,sottomodello,db /solu –––––––––– per applicare gli spostamenti interpolati al sottomodello: –––––––––– /input,sottomodello,cbdo finish –––––––––– File input completo. –––––––––– finish /clear 96

/filename,coarse dp=30 $ da=40 $ bc=10 $ lp=30 $ la=200 $ f=4000 /prep7 et,1,plane42 et,2,solid45 mp,ex,1,206000 mp,nuxy,1,0.3 k,1,0,(da-dp)/2 k,2,lp,(da-dp)/2 k,3,lp,0 k,4,la/2,0 k,5,la/2,da/2 k,6,lp,da/2 k,7,0,da/2 a,1,2,6,7 a,2,3,4,5,6 esize,2.5 amesh,1,2 esize„10 vrotat,all„„„7,6,180 /solu ! VINCOLI SIMMETERIA nsel,r,loc,z,0 $ dsym,symm,z $ nall nsel,r,loc,x,la/2 $ dsym,symm,x $ nall ! VINCOLO APPOGGIO nsel,r,loc,x,bc $ nsel,r,loc,z,0 $ nsel,r,loc,y,da/2 $ d,all,uy $ nall ! FORZA APPLICATA nsel,r,loc,x,la/2 $ nsel,r,loc,z,0 $ nsel,r,loc,y,da $ f,all,fy,-f/4 $ nall

97

save solve finish /clear /filename,submodel dp=30 $ da=40 $ bc=10 $ lp=30 $ la=200 $ r=2 /prep7 et,1,plane42 et,2,solid45 mp,ex,1,206000 mp,nuxy,1,0.3 k,1,2*lp/3,(da-dp)/2 k,2,lp-r,(da-dp)/2 k,3,lp,(da-dp)/2-r k,4,lp,0 k,5,lp+10,0 k,6,lp+10,da/2 k,7,2*lp/3,da/2 l,1,2 $ larc,2,3,1,r a,1,2,3,4,5,6,7 esize,1 amesh,1 esize„10 vrotat,all„„„7,6,180 ! I nodi delle superfici di ‘taglio’ del sottomodello, devono essere selezionati e scritti su un file nsel,r,loc,x,2*lp/3 $ nsel,a,loc,x,lp+10 nwrite $ nall ! crea il file submodel.node

98

save finish /clear resume,coarse,db /post1 file,coarse,rst set,1,1 cbdof,submodel,node„submodel,cbdo ! effettua l’interpolazione per i nodi salvati nel filecbdof finish /clear resume,submodel,db /solu antype,static,new /input,submodel,cbdo ! applica gli spostamenti interpolati al sottomodello ! VINCOLI SIMMETERIA nsel,r,loc,z,0 $ dsym,symm,z $ nall solve finish /post1 /dscale,1,1 plnsol,s,eqv ––––––––––

99

11

Esercitazione 10: analisi termo-strutturale

Si considera una lamina in acciaio vincolata a due strutture solide ad una temperatura di riferimento di 0◦ C (273 K) (Figura 45). Una delle strutture solide viene riscaldata ad una temperatura di 75◦ C (348 K). Quando il calore viene trasferito dalla struttura solida alla lamina, questa tenterà di espandersi. Tuttavia, questo non può verificarsi poichè la lamina è vincolata tra le due strutture, quindi si creeranno delle tensioni interne alla stessa. Si chiede di determinare la soluzione allo steady-state.

Figure 45: Lamina vincolata a due blocchi a differenti temperature. Non c’è alcun carico applicato alla lamina, ma solo una variazione di temperatura di 75 gradi Celsius. La lamina è in acciaio con un modulo di elasticità di 200 GPa, conduttività termica di 60, 5 [ mWK ] e coefficiente di dilatazione termica pari a 12 · 10−6 [ K1 ]. File input completo. –––––––––– finish /clear /prep7 100

k,1,0,0 k,2,1,0 l,1,2 et,1,link33 r,1,4e-4, mp,kxx,1,60.5 esize,0.1 lmesh,all ! QUI VIENE CREATO UN AMBIENTE FISICO CHIAMATO THERMAL IN CUI SONO IMMAGAZZINATE LE INFORMAZIONI RELATIVE AL PROBLEMA TERMICO physics,write,thermal physics,clear etchg,tts

! TALE COMADO SERVE PER

PASSARE DALL’ELEMENTO TERMICO link33 A QUELLO STRUTTURALE link180 mp,ex,1,200e9 mp,prxy,1,0.3 mp,alpx,1,12e-6 ! QUI VIENE CREATO UN AMBIENTE FISICO CHIAMATO STRUCT IN CUI SONO IMMAGAZZINATE LE INFORMAZIONI RELATIVE AL PROBLEMA STRUTTURALE physics,write,struct physics,clear finish /solu antype,0

101

physics,read,thermal

! VIENE RICHIAMATO IL PROB-

LEMA TERMICO dk,1,temp,348 solve finish /solu physics,read,struct

! VIENE RICHIAMATO IL PROB-

LEMA STRUTTURALE ldread,temp„„„rth

! VENGONO RICHIAMATE LE

TEMPERATURA OTTENUTE DALLA SOLUZIONE DEL PROBLEMA TERMICO CHE DIVENTANO IL CARICO IN INPUT PER IL PROBLEMA STRUTTURALE tref,273

! VIENE ASSEGNATA LA TEM-

PERATURA DI RIFERIMENTO DELLA LAMINA RISPETTO ALLA QUALE CALCOLARE LA VARIAZIONE DI TEMPERATURA dk,1,all,0 dk,2,UX,0 solve finish /post1 ! CREAZIONE DI UNA TABELLA CON GLI SFORZI DI COMPRESSIONE E SUA VISUALIZZAZIONE A VIDEO etable,CompStress,LS,1 PRETAB,CompStress

102

––––––––––

103

12

Esercitazione 11: Contatto hertziano

Di seguito viene mostrato come risolvere il problema di contatto tra il cilindro e il semipiano schematizzati in Figura 46.

Figure 46: Esempio di contatto Hertziano. File input completo. –––––––––– FINISH /CLEAR /UNITS,SI /FILENAME,CONTACT /TITLE,HERTZ CONTACT /PREP7 ! DEFINIZIONE PARAMETRI GEOMETRIA 104

R=0.002 ! RAGGIO CILINDRO LY=0.01 ! ALTEZZA SEMIPIANO LX=0.01 ! LARGHEZZA SEMIPIANO ! DEFINIZIONE KEYPOINTS K,1,0,0 $ K,2,LX,0 $ K,3,LX,-LY $ K,4,0,-LY $ K,5„R $ K,6,R,2*R $ K,7„2*R ! DEFINIZIONE LINEE L,1,2 $ L,2,3 $ L,3,4 $ L,4,1 $ LARC,5,6,7,R $ L,6,7 $ L,7,5 ! DEFINIZIONE AREE AL,1,2,3,4 AL,5,6,7 ! MATERIALE E TIPO ELEMENTO ET,1,PLANE82„,2 MP,EX,1,206E9 MP,NUXY,1,0.3 ! CONTROLLI DI MESH LESIZE,1„,40,20 $ LESIZE,2„,40,20 $ LESIZE,3„,40,0.05 $ LESIZE,4„,40,0.05 LESIZE,5„,80 $ LESIZE,6„,80 $ LESIZE,7„,80 ! MESH TYPE,1 $ MAT,1 $ AMESH,ALL ! ELEMENTI DI CONTATTO ET,3,TARGE169 ET,4,CONTA172 KEYOPT,4,1,0 ! GRADI DI LIBERTÀ KEYOPT,4,7,2 ! TIME INCREMENTATION (reasonable increment for next substep) ! GENERATE THE TARGET SURFACE LSEL,S„,1 CM,_TARGET,LINE

105

TYPE,3 NSLL,S,1 $ ESLN,S,0 ESURF,ALL $ ALLSEL,ALL ! GENERATE THE CONTACT SURFACE LSEL,S„,5 CM,_CONTACT,LINE TYPE,4 NSLL,S,1 $ ESLN,S,0 ESURF,ALL $ ALLSEL,ALL ! VINCOLI ! SPOSTAMENTI NULLI ALLA BASE LSEL,R,LINE„3 $ NSLL,R,1 $ D,ALL,UY $ ALLSEL,ALL KSEL,R,KP„3 $ NSLK,R $ D,ALL,UX $ ALLSEL,ALL ! VINCOLI DI SIMMETRIA LSEL,R,LINE„4,7,3 $ NSLL,R,1 $ DSYM,SYMM,X,0 $ ALLSEL,ALL ! SOLUZIONE /SOLU AUTOTS,ON KBC,0 ! PRIMO STEP DI CARICO LSEL,R,LINE„6 $ NSLL,R,1 $ D,ALL,UY,-R $ ALLSEL,ALL LSWRITE,1 ! SECONDO STEP DI CARICO *DO,I,1,2 LSEL,R,LINE„6 $ NSLL,R,1 $ DCUM,ADD $ D,ALL,UY,-R/100 $ ALLSEL,ALL LSWRITE,I+1 *ENDDO LSSOLVE,1,3

106

/POST1 /DSCALE,1,1 PATH,PLOT,17 PPATH,1,4962 *DO,I,2,17 PPATH,I,4963+2*(I-1) *ENDDO PDEF,P,CONT,PRES PDEF,PENE,CONT,PENE PLPATH,P $ PCALC,INTG,CARICO,P,S,1 ––––––––––

107

13

Esercitazione 12: Eigenvalue Buckling Analysis

Il fenomeno del buckling può essere affrontato attraverso 2 tipi di analisi: 1. Lineare 2. Non lineare, i cui risultati sono più attendibili. Si consideri la Figura 47:

Figure 47: Fenomeno del buckling rappresentanto nel diagramma forza-spostamento. Con riferimento alla Figura 47a ad un incremento della forza corrisponde un aumento dello spostamento. Ad un determinato valore di carico la curva raggiunge un massimo e successivamente il carico si riduce con lo spostamento indicando che la rigidezza della struttura è cambiata e, in particolare, si è ridotta notevolmente. Oltre un determinato valore di spostamento, la struttura recupera le caratterisitche di rigidezza. Si osservi che l’andamento della curva inizialmente è lineare, successivamente avvicinandosi al punto di biforcazione (o punto di instabilità) ci si discosta dalla linearità. Il metodo non lineare consiste nell’attivare il solutore non-lineare e incrementare il carico applicandolo a piccoli step, in maniera tale da ricostruire la curva forzaspostamento. Quando il solutore non sarà più in grado di ottenere la convergenza della soluzione si sarà giunti al punto di biforcazione della curva e quindi si sarà individuata la condizione di carico che porta ad instabilità. 108

Tuttavia, si potrebbe verificare la non convergenza della soluzione per altre ragioni anche in campo lineare. In quest’ultimo caso, le ragioni di tale fenomeno devono essere cercate altrove poichè è improbabile che l’instabilità elastica si verifichi in campo lineare. In generale, il punto di instabilità è preceduo da un tratto di comportamento non lineare. L’analisi lineare (schematizzata in Figura 47b) è basata sull’ipotesi che il sistema si comporta sempre linearmente sino al punto di biforcazione. Ciò, porta a sovrastimare leggermente il carico critico d’instabilità. Il metodo lineare è basato sulla definizione di un problema agli autovalori-autovettori la cui soluzione fornisce le indicazioni sui carichi di instabilità. In generale, prima di effettuare un’analisi non lineare è sempre buona norma condurre una preliminare analisi lineare per individuare il range in cui ci si aspetta l’occorrenza del fenomeno di instabilità.

13.1

Analisi lineare del fenomeno del buckling

Per affrontare il problema del buckling attraverso l’analisi lineare si deve definire il modello FEM, e risolvere il problema staticamente. In seguit, si effettua l’analisi di buckling, dalla quale si possono estrapolare i diversi modi di instabilità con i rispettivi valori di carico. Affinchè il problema sia considerato lineare, sia le equazioni costitutive che di congruenza devono essere lineari. Inoltre, è interessante notare come gli autovalori rappresentino dei fattori di scala dei carichi. Dalla Figura 48 si evince che applicando un carico A = 1.0 il primo autovalore calcolato è pari a 100. Questo significa che il carico critico d’instabilità risulta 100 volte il carico applicato: Fcr = 100 (A + W0 ), dove W0 è il peso proprio della struttura. Modificando il carico in A = 100 si trova un autovalore λ = 1.1, da cui Fcr = 1.1 (A + W0 ) e cosi via fino a determinare λ = 1.

109

Figure 48: Autovalori e carico applicato. 13.1.1

Carico critico euleriano

E’ noto che il carico critico Euleriano per una trave incastrata ad un estremo e lbera all’altro è Fcr,E =

π 2 EJmin 4 l2

(3)

con ovvio significato dei simboli. Si voglia determinare il carico critico per la trave mostrata in Figura 49.

Figure 49: Struttura sottoposta a un carico di punta F. Si effettua un’analisi di buckling lineare. Si parte con la definizione dei parametrici 110

geometrici, del materiale e di carico: –––––––––– L=1 r=0.05 E=210e+9 ni=0.3 F=1 –––––––––– Successivamente si entra in ambiente pre-processing, in modo tale da definire i keypoints, le linee e la tipologia di elemento da utilizzare per il modello FEM. –––––––––– k,1 k,2„L l,1,2 et,1,beam188„,3 sectype,1,beam,csolid secdata,r mp,ex,1,E mp,nuxy,1,ni –––––––––– Si procede con la realizzazione della mesh dell’elemento. –––––––––– esize„10 lmesh,1 –––––––––– Dopodichè, si inseriscono rispettivamente i vincoli e i carichi. –––––––––– d,1,all

111

f,2,fy,-F –––––––––– Prima di effettuare l’analisi si devono attivare gli effetti "prestress" attraverso la sintassi: PSTRES, Key. Nel campo KEY si inserisce: • OFF, per non includere gli effetti di "prestress" (default); • ON, per includere gli effetti di "prestress". Quindi, –––––––––– Pstres,on –––––––––– Per calcolare i diversi modi di buckling si deve effettuare una nuova analisi richiamando l’algoritmo BUCKLE attraverso il comando ANTYPE. –––––––––– antype,buckle bucopt,lanb,1 mxpand,1„,yes outress,all,all outpr,nsol,all –––––––––– Con il comando BUCOPT è possibile specificare quale algoritmo utilizzare per risolvere il problema agli autovalori e agli autovettori. La sintassi del comando BUCOPT è: BUCOPT, Method, NMODE, SHIFT, LDMULTE, RangeKey che richiede in ordine la tipologia di algortimo e il numero dei modi di buckling da calcolare. 112

Per ottenere la deformata associata ai modi di buckling calcolati precedentemente, si deve procedere con la loro espansione. Ciò è possibile attraverso il comando MXPAND. La sintassi di tale comando è: MXPAND, NMODE, FREQB, FREQE, Elcalc, SIGNIF, MSUPkey, ModeSelMethod che richiede nel campo NMODE il numero dei modi da espandere, e nel campo Elcalc se si vogliono calcolare anche gli sforzi associati a quei modi di buckling. Di seguito si riporta il codice sorgente completo: –––––––––– finish /clear !PARAMETRI GEOMETRICI L=1 r=0.05 E=210e+9 ni=0.3 F=1 !PARAMETRI DEL MATERIALE, SEZIONE, ELEMENTO /prep7 k,1 k,2„L l,1,2 et,1,beam188„,3 sectype,1,beam,csolid secdata,r mp,ex,1,E mp,nuxy,1,ni !PARAMETRI DI MESH

113

esize„10 lmesh,1 !PARAMETRI DI CARICO E VINCOLI d,1,all f,2,fy,-F /solu Pstres,on Solve finish !ANALISI DEL BUCKLING /solu antype,buckle bucopt,lanb,1 mxpand,1„,yes outress,all,all outpr,nsol,all solve /post1 /dscale,1,1 set,1,1 pldisp,1 –––––––––– La Figura 50 riporta la prima forma modale di buckling (instabilità flessionale nel piano xz ). In Figura 51 è, invece riportato l’autovalore corrispondente. Avendo applicato un carico unitario, esso coincide proprio con il carico critico d’instabilità. Applicando la formula (3) si ottiene Fcr,E =

π 3 Er4 = 2.543 × 106 16 L2 114

(4)

Figure 50: Deformata riguardante il primo modo di instabilità flessionale. essendo Jmin = πr4 /4. Tale valore è molto vicino a quello calcolato in ANSYS. Si tenga conto che la formula (3) è ricavata con riferimento alla teoria di Eulero-Bernoulli della trave, mentre in ANSYS è stato utilizzato l’elemento beam188 basato sulla teoria di Timoshenko.

Figure 51: Valori estrapolati. Risultati identici si ottengono per il modo di instabilità flessionale nel piano yz. Per estrarre quest’ultimo, è necessario modificare il codice sorgente specificando di dover calcolare i primi due modi di instabilità: –––––––––– 115

antype,buckle bucopt,lanb,2 mxpand,2„,yes outress,all,all outpr,nsol,all –––––––––– Se si è interessati anche al terzo modo di buckling, occorre specificare di estrarre i primi 3 modi: –––––––––– antype,buckle bucopt,lanb,3 mxpand,3„,yes outress,all,all outpr,nsol,all –––––––––– La deformata corrispondente al terzo modo di instabilità flessionale è riportata in Figura 52 Gli autovalori sono mostrati in Figura 53.

14

Estrazione della matrice di rigidezza

Quando si abbia necessità di conoscere i coefficienti della matrice di rigidezza, è possibile estrarn il valore facendo eseguire al software la seguente soluzione: –––––––––– /solu antype,substr seopt,file,1,1 m,all,all

116

Figure 52: Deformata riguardante il secondo modo di instabilità flessionale.

Figure 53: Valori estrapolati.

117

/output,filename

!filename è il nome del file dove verrà memo-

rizzata la matrice solve /output

118