Teoria 7

Da CdM_unimore.

Rho tubo elastoplastico con metodo di Newton delle tangenti


Si considera l'applicazione del metodo di Newton delle tangenti al caso di tubo elastoplastico. In generale, la terminologia prevede che si parli di "metodo di Newton" (Newton.pdf) se si va a trattare una singola equazione, "metodo di Newton-Rapson" se invece si va a trattare un sistema di equazioni.


  • Il metodo di Newton è impiegato per approssimare la soluzione di una equazione nella forma , pertanto quando una generica curva interseca l'asse delle ascisse. Il metodo viene applicato dopo aver determinato un intervallo che contiene una sola radice. Il metodo prevede di considerare un punto della curva che abbia un ascissa compresa nell'intervallo individuato (può essere anche un estremo) da cui far partire l'algoritmo. Si considera la tangente alla curva nel punto della curva di ascissa considerata, e se questa interseca l'asse delle ascisse in un punto compreso tra gli estremi a e b dell'intervallo, tale punto sarà una soluzione di prima approssimazione. Procedendo con il metodo, si traccia la tangente alla curva nel punto di ascissa pari alla soluzione appena trovata (si determina dunque una tangente diversa dalla precedente) e si considera come nuova soluzione il punto di intersezione con l'asse x di questa nuova tangente alla curva. Questo processo viene iterato più volte, cercando di raggiungere la convergenza del calcolo. Esiste un criterio di convergenza globale del metodo di Newton che dimostra, sotto precise ipotesi, la convergenza del metodo stesso.


Si considera una funzione , dove , variabile indipendente, è il raggio di frontiera plastica nel tubo, cioè il raggio, compreso tra il raggio interno e il raggio esterno, che delimita la zona del tubo (compresa quindi tra raggio interno e raggio di frontiera) plasticizzata sotto l'azione di una pressione interna sufficientemente elevata.


Grafico metodo di newton.png

L'equazione della retta tangente alla curva in è:

dove è il nostro starting point.

Infatti indicando con b l'ipotenusa del triangolo rettangolo

possiamo scrivere





quindi



In seconda iterazione si calcolerà



Il processo di iterazione finirà quando verrà trovato un valore contenuto in un intervallo piccolo a piacere da noi predefinito.

La funzione del tubo elastoplastico è:


vogliamo all'interno del programma annullarla in , per far questo nel codice useremo:


Le condizioni per le quali il metodo di Newton converge globalmente sono:

Sia ossia funzione continua nell'intervallo [a,b], con derivata prima e derivata seconda continua nel medesimo intervallo

1) sempre soddisfatta perché Ri è sempre minore di Re

2) sempre soddisfatta perché la funzione è derivabile in tutto il dominio

3) sempre soddisfatta perché la funzione è concava in tutto il dominio

4)


Tale condizione non è soddisfatta per tutti i valori di primo tentativo, ad esempio se scegliessimo come valore di primo tentativo ,la tangente non intersecherebbe mai l’asse all'interno del dominio, in quanto . Il punto d'intersezione della retta tangente con l'asse delle ascisse deve cadere all'interno dell'intervallo in cui è definita la funzione; se poniamo , il punto d'intersezione coincide con il punto di massimo relativo, per cui la derivata prima è nulla e la tangente è orizzontale. Bisogna quindi ricordare che se la tangente nel punto è nulla allora la funzione risulta essere divergente e la condizione 4) non è verificata. Per noi l'estremo coincide con e con .


Programma Fortran



Il programma fortran per il calcolo di con il metodo di newton è:


Prog metodo di newton.PNG


  CODICE
     Program Calcolo_rho_con_Newton
     
     Print*,'inserisci Rs, Ri, Re, PFopt,Nø iterazioni max'
     Read*, Rs, Ri, Re, Pfopt, n
     Epsi=1.E-4
     
     C poniamo un rho di primo tentativo = Ri per essere sicuro che il metodo converga. Se scegliessimo invece la media tra 
     C i raggi c’è il rischio di non cadere all’interno del dominio della funzione. Se addirittura prendessimo Re il metodo 
     C non convergerebbe perché, come è stato precedentemente detto, la condizione 4) non è rispettata.
     
     Rho=Ri
     
     C ora procediamo con la diagnostica interna:
     
     if (Pfopt.GE.Rs*log(Re/Ri)) then
         print*,'WARNING:BOOM'
         stop
         elseif (Pfopt .LE. (Rs/2*(1-(Ri/Re)**2))) then
         print*,'Il tubo non si deforma plasticamente'
         stop
     endif
     do i=1,n
       f=(Rs/2)*(2*log(rho/Ri)+1-(rho/Re)**2)-Pfopt
       if (abs(f) .LT. epsi) then
       go to 5
       else
       df= ((2/rho)-((2*rho)/Re**2)*Rs)/2
       rho= rho-(f/df)
       endif
     enddo
     print*,'WARNING! n_max raggiunto'
     stop
  5  print*,'raggio (rho)= '
     Stop
     End

Portale incastrato con coppia agente sulla traversa


Valutiamo l'andamento del Momento Flettente considerando le seguenti regole:

  • conservazione dell'angolo retto
  • conservazione della curvatura

Struttura iniziale.png Struttura con carichi.png

Dalla figura appena rappresentata, si vede chiaramente che il problema è antisimmetrico.
L'andamento del Momento Flettente corretto sarà da ricercarsi in una delle due soluzioni che verranno ora presentate:

Casi mf.png

Il 1° Caso però non è corretto perchè non è antisimmetrico. Inoltre nella figura sottostante, che rappresenta la deformata della trave nel 1° Caso, a cavallo dell'angolo retto di sinistra, i due flessi sono diversi e quindi questa non potrà essere una soluzione valida.

Caso 1 flesso.png Caso 3 flessi.png

In arancio l'andamento qualitativo del momento flettente del portale:

Mf qualitativo.png


Programma agli elementi finiti (pag. 133 PA)



Una traccia del MAIN:problema principale.

Il MAIN è rapido perché abbiamo costruito molte routine in precedenza e ci basterà richiamarle.

Per la lettura dei dati uso una subroutine chiamata readin, → CALL READIN a cui farò leggere i valori di
E, v,
xy (2,nodes),
NVERT(3,nelems)(NVERT dice il numero del nodi che sono i 3 vertici di ogni elemento. Indica la topologia dell’elemento),
Fxy(2, IFMAX)(Fxy sono le forze nodali lungo x e lungo y, IFMAX è ,invece, il numero dei nodi caricati),
IFxy(IFMAX)(puntatore, se il primo nodo a essere caricato è il 27, il primo numero dentro IFMAX è il 27),
[dopo i nodi caricati ho bisogno dei nodi vincolati dandogli gli spostamenti imposti i quali spesso sono 0 in quanto sono incastri],
CNxy (2,ICNMAX)(ICNMAX sono il numero massimo dei nodi vincolati. L'indice 2 fa si che se al primo posto compare 1 lo spostamento (vincolo) è imposto lungo x, se 2 lo spostamento(vincolo) è imposto lungo y, se 3 lungo x e y),
ICNxy(2,ICNMAX)(2:la prima riga indica il nodo considerato, la seconda riga indica il tipo di vincolo)

Si ricorda la coincidenza fra CNXY e CNST come anche fra ICNXY e ICNST (notazioni usate nella scrittura delle subroutine durante la lezione di Teoria_6 )

 CALL READIN(YM, PR, ITDP, NODES, XY, NELEMS, NVERT, IFMAX
1IFXY, FXY, ICNMAX, ICNXY, CNXY)
CALL ECHO (YM, PR, NODES, XY, NELEMS, NVERT, IFMAX 1IFXY, FXY, ICNMAX, ICNXY, CNXY)

Ora dobbiamo inizializzare a 0 tutto quello che successivamente verrà riempito ma solo parzialmente.
DMAT ha un bel po’ di zeri perché il legame tra e è solo diagonale.
La matrice B andrà inizializzata a zero, in quanto come è possibile vedere a pagina 118 della dispensa del Professor Strozzi, B è riempita di zero per quasi la metà.
Il vettore dei termini noti FORCE andrà inizializzato.
Allora richiamiamo la subroutine clear che inizializzava a zero vettori e matrici

CALL CLEAR(2*NODES, 2*NODES, TOTK)
CALL CLEAR(2*NODES, 1, FORCE)
CALL CLEAR(3, 6, B)
CALL DMAT(YM, PR, ITDP,D)

Ora bisogna creare la matrice di rigidezza della struttura, la si costruisce con un "DO" che spazzola gli elementi. La matrice di rigidezza va a raggruppare le matrici di rigidezza di ogni elemento.

DO 10, I10= 1, NELEMS
XI= XY(1,NVERT(1,I10))
YI= XY(2,NVERT(1,I10))
XJ= XY(1,NVERT(2,I10))
YJ= XY(2,NVERT(2,I10))
XK= XY(1,NVERT(3,I10))
YK= XY(2,NVERT(3,I10))

XY scorre sui nodi. NVERT centra con gli elementi, ha dentro solo dei numeri che sono i 3 nodi che arrivano ai 3 vertici di ogni elemento. I10 è l’indice corrente e scorre sugli elementi. Il DO fondamentale è sull’elemento così si trova la matrice di rigidezza dell’elemento.
Ricordando che la matrice di rigidezza dell'elemento è:




e siccome D è uguale per tutti gli elementi, allora non la metto dentro il DO.
richiamo la subroutine che mi da la matrice di rigidezza dell'elemento:

CALL ELKMAT


Ora dobbiamo fare l’assemblaggio, quindi richiamiamo le subroutine del puntatore e dellì'assemblaggio:

CALL POINTER(Nelems,Nelem,Nvert, IPOINT)
CALL ASSEMBL(Nodes, ELK, IPOINT, TOTK)
10 CONTINUE

Ho costruito, quindi, la K della struttura. (Se la struttura ha 50 nodi la K è una matrice 100*100)

Ora dobbiamo impostare la relazione


.


Questo è un problema agli spostamenti quindi le forze sono note e gli spostamenti sono incogniti. Fino a questo punto tutti gli spostamenti sono incogniti quindi la struttura è labile e quindi un programma serio non girerebbe.

CALL FORCES(Nodes, IFMAX, IFXY, FXY, FORCE)
CALL CNSTNG(Nodes,ICNMAX,ICNST,CNST,FORCE,TOTK)

messi a posto i vincoli ora posso risolvere il problema

CALL GAUSS

Gauss parte dall’equazione ricavando (gli spostamenti nodali). Da li in poi c’è un post processo durante il quale non è più un incognita ma è nota.
Se mi interessano le tensioni su un solo elemento è considerato noto e avrò bisogno di un altro "DO" che scorre sugli elementi e, elemento per elemento, attraverso la relazione :




Calcola le 3 sigma dentro l’elemento.
Il Mark cerca di tirarti fuori un output più gradevole possibile smoothando tutto facendo interpolazioni tra un elemento e l’altro. Questo smooth è possibile sia metterlo che toglierlo di default.