Come ovvio, una volta abbondonata la soluzione elastica lineare, la soluzione delle equazioni di equilibrio della struttura si traduce in un set di equaioni non lineari (tralasciamo qui la fonte di tali non linearità) nelle quali, nella stragrande maggiornza dei casi la fonte delle non-linearità risiede nella variazione della rigidezza della struttura e quindi nella matrice di rigidezza K della struttura stessa [i].
Un modo per risolvere il problema consiste nell'approssimare la curva di carico della struttura attraverso una suddivisione in passi di carico e monitornando all'interno di ogni singolo passo le variazioni della rigidezza della struttura.
P= Σ dPi
U= Σ dUi
K dUi= dPi
Il modo con cui le equazioni di equilibrio incrementale vengono formulate ed il modo con cui viene conseguito l'equilibrio all'interno di ogni passo di carico caraterizza il tipo di analisi eseguita. Tralasciando i metodi che fanno uso della matrice di rigidezza secante (per altro molto utili in altri campi quali il calcolo del moltiplicatore di collasso di sezioni in c.a., ecc...) ed i metodi single-step in cui la matrice di rigidezza viene formulata una sola volta all'inzio del passo di carico, il metodo di analisi più utilizzato, ed impiegato nel codice di calcolo non-lineare WinStrand, è noto come metodo iterativo.
Utilizzando questo metodo, al solito, il carico esterno viene suddiviso in n parti o 'passi di carico' ed all'interno di ogni passo di carico si cerca iterativamente, utilizzando le propietà meccaniche aggiornate della struttura, di conseguire l'equilibrio. Pertanto gli spostamenti all'interno di ogni singolo step di carico i all'interno di ogni singola iterazione j potranno scriversi come:
Δi= Δi-1 + Σj dΔij , con j= 1, 2, ..., m
dove il campo di spostamenti dΔij è ottenuto risolvendo, iterativamente, un'equazione del tipo:
K-1 dΔij = dPij + Rij-1
essendo
K-1 : matrice di rigidezza della struttura valutata in configurazione deformata tenendo conto dello stato di sollecitazione presente nelle membrature
Rij-1 : vettore delle forze squilibrate ottenuto controllando l'equilibrio fra le forze esterne applicate e le forze presenti nella struttura oggetto di analisi
ossia
Rij-1 = Pij-1 - Fij-1
Ai fini di una più generale descrizione del metodo è conveniente fare le seguenti posizioni:
per mezzo delle quali è intuitivo comprendere lo sviluppo del processo itereativo in un ciclo di calcolo.
In sostanza l'algoritmo risolutivo, all'i-esimo passo di carico, si può sintetizzare nei seguenti passaggi:
- Calcolo della matrice di rigidezza tangente noti che siano lo stato deformativo e di sollecitazione della struttura al passo i-1
- Calcolo dell'incremento di spostamenti dΔij dovuto all'incremento di carico dPij ed alle azioni squilibrate dRij-1 (nulle al primo ciclo del processo iterativo)
- Determinazione dello stato di deformazione e sollecitazione nell'intera struttura
- Determinazione delle azioni squilibrate dRij
- Controllo del raggiungimento della convergenza (valutata sugli spostamenti/azioni squilibrate/ecc...). Se la convergenza richiesta è stata raggiunta stop, altrimenti si ritorna al passo 1 per un nuovo ciclo.
Questo modo di procedere è noto come algoritmo di Newton-Raphson e prevede l'aggiornamento della matrice di rigidezza tangente ad ogni ciclo di iterazione j all'interno del passo di carico i. Questa operazione risulta essere computazionalmente tra le più onerose e pertanto, oltre all'algoritmo originale, viene comunemente impiegato con successo il cosiddetto algoritmo di Newton-Raphson modificato impiegando il quale la matrice di rigidezza tangente viene assemblata e fattorizzata una sola volta all'inizio del passo di carico i per poi essere reimpiegata nelle successive m iterazioni interne allo stesso passo di carico i.
Questo modo di operare, in generale, ha una velocità di convergenza minore rispetto all'algoritmo originale ma ciò nonostante, data la natura dei problemi usualmente trattati nell'ingegneria civile, consente di compendiare due esigenze e cioè:
- Una buona velocità di convergenza assicurata solitamente dalle buone caratteristiche della matrice di rigidezza (prevalenza della diagonale principale rispetto ai termini fuori diagonale, buon condizionamento della matrice di rigidezza stessa) perlomeno per le usuali strutture ed al di fuori delle zone critiche (punti di criticità per instabilità di IIa specie, ecc..)
- Una economicità computazionale derivante dal non dover continuamente rivalutare le matrici di rigidezza dei singoli elementi e soprattutto riassemblare e fattorizzare la matrice di rigidezza globale.
Nella versione attuale del condice di calcolo non lineare è implementata, per le ragioni sopra menzionate, questa seconda versione del metodo di Newton-Raphson. La matrice di rigidezza tangente viene assemblata e fattorizzata una sola volta all'inizio dell'i-esimo step di carico in forma LDLT ed in output vengono forniti i termini diagonali kii min/max pre-fattorizzati, post-fattorizzati oltre al decadimento massimo dei termini diagonali.
Per quanto riguarda il singolo step di carico, attualemnte è implementato il cosiddetto 'Load Control Method'. Questo metodo consiste nell'applicare al primo ciclo iterativo (j=1) del singolo passo di carico l'incremento di carico pertinente mentre nei successivi cicli iterativi (j>1) vengono applicati i soli squilibri fino al raggiungimento dell'equilibrio.
è evidente che questo metodo pone grossi problemi in vicinanza dei punti critici in quanto non consente di intervenire ne sul livello di carico ne sull'entità degli spostamenti.
Ciononostante, questa tecnica, da sempre base dei metodi iterativi, nelle situazioni usuali garantisce buoni risultati senza dover intervenire dall'esterno fornendo ulteriori parametri di convergenza quali il massimo campo di spostamenti accettabile, il lavoro massimo compiuto dai carichi esterni, ecc..
Nelle successive revisioni del software sono pianificate le implementazioni di uno o più di questi metodi basati sul cosidetto 'Arc lenght method' [ii].
Infine una nota sui metodi attraverso i quali valutare la raggiunta convergenza del ciclo iterativo. Considerando il fatto che il campo di deformazioni manifestato dalla struttura è il cuore del meccanismo di aggiornamento, la letteratura tecnica suggerisce vari metodi per valutare la raggiunta convergenza quali:
Norma degli spostamenti modificata:
|e|= 1/N Σk=1N |dΔk/dΔref|
Norma euclidea modificata:
|e|= 1/N Σk=1N (dΔk/dΔref)2
Norma massima:
|e|= max1<=k=N{ |dΔk/dΔref|}
In queste equazioni N è il numero di incognite spostamento del problema, dΔk è la generica componente di spostamento incrementale e dΔref è lo spostamento di riferimento (traslazione/rotazione) assunto come parametro di riferimento (solitamente pari alla massima componente di traslazione/rotazione presente nel vettore degli spostamenti totali). In tutti questi casi la norma dello spostamento |e| va confrontata con il massimo errore ammesso e tipicamente compreso fra 10-2 e 10-6 che và predefinito dall'utilizzatore. Attualmente, nel codice viene impiegata la Norma Euclidea Modificata. Future revisioni del software consentiranno l'impiego di altri criteri di convergenza.
[i] In ingegneria civile si hanno pochi casi di sistemi in cui la fonte di non linearità è dovuta ai carichi. Fra questi ricordiamo le azioni di precompressione nei conci dei ponti 'Dividag' e similari soprattutto in fase transitoria, ecc.
[ii] M.J. Clarke & G.J. Hancook "A study of incremental-Iterative
strategies for nonlinear analysis" Int J. Numerical Methods In Engr. Vol 29 1990 pp. 1365-1391
O. Zienkiewicz. "Incremental
displacement in Non Linear Analysis" Int. J. Num. Meth. In Egrr. Vol 3, 1971, pp 587-592
M.A. Crisfield "NonLinear
Finite Element Analysis of solid and structures. Volume I essentials Volume II Advanced Topics" John Wiley and Sons
New York 1991-1997
K.J, Bathe and E.L.Wilson "Numerical methods in finite element nalysis" Prentice-Hall Englewood
cliffs N.J. 1976