Teoria 4

Da CdM_unimore.

Determinazione delle forze nodali

Nella lezione precedente abbiamo visto come determinare le forze sui nodi di un elemento triangolare sfruttando il teorema dei lavori virtuali. Si può giungere allo stesso risultato , seguendo una diversa via, applicando il teorema di Castigliano. Per prima cosa, si parte dalla definizione di energia interna

Per il teorema di Castigliano si può scrivere

A questo punto, si sfrutta la definizione della derivata di uno scalare rispetto a un vettore

dove x è un vettore e K una matrice. Applicando la formula al caso in esame si ottiene

Si ottiene esattamente lo stesso risultato che si è ottenuto con l'approccio classico.

Per avere una conferma della formula utilizzata, è possibile fare un esempio considerando il caso piano, con matrice e vettori di dimensione 2. Chiaramente, non si tratta di una vera e propria dimostrazione, ma è comunque un metodo valido per capire il significato del procedimento utilizzato.
Siano x un vettore colonna 2x1 e K una matrice 2x2

A questo punto, procedo con le derivate parziali:

Si osservi che la matrice K di rigidezza dell'elemento triangolare è simmetrica, per cui risulta che . Le formule precedenti diventano dunque

Scritto in forma vettoriale:

Subroutine PRODMAT

Andiamo a definire il prodotto di matrici tramite un algoritmo Fortran:

SUBROUTINE PRODMAT (K1,K2,K3,A,B,AB)

ove ,, e sono numeri interi definiti grazie alla lettera K ( le lettere I, J, K, L, M ed N sottintendono numeri interi,quindi ci permettono di omettere integer,comando che ricordiamo definisce variabili intere in for77), A e B sono le matrici che utilizzeremo come variabili di ingresso. Infine AB è la matrice che verrà fornita in output.
Predefiniamo le dimensioni delle matrici A, B, AB con in comando dimension:

DIMENSION A(K1,K2), B(K2,K3), AB(K1,K3)

Il nostro scopo è riempire la matrice AB, quindi partiamo subito dalla variabile di uscita:

DO 10,I10=1,K1
DO 30,I30=1,K3
AB(I10,I30)=0

In questo modo si inizializza la matrice AB, come una matrice di soli zeri. Alla matrice di output piena di zeri andrà poi sommata la matrice risultante dal prodotto della matrice A per la matrice B, ottenendo così la matrice output di prodotto, si ha bisogno quindi di un valore iniziale definito che crescerà ciclicamente.

DO 20,I20=1,K2
AB(I10,I30)=AB(I10,I30)+A(I10,I20)*B(I20,I30)

Il comando permette la moltiplicazione riga per colonna della matrice.

20 CONTINUE
30 CONTINUE
10 CONTINUE
RETURN
END

NB: è importante non intrecciare i cicli DO; si iniziano in sequenza i DO 10, 30 e 20 e si chiudono nell'ordine 20,30 e 10. Il comando RETURN serve per tornare al MAIN.
In questo modo si è creata una subroutine che può essere richiamata nel main ogniqualvolta ci sia bisogno di effettuare un prodotto tra matrici. Sarà sufficiente scrivere: call nome_subroutine(var_ingresso, var_uscita). In tal caso si avrebbe:

call PRODMAT(K1, K2, K3, A, B, AB)

Chiaramente, tutte le varibili in ingresso devono essere definite e tutte quelle di uscita devono essere determinate all'interno della subroutine, altrimenti si registrerà un errore.

Subroutine KELMAT

Studiamo ora come costruire la matrice di rigidezza K, data dal prodotto fra uno scalare A, la matrice B trasposta (6x3), la matrice D(3x3) e la matrice B(3x6):

SUBROUTINE KELMAT(XI,YI,XJ,YJ,XK,YK,YM,PR,3,3,6)

dove KELMAT sta per K Elements Matrix, YM sta per Young Modulus, PR sta per Poisson Ratio.

DIMENSION B(3,6),BT(6,3),D(3,3),ELK(6,6)

Abbiamo sostituito ELK a KEL per evitare che il programma interpreti la scrittura come numero intero a causa della lettera K iniziale.
Ora procediamo richiamando le subroutine definite nelle scorse lezioni per ottenere il risultato finale:

CALL BMAT(XI,YI,XJ,YJ,XK,YK,B)

che calcola la matrice B che collega gli spostamenti nodali alle deformazioni dell'elemento, distinguendo inoltre se l'area è positiva o negativa.

CALL TRSP(3,6,B,BT)

La subroutine richiamata sopra traspone la matrice.

CALL DMAT(YM,PR,D)

infine, con il comando sopra si calcolano gli elementi della matrice D; in alternativa potrei inserire la matrice D nell'elenco delle variabili per evitare che venga ricalcolata ogni volta.

CALL PRODMAT (3,3,6,D,B,DB)
CALL PRODMAT (6,3,6,BT,DB,ELK)

Il prodotto di D per B comporta le variabili 3,3,6 poichè la prima matrice è una 3x3 e la seconda una 3x6. Il numero di colonne della prima e il numero di righe della seconda coincidono, quindi possono essere assimilati ad una sola variabile. Da qui la presenza di 3 sole variabili anzichè 4. Stesso discorso per il prodotto di BT(6X3) per DB(3X6).

AREA2 =( XJ*YK + XK*YI + XI*YJ - XK*YJ - XI*YK - XJ*YI )
AREA=AREA2/2
DO10,I10=1,6
DO20,I20=1,6
ELK(I10,I20)=ELK(I10,I20)*AREA
20 CONTINUE
10 CONTINUE
RETURN
END

Esempio

Vediamo ora un esempio applicativo del programma Fortran al metodo degli elementi finiti. Lo scopo è quello di indicizzare vettori appartenenti ai vari elementi finiti oggetto di calcolo. Nella prima fase abbiamo il procedimento detto 'Assemblaggio' nel quale si associano le informazioni di ambito locale alla matrice di struttura globale. In figura è rappresentato un sistema cartesiano di sei elementi finiti triangolari costituiti da 8 nodi. I numeri cerchiati rappresentano la numerazione degli elementi, i numeri in nero rappresentano i nodi.
Immagine sistema cartesiano.png


Come dati abbiamo le coordinate nodali

1=0.,0.              5=0.,1.
2=1.,0.              6=1.,1.
3=2.,0.              7=2.,1.
4=3.,0.              8=3.,1.

e gli elementi con i rispettivi nodi interessati

1=>1,2,5
2=>5,2,6
3=>2,3,6
4=>6,3,7
5=>3,4,7
6=>7,4,8

NB: i nodi degli elementi vanno sempre contati in senso antiorario poichè per convenzione si considera che un'area sia positiva se calcolata nel suddetto senso. Questa convenzione è valida solo nel caso piano e non è trasferibile nello spazio (non esiste la concezione di senso orario e antiorario nello spazio).
Come primo esempio prendiamo l'elemento 1.
Scriviamo (in parte) la matrice rigidezza dell'elemento:

nel sistema del singolo elemento chiameremo il nodo 1 i, il nodo 2 j e il nodo 5 k. Ogni nodo ha bisogno di due righe della matrice per essere descritto, quindi i occuperà la prima e la seconda riga, j la terza e la quarta, k la quinta e la sesta.
Passando dalla matrice rigidezza dell'elemento alla matrice di rigidezza della struttura globale, bisogna rinominare adeguatamente i nodi, utilizzando la formula n matrice elemento=(2n-1, 2n) matrice struttura. Spiegando meglio tramite un esempio, l'indice k associato al nodo 5 occuperà le righe 5*2-1 e 5*2 (9 e 10) della matrice dell'intera struttura.
Detto ciò, se si prende ad esempio l'elemento K(3,5) della matrice di rigidezza dell'elemento, esso corrisponderà a K(3,9) nella matrice di rigidezza della struttura globale.




                   1     i     K1,1   K1,2   K1,3   K1,4   K1,5   K1,6
2 K2,1
3 j K3,1 K3,9
4 K4,1
9 k K5,1
10 K6,1

Immagine matrice elemento 1.png

Immagine elemento 1.png

Di seguito viene riportato un ulteriore esempio, che ripete il procedimento visto sopra per l'elemento 5, dove i,j,k rappresentano rispettivamente i nodi 3,4,7.


                   1   i  K1,1   K1,2   K1,3   K1,4   K1,5   K1,6
2 K2,1
3 j K3,1 K7,13
4 K4,1
9 k K5,1
10 K6,1


Immagine elemento 2.png

ora scriviamo la subroutine che permette lo svolgimento di questo procedimento:

SUBROUTINE POINTVR (NELEMS,NELEM)

dove POINTVR sta per point vector, NELEMS sta per numero di elementi e NELEM è l'elemento in questione.

DIMENSION NVERT (3,NELEMS), IPOINT(6)

IPOINT è una notazione implicita affinchè il programma veda il numero come intero, esso memorizza nelle sue sei caselle corrispondenti ai sei indici locali i sei indici globali. Mentre il vettore NVERT serve per determinare a quale nodo, nella numerazione globale, arrivano i tre vertici dell'elemento finito triangolare.

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
STOP
END



Pressione di frontiera di un tubo elastoplastico (esercizio pagina 19)

PROGRAMMA TUBO ELASTOPLASTICO

PRINT*, 'battibrs,ri,reb'
READ*, rs,ri,re
PRINT*, 'rhoditentativob'
READ*, rho
Pf = rs/2[2*log(rho/ri)+1-(rho/re)**2]
PRINT*,'Pf=b',Pf


VARIANTE DELL'ESERCIZIO PROPOSTO → PROBLEMA INVERSO

PRINT*,'battibrs,ri,re,Pfottimale'
READ*,rs,ri,re,Pfottimale
10 PRINT*, 'rhoditentativob'
READ*, rho
Pf = rs/2[2*log(rho/ri)+1-(rho/ri)**2]
PRINT*,'Pf=b', Pf
IF (Pf.LT.Pfottimale) THEN
PRINT*,'devi aumentare rho'
ELSE
PRINT*, 'devi diminuire rho'
ENDIF
GOTO10
STOP
END



Si tratta di un algoritmo iterativo dove ad ogni iterazione si inserisce il valore di rho(variabile raggio di frontiera) di tentativo e viene calcolata la nuova pressione di frontiera fino alla convergenza tra quest'ultima e la pressione ottimale(fissata dal progettista). Ovviamente il costrutto IF è stato impostato tenendo conto(come si nota dalla formula di Pf sopra)che la pressione di frontiera Pf è monotonicamente crescente con il raggio ρ .

N.B.
Nella stesura del programma sono stati inseriti alcuni spazi e alcuni caratteri maiuscoli ai fini di una più semplice lettura degli appunti.