======Programma agli elementi finiti triangolari======
=====Dimensionamento=====
La prima parte del programma si occupa della __dichiarazione__ delle variabili necessarie per definire la struttura;
PARAMETER (NODES=9)
PARAMETER (NELEMS=8)
PARAMETER (IFMAX=3,ICNMAX=3)
DIMENSION XY(2,NODES),NVERT(3,NELEMS),IFXY(IFMAX),FXY(2,IFMAX),
ICNXY(2,ICNMAX),CNXY(2,ICNMAX),STRUTK(2*NODES,2*NODES),
FORCE(2*NODES),UV(2*NODES)
**NODES**: numero di nodi\\
**NELEMS**: numero degli elementi\\
**IFMAX**: numero di nodi caricati\\
**ICNMAX**: numero di nodi vincolati\\
La stringa "DIMENSION" permette al programma di allocare le matrici necessarie per i calcoli, dedicandogli un appropriato spazio in memoria. In Fortran 77 (linguaggio di questo codice) è importante conoscere la dimensione di ogni array di tipo matrice o vettore perchè l'allocazione di memoria per quest'ultimo avvenga correttamente; ad esempio NODES deve essere una costante poichè è usata per definire il numero di colonne della matrice XY, agendo diversamente non conoscerei la dimensione di XY e quindi quanta memoria dedicargli.\\
**XY**: matrice formata da due righe e NODES colonne, contiene nella prima riga l'ascissa di ogni nodo e nella seconda l'ordinata; si rammenta, servendosi del caso specifico, che in Fortran la numerazione parte, salvo diversa segnalazione, sempre da 1 e arriva al NODES-esimo elemento mentre se si partisse da 0 si arriverebbe al NODES-esimo -1 elemento\\
{{ :wikipaom2016:img_1462.jpg?300 |}}
**NVERT**: matrice formata da tre righe e NELEMS colonne, contiene in ogni colonna il numero (label o etichetta) di ogni nodo che compone l'elemento; è detta matrice di connettività dell'elemento\\
{{ :wikipaom2016:img_1465.jpg?300 |}}
Con i dati indicati in precedenza è possibile disegnare sul piano la struttura per capire meglio di cosa stiamo parlando.
{{ :wikipaom2016:img_1463.jpg?300 |}}
**IFXY**: vettore formato da IFMAX componenti, contiene i nomi (label) dei nodi caricati\\
**FXY**: matrice formata da due righe e IFMAX colonne, contiene nella prima riga le componenti lungo x delle forze applicate al nodo e nella seconda le componenti lungo y\\
L'impossibilità di salvare all'interno dello stesso array un reale ed un intero obbliga a creare IFXY e FXY per distinguere tra l'array che indica quale nodo è caricato e l'array che indica come è caricato quello stesso nodo (si è ragionato nello stesso modo per ICNXY e CNXY).\\
{{ :wikipaom2016:img_1464.jpg?300 |}}
**ICNXY**: matrice formata da due righe e ICNMAX colonne, contiene nella prima riga i nomi (label) dei nodi vincolati e nella seconda il tipo di vincolamento secondo il seguente schema:\\
* spostamento in direzione x associato al numero 1
* spostamento in direzione y associato al numero 2
* spostamento in direzione x e y associato al numero 3
{{ :wikipaom2016:img_1466.jpg?300 |}}
{{ :wikipaom2016:img_1467.jpg?300 |}}
**CNXY**: matrice formata da due righe e ICNMAX colonne, contiene nella prima riga il valore dello spostamento imposto lungo x e nella seconda lungo y\\
{{ :wikipaom2016:img_1468.jpg?300 |}}
Notare che la scrittura "0." indica un numero reale.\\ Ora la struttura può essere disegnata inserendo i vincolamenti ed i caricamenti ai nodi interessari.
{{ :wikipaom2016:img_1469.jpg?300 |}}
**STRUTK**: matrice di rigidezza dell'intera struttura\\
**FORCE**: vettore delle forze imposte\\
**UV**: vettore spostamenti nodali\\
DIMENSION D(3,3),B(3,6),BT(6,3),DB(3,6),ELK(6,6),IPOINT(6),
DELTAEL(6),SIGMAEL(3)
Dopo il dimensionamento degli elementi variabili si termina l'allocazione di matrici e vettori defininendo le grandezze dipendenti dalla tipologia di risoluzione scelta (tria 3 nodi).\\
**D**: matrice di legame costitutivo (legge di Hooke), lavorando in 3D diventa una matrice 6x6\\
**B**: matrice che lega spostamenti e deformazione\\
**BT**: trasposta di B\\
**DB**: prodotto tra matrice D e B (matrice temporanea)\\
**ELK**: matrice di rigidezza degli elementi (aggiornata per ogni elemento, evitando così di avere un array in tre dimensioni 6x6xNELEMS)\\
**IPOINT**: mappatura dei gdl globali e locali per ogni elemento. Ad esempio per l'elemento 3 si ha:\\
{{ :wikipaom2016:img_1470.jpg?300 |}}
**DELTAEL**: vettore di 6 spostamenti locali dei 3 vertici dell'elemento\\
**SIGMAEL**: vettore tre elementi che contiene la σx, la σy e la τxy
=====MAIN=====
OPEN(1,FILE='mesh_tria_comp.dat',STATUS='OLD')
OPEN(2,FILE='mesh_tria_comp.3.txt',STATUS='UNKNOWN')
CALL READIN(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,FXY,
ICNMAX,ICNXY,CNXY)
CALL ECHO(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,FXY,ICNMAX,
ICNXY,CNXY)
CALL CLEAR(2*NODES,2*NODES,STRUTK)
CALL CLEAR(2*NODES,1,FORCE)
C CALL CLEAR(3,6,B)
C CALL CLEAR(6,6,ELK)
CALL DMAT(YM,PR,ITDP,D)
**mesh_tria_comp.dat**: file __esistente__ (STATUS=OLD) in cui si trovano tutte le informazioni note del problema (coordinate nodali, Modulo di Young, ecc...)\\
**mesh_tria_comp.3.txt**: file __che potrebbe esistere__ (STATUS=UNKNOWN) in cui riporto tutte le informazioni che stampo a video durante l'esecuzione del programma (il file viene sovrascritto se già esistente)\\
Se il comando OPEN è dato nel modo descritto sopra, il codice cerca il file all'interno della cartella corrente; in alternativa si può specificare il percorso del file qualora si trovi in una directory remota.\\
====Subroutine Readin=====
Sfruttando più cicli DO per estendere la lettura ad ogni nodo, questa subroutine va a leggere dalla prima riga del file precedentemente aperto e definito come 1 (mesh_tria_comp.dat) rispettivamente il modulo di Young (YM), il modulo di Poisson (PR) e l'indice discriminante, 0 o 1 rispettivamente, tra lo stato di tensione o deformazione piana (ITDP). Legge inoltre dalle successive righe le coordinate nodali, le coordinate dei vertici degli elementi triangolari, le forze nodali e le condizioni di vincolamento.\\
SUBROUTINE READIN(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,
FXY,ICNMAX,ICNXY,CNXY)
DIMENSION XY(2,NODES),NVERT(3,NELEMS),IFXY(IFMAX),FXY(2,IFMAX),
ICNXY(2,ICNMAX),CNXY(2,ICNMAX)
C Lettura di YM, PR e ITDP:
READ(1,*)YM,PR,ITDP
C Lettura delle coordinate nodali (x,y):
DO 10,I10=1,NODES
READ(1,*)NODE,XY(1,I10),XY(2,I10)
IF(I10.NE.NODE)THEN
write(*,*) 'errore nella numerazione dei nodi sul file di input'
stop
END IF
10 CONTINUE
C Lettura dei vertici degli elementi triangolari:
DO 20,I20=1,NELEMS
READ(1,*)NELEM,NVERT(1,I20),NVERT(2,I20),NVERT(3,I20)
IF(I20.NE.NELEM)THEN
write(*,*) 'errore nella numerazione degli elementi'
stop
END IF
20 CONTINUE
C Lettura delle forze nodali:
DO 30,I30=1,IFMAX
READ(1,*)IFXY(I30),FXY(1,I30),FXY(2,I30)
30 CONTINUE
C Lettura delle condizioni di vincolamento:
DO 40,I40=1,ICNMAX
READ(1,*)ICNXY(1,I40),ICNXY(2,I40),CNXY(1,I40),CNXY(2,I40)
40 CONTINUE
RETURN
END
Per ogni ciclo DO presente nella subroutine il programma esegue un controllo di sequenzialità dei parametri (coordianate nodali, spostamenti e forze applicate) dei file di input.\\ Si dice che durante questa operazione di lettura le variabili sono "passate per riferimento" ovvero che il valore di YM, ad esempio, non viene salvato per intero ma viene salvata solo la posizione della prima cella dell'indirizzo che contiene, scritto in byte, YM.\\
====Subroutine Echo====
Subroutine che scrive nel file 2 (mesh_tria_comp.3.txt) tutto quello che è stato eseguito nella subroutine READIN, si nota infatti che la struttura delle due è identica con l'accortezza di modificare il comando READ(1,*) con WRITE(*,2).\\
SUBROUTINE ECHO(YM,PR,ITDP,NODES,XY,NELEMS,NVERT,IFMAX,IFXY,FXY,
ICNMAX,ICNXY,CNXY)
DIMENSION XY(2,NODES),NVERT(3,NELEMS),IFXY(IFMAX),FXY(2,IFMAX),
ICNXY(2,ICNMAX),CNXY(2,ICNMAX)
C Controllo di YM e PR:
WRITE(2,*)'MODULO DI YOUNG: ',YM
WRITE(2,*)'COEFFICIENTE DI POISSON: ',PR
C Controllo di ITDP:
IF(ITDP.EQ.0)THEN
WRITE(2,*)'STATO DI TENSIONE PIANA'
END IF
IF(ITDP.EQ.1)THEN
WRITE(2,*)'STATO DI DEFORMAZIONE PIANA'
END IF
C Controllo delle coordinate nodali:
WRITE(2,*)'NODO: X: Y: '
DO 10,I10=1,NODES
WRITE(2,*)I10,XY(1,I10),XY(2,I10)
10 CONTINUE
C Controllo dei vertici degli elementi trangolari:
WRITE(2,*)'ELEMENTO: Vi: Vj: Vk: '
DO 20,I20=1,NELEMS
WRITE(2,*)I20,NVERT(1,I20),NVERT(2,I20),NVERT(3,I20)
20 CONTINUE
C Controllo delle forze:
WRITE(2,*)'NODO: FX: FY: '
DO 30,I30=1,IFMAX
WRITE(2,*)IFXY(I30),FXY(1,I30),FXY(2,I30)
30 CONTINUE
C Controllo delle condizioni di vincolamento:
WRITE(2,*)'NODO: TIPO VINCOLAMENTO: DELTAX: DELTAY: '
DO 40,I40=1,ICNMAX
WRITE(2,*)ICNXY(1,I40),ICNXY(2,I40),CNXY(1,I40),CNXY(2,I40)
40 CONTINUE
RETURN
END
====Subroutine Clear====
Inizializzazione a zero degli elementi di una matrice servendosi di due cicli DO nidificati (uno all'interno dell'altro).
SUBROUTINE CLEAR(K1,K2,A)
DIMENSION A(K1,K2)
DO 10,I10=1,K1
DO 20,I20=1,K2
A(I10,I20)=0.
20 CONTINUE
10 CONTINUE
RETURN
END
====Subroutine Dmat====
Costruzione della matrice D di applicazione del legame Hookeano tra tensioni e deformazioni: **σ** = D * **ε**. la matrice D è definita come un array 3x3 a partire da YM, PR, ITDP.
SUBROUTINE DMAT(YM,PR,ITDP,D)
DIMENSION D(3,3)
IF(ITDP.EQ.0)GO TO 100
IF(ITDP.EQ.1)GO TO 200
C Per Stato di Tensione Piana:
100 CONTINUE
CTP=YM/(1-PR**2)
G=YM/(2*(1+PR))
C Struttura della matrice D per Stato di Tensione Piana:
C | 1*CTP PR*CTP 0 |
C D = | PR*CTP 1*CTP 0 |
C | 0 0 G |
D(1,1)=1.
D(1,2)=PR
D(2,1)=PR
D(2,2)=1.
DO 10,I10=1,2
DO 20,I20=1,2
D(I10,I20)=D(I10,I20)*CTP
20 CONTINUE
10 CONTINUE
D(1,3)=0.
D(2,3)=0.
D(3,1)=0.
D(3,2)=0.
D(3,3)=G
GOTO 300
C Per Stato di Deformazione Piana:
200 CONTINUE
CDP=YM*(1-PR)/((1+PR)*(1-2*PR))
CPR=PR/(1-PR)
G=YM/((1+PR)*2)
C Struttura della matrice D per Stato di Deformazione Piana:
C | 1*CDP CDP*CPR 0 |
C D = | CDP*CPR 1*CDP 0 |
C | 0 0 G |
D(1,1)=1.
D(1,2)=CPR
D(2,1)=CPR
D(2,2)=1.
DO 30,I30=1,2
DO 40,I40=1,2
D(I30,I40)=D(I30,I40)*CDP
40 CONTINUE
30 CONTINUE
D(1,3)=0.
D(2,3)=0.
D(3,1)=0.
D(3,2)=0.
D(3,3)=G
GOTO 300
300 CONTINUE
RETURN
END
Per decidere il tipo di legame costitutivo da applicare per la generazione della matrice D, la subroutine legge il valore di ITDP attraverso due costrutti IF: questo genere di approccio potrebbe portare ad errore nel caso in cui ITDP fosse uguale a qualunque altro numero che non sia zero o uno, in quanto la subroutine continuerebbe con l'esecuzione delle istruzioni nell'ordine in cui sono fornite; assocerebbe quindi ad uno stato tensionale sconosciuto (ITDP diverso da 0 e 1) lo stato di tensione piana (poiché è il primo che legge).
=====MAIN=====
Al termine della subroutine precedente ho definito la matrice D e torno al programma principale e attivo un DO che va da 1 a NELEMS (numero elementi). Queste istruzioni vengono ripetute una volta per elemento e I10 è lo specifico numero dell'elemento su cui opero o per il quale opero.
DO 10,I10=1,NELEMS
XI=XY(1,NVERT(1,I10))
YI=XY(2,NVERT(1,I10)) !CALL AREA, KELMAT
XJ=XY(1,NVERT(2,I10)) !POINTVT eASSEMBL
YJ=XY(2,NVERT(2,I10)) !sono dentro lo stesso
XK=XY(1,NVERT(3,I10)) !ciclo DO 10, cos
YK=XY(2,NVERT(3,I10)) !vengono calcolati per
CALL AREA(XI,YI,XJ,YJ,XK,YK,I10,TWODEL) !ogni elemento.
CALL KELMAT(XI,YI,XJ,YJ,XK,YK,TWODEL,I10,D,ELK)
C *****Assemblaggio***********
CALL POINTVT(NELEMS,I10,NVERT,IPOINT)
CALL ASSEMBL(ELK,NODES,IPOINT,STRUTK)
10 CONTINUE
Per ogni elemento prendo il nome dei nodi che gli appartengono dalla matrice NVERT. Scelto il nodo, prendo le coordinate x e y all'interno della matrice XY (scegliendo tra la riga 1 o la riga 2).
====Subroutine Area====
Calcolo della superficie DEL di ciascun elemento finito triangolare a partire dalle coordinate nodali dei vertici. Controllo di dimensioni e segno dell'area calcolata. Le coordinate nodali vengono prese, come di consueto, in verso antiorario. Con TWODEL si identifica il determinante della matrice, ovvero il doppio della superficie.
SUBROUTINE AREA(XI,YI,XJ,YJ,XK,YK,NELEM,TWODEL)
TWODEL=XJ*YK+XK*YI+XI*YJ-(XK*YJ+XI*YK+XJ*YI)
AREAMIN=1.E-6
C Controllo sulle dimensioni dell'area calcolata:
IF(ABS((TWODEL/2)).LE.AREAMIN)THEN
WRITE(*,*)'ATTENZIONE! AREA DELL*ELEMENTO TROPPO PICCOLA'
WRITE(*,*)NELEM
WRITE(*,*)'AREA=',TWODEL/2
END IF
C Controllo sul segno dell'area calcolata:
IF((TWODEL/2).LT.-AREAMIN)THEN
WRITE(*,*)'ATTENZIONE! AREA DELL*ELEMENTO NEGATIVA'
WRITE(*,*)NELEM
WRITE(*,*)'AREA=',TWODEL/2
END IF
RETURN
END
====Subroutine Kelmat====
Costruzione della matrice di rigidezza ELK (6x6, reale) di ciascun elemento finito triangolare attraverso la formula: F = ELK * DELTA => ELK = A * (BT * D * B). Le coordinate nodali vengono prese, come di consueto, in verso antiorario. Gli input sono le coordinate nodali dei nodi, il determinante, l'indice di elemento (utile solo per un eventuale warning) e la matrice di legame costitutivo (D comune a tutti gli elementi).
SUBROUTINE KELMAT(XI,YI,XJ,YJ,XK,YK,TWODEL,NELEM,D,ELK)
DIMENSION B(3,6),BT(6,3),D(3,3),DB(3,6),ELK(6,6)
CALL CLEAR(3,6,B)
CALL BMAT(XI,YI,XJ,YJ,XK,YK,TWODEL,B)
CALL PRODMAT(3,3,6,D,B,DB)
CALL TRSPMAT(3,6,B,BT)
CALL PRODMAT(6,3,6,BT,DB,ELK)
DEL=TWODEL/2
DO 10,I10=1,6
DO 20,I20=1,6
ELK(I10,I20)=DEL*ELK(I10,I20)
20 CONTINUE
10 CONTINUE
RETURN
END
La matrice D, passata dal main, viene allocata in memoria al lancio della subroutine, mentre la matrice B (locale) richiede un allocamento all'inizio della subroutine, con relativa chiamata della subroutine CLEAR per inizializzazione della matrice stessa.
====Subroutine Pointvt====
Costruzione del vettore puntatore IPOINT necessario ad effettuare l'assemblaggio della matrice ELK di rigidezza di ciascun elemento finito triangolare (NELEM) nella matrice STRUTK di rigidezza globale della struttura.\\
Nelle celle di IPOINT viene scritto quale è il relativo grado di libertà in numerazione globale, viene riportato l'esempio relativo all'elemento 3.
SUBROUTINE POINTVT(NELEMS,NELEM,NVERT,IPOINT)
DIMENSION NVERT(3,NELEMS),IPOINT(6)
IPOINT(1)=NVERT(1,NELEM)*2-1
IPOINT(2)=NVERT(1,NELEM)*2
IPOINT(3)=NVERT(2,NELEM)*2-1
IPOINT(4)=NVERT(2,NELEM)*2
IPOINT(5)=NVERT(3,NELEM)*2-1
IPOINT(6)=NVERT(3,NELEM)*2
RETURN
END
Come input prende il numero di elementi (per sapere la lunghezza di NVERT), lo specifico elemento (per sapere la colonna di NVERT da consultare), NVERT e IPOINT come risultato.\\
Viene trovato poi il grado di libertà locale associato che riempie il vettore IPOINT che viene poi passato alla chiamante con il comando RETURN. La chiamante lo passa poi pari pari alla subroutine ASSEMBL.
====Subroutine Assembl====
Assemblaggio della matrice ELK di rigidezza di ciascun elemento finito triangolare nella matrice STRUTK di rigidezza globale della struttura.\\
La matrice STRUTK viene inzializzata a zero nel MAIN PROGRAM, pertanto si ritiene superflua la sua inizializzazione in questa subroutine.
SUBROUTINE ASSEMBL(ELK,NODES,IPOINT,STRUTK)
DIMENSION ELK(6,6),IPOINT(6),STRUTK(2*NODES,2*NODES)
DO 10,I10=1,6
DO 20,I20=1,6
IND10=IPOINT(I10)
IND20=IPOINT(I20)
STRUTK(IND10,IND20)=STRUTK(IND10,IND20)+ELK(I10,I20)
20 CONTINUE
10 CONTINUE
RETURN
END
Nota: il vettore NODES viene utilizzato unicamente per il dimensionamento della matrice STRUTK.
La riga fondamentale della sub è quella che associa ad un valore della matrice di rigidezza globale della struttura il suo vecchio valore + il valore da aggiungere dovuto al nodo considerato, dato che al valore del termine STRUTK(I10,I20) della matrice di rigidezza globale contribuiscono in genere più elementi finiti.
=====MAIN=====
C Per il Caricamento si richiama questa sub
CALL FORCES(IFMAX,IFXY,FXY,NODES,FORCE)
====Subroutine Forces====
Riempimento del vettore globale delle forze nodali FORCE.
SUBROUTINE FORCES(IFMAX,IFXY,FXY,NODES,FORCE)
DIMENSION IFXY(IFMAX),FXY(2,IFMAX),FORCE(2*NODES)
DO 10,I10=1,IFMAX
INDX=IFXY(I10)*2-1
INDY=IFXY(I10)*2
FORCE(INDX)=FXY(1,I10)
FORCE(INDY)=FXY(2,I10)
10 CONTINUE
RETURN
END
Dove la sub importa dal main il numero di nodi caricati (IFMAX), il vettore contenente il nome dei nodi caricati (IFXY) e la matrice contenente i valori delle forze lungo x e y per ogni singolo nodo. Il ciclo DO legge i nodi caricati e salva i valori delle forze nel vettore dei termini noti.
=====MAIN=====
====Subroutine Bmat====
Costruzione della matrice B di collegamento tra gli spostamenti nodali e le deformazioni di ciascun elemento finito triangolare attraverso la formula: EPSILON = B * DELTA.\\
Le coordinate nodali vengono prese, come di consueto, in verso antiorario.\\
Infine, poiché i termini nulli, sempre nelle stesse posizioni, vengono gi
inizializzati a zero inizialmente, se ne riportano le rispettive assegnazioni in commento.
SUBROUTINE BMAT(XI,YI,XJ,YJ,XK,YK,TWODEL,B)
DIMENSION B(3,6)
AI=XJ*YK-XK*YJ !
AJ=XK*YI-XI*YK ! Questi valori servono per il calcolo dell'area.
AK=XI*YJ-XJ*YI !
BI=YJ-YK
BJ=YK-YI
BK=YI-YJ
CI=XK-XJ
CJ=XI-XK
CK=XJ-XI
B(1,1)=BI
B(1,2)=0.
B(1,3)=BJ
B(1,4)=0.
B(1,5)=BK
B(1,6)=0.
B(2,1)=0.
B(2,2)=CI
B(2,3)=0.
B(2,4)=CJ
B(2,5)=0.
B(2,6)=CK
B(3,1)=CI
B(3,2)=BI
B(3,3)=CJ
B(3,4)=BJ
B(3,5)=CK
B(3,6)=BK
TWODEL=AI+AJ+AK
DO 10,I10=1,3
DO 20,I20=1,6
B(I10,I20)=B(I10,I20)/TWODEL
20 CONTINUE
10 CONTINUE
RETURN
END
====Subroutine Prodmat====
Esecuzione del prodotto AB di due matrici A e B; prende in input il numero di righe della prima matrice, il numero di colonne della prima matrice (uguale alle righe della seconda) e numero di colonne della seconda matrice.
SUBROUTINE PRODMAT(K1,K2,K3,A,B,AB)
DIMENSION A(K1,K2),B(K2,K3),AB(K1,K3)
DO 10,I10=1,K1
DO 30,I30=1,K3
AB(I10,I30)=0. !Inizializzazione a zero della matrice prodotto.
DO 20,I20=1,K2
AB(I10,I30)=AB(I10,I30)+A(I10,I20)*B(I20,I30) !Prodotto.
20 CONTINUE
30 CONTINUE
10 CONTINUE
RETURN
END
Il prodotto è una operazione computazionalmente molto costosa: date due matrici quadrate di ordine "n", la complessità computazionale del calcolo del loro prodotto risulta di n^3.\\
I due cicli più esterni (I10 e I30) scorrono sugli elementi della matrice prodotto e, per ogni elemento della matrice prodotto si fa una sommatoria fra la I10-esima riga della matrice del primo fattore e la colonna I30 del secondo fattore di cui fai il prodotto scalare.\\ La sommatoria ha bisogno di una variabile di accumulo che appunto è la cella della matrice di destinazione, AB (I10, I30) che deve essere all'inizio azzerata. Per ogni I20 viene ridefinita la variabile di accumulo con un nuovo contributo.
{{:wikipaom2016:main_tria_base3_v002.for|versione a fine lezione}}
~~DISCUSSION~~