Modelli numerici elasto-plastici

In questa analisi vengono introdotti alcuni modelli computazionali per l’analisi agli elementi finiti di comportamenti di hardening isotropo e cinematico. Viene analizzata la risposta ad un carico ciclico per alcuni acciai inossidabili. Queste simulazioni sfruttano un elemento finito quadratico lineare con comportamento elasto plastico con hardening isotropo e cinematico. Tali modelli risultano fondamentali per analizzare la risposta materiale a seguito di un carico ciclico ripetuto.

Background

Il comportamento materiale più noto è il comportamento elastico in cui la fenomenologia osservata è tale per cui se consideriamo tensione e deformazione, applichiamo una deformazione e osserviamo una certa tensione e una volta rimossa si torna allo stato iniziale.

Elasto-plasticità

Un’estensione è data dalla possibilità di portare in conto anche un modello di plasticità per il quale, superata una certa tensione σ_y si attiva un meccanismo di plasticità. Ovvero la tensione verrà governata anche dalla deformazione plastica:

\sigma =\mathbb{C}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)

Questo richiede, per poter essere risolto, la conoscenza dell’evoluzione dei meccanismi di plasticità e quindi della variazione di ε_p . Ovvero richiede l’introduzione di una legge evolutiva che permetta di esprimere i fenomeni di plasticità e allo stesso tempo rispettare le leggi termodinamica. per farlo si definisce una funzione di snervamento, espressa nello spazio delle tensioni, al cui interno saremo in regime di elasticità mentre al di fuori si saranno innescati dei meccanismi plastici.

Una scelta di energia di deformazione che soddisfa al contempo il comportamento materiale e le leggi della termodinamica, sarà:

\psi\left(\underline{\varepsilon}, \underline{\varepsilon}_{p}\right)=\frac{1}{2}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)^{t} \mathbb{C}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)

Allora l’importanza risiede nell’introduzione di una legge evolutiva che soddisfi i requisiti termodinamici e la scelta tipica consisiste in una legge associativa. Ovvero una legge tale da essere assciata alla funzione di snervamento tramite un moltiplicatore incognito, dove di sceglie f come una funzione convessa.

 \dot{\varepsilon}_{p}=\dot{\lambda} \frac{\partial f}{\partial \sigma}

Questo dice che la deformazione plastica è un vettore con direzione normale alla superficie di snervamento. Seguendo l’evoluzione della tensione questa dovrà, una volta raggiunto il valore di snervamento, rimanere a questo livello di tensione quindi una volta attivata la plasticità rimarrà sulla superficie muovendosi muoverò sul bordo.

Allora questo, dal punto di vista computazionale, intoduce una nuova funzione (oltre l’equilibrio elastico) che dovrà essere risolta dal medesimo stato tensionale tramite un procedimento iterativo.

Hardening

Un modello ancora più sofisticato include la possibilità di incrudimento, ovvero la deformazione plastica non sarà più lineare ma avrà una certa pendenza, espressa appunto dal modulo di hardening.

Dunque un primo modello di incrudimento è quello isotropo dove, rimosso il carico, si torna indietro lungo lo stesso ramo elastico. Rimarrà quindi una certa deformazione plastica. Continuando il processo di carico nel senso opposto si attiverrà la plasticità opposta. L’incrudimento avrà determinato una maggior resistenza del materiale in tutte le direzioni.

Invece nell’incrudimento cinematico, variando la direzione potrà variare il nuovo livello di tensione di snervamento a cui si attiva la plasticità.

Nel caso ideale si misura lo stato tensionale verificando che sia all’interno della funzione di snervamento

f(\underline{\sigma})=\|\underline{\sigma}\|-\sigma_{y_{0}}

Mentre nell’incrudimento isotropo sarà la funzione di snervamento a variare:

f(\underline{\sigma}, q)=\|\underline{\sigma}\|-\left(\sigma_{y_{0}}-q\right)

E un’energia libera funzione di un campo scalare α che porta in conto il modulo di Hardening, H, ovvero:

\psi\left(\underline{\varepsilon}, \underline{\varepsilon}_{p}, \alpha\right)=\frac{1}{2}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)^{T} \mathbb{C}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)+\frac{1}{2} H \alpha^{2}

Invece nell’incrudimento cinematico invece verrà semplicemente traslata:

f(\underline{\sigma}, \underline{q})=\|\underline{\sigma}-\underline{q}\|-\sigma_{y_{o}}

E un’energia libera analoga alla precedente ma che porta in conto una grandezza tensoriale:

\psi\left(\underline{\varepsilon}, \underline{\varepsilon}_{p}, \underline{\alpha}\right)=\frac{1}{2}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)^{T} \mathbb{C}\left(\underline{\varepsilon}-\underline{\varepsilon}_{p}\right)+\frac{1}{2} \underline{\alpha}^T \mathbb{H} \underline{\alpha}

Questo porta a considerare la tensione di incrudimento che prende il nome di backstress a cui viene riferita la legge associativa:

\underline{q}=-\frac{\partial \psi}{\partial \underline{\alpha}}=-\mathbb H\underline{\alpha}
\\\:\\
\dot{\alpha}=\dot{\lambda} \frac{\partial f}{\partial q}

Formulazione FEM

Questo modello di hardening deve fondersi con la formulazione agli elementi finiti per risolvere il problema dell’equilibrio elastico. Partendo dalla forma debole si va a risalire alla geometria parente e intrdurre l’interpolazione sulle coordinate nodali.

\int_{V} \underline{\sigma} \cdot \delta \underline{\varepsilon} d V-\int_{V} \underline{f} \cdot \delta \underline{u} d V-\int_{A} t \cdot \delta u d A=0

Ovvero si risale all’i-esima componente del residuo

[R]_i=\sum_{g} B_{I}^{T}(\xi, \eta) \mathbb{C}\left(\sum_{k} B_{k}\left(\xi_{g}, \eta_{g}\right) \underline{u}_{k}-\left.\underline{\varepsilon}\right|_{g}\right) J_{e, g} d V

Questo richiede le coordinate nodali e la deformazione plastica regolata dalla flow rule, variabile nel tempo, per la quale si ricorre ad una discretizzazione in passi temporali.

\underline{\varepsilon}_{p}=\underline{\varepsilon}_{p_{n}}+\Delta t \cdot \underline{\dot{\varepsilon}}_{p_{n}}

Cioè:

\dot{\varepsilon}_{p}=\frac{\dot{\varepsilon}_{p}-\underline{\varepsilon}_{p_{n}}}{\Delta t}=\dot{\lambda} \cdot \underline{n}_{\sigma}=\frac{\lambda-\lambda_{n}}{\Delta t} \underline{n}_{\sigma}

E da questa relazione è possibile considerare le leggi evolutive rappresentanti di una plasticità tempo indipendente che non risulta legata alla velocità con cui viene applicato il carico:

\left\{\begin{array}{l}\underline{\varepsilon}_{p}=\underline{\varepsilon}_{p_{n}}+\left(\lambda-\lambda_{n}\right) \underline{n}_{\sigma}(\underline{\sigma}, q) \\ \underline{\alpha}=\underline{\alpha}+\left(\lambda-\lambda_{n}\right) \underline{n}_{q}(\underline{\sigma}, q) \\ f(\underline{\sigma} \cdot \underline{q})=0\end{array}\right.

Ovvero queste relazioni costituiscono un secondo sistema algebrico che dovrà essere risolto insieme a quello rappresentante l’equilibrio elastico.

R(\underline{u}, \underline{h})=0
\\\:\\
\underline{Q}(\underline{u}, \underline{h})=0

Nel risolverlo si sfrutta un algortimo predittore-correttore che calcola le tensioni assumendo inizialmente una deformazione elastica e varificando che lo stato tensionale sia entro i limiti stabiliti dalla funzione di snervamento. Se non viene soddisfatta la condizione si va a corregere la predizione trovando la soluzione dal sisitema Q che riporta sul bordo della funzione di snervamento. Quindi si trova un nuovo stato tensionale che dovrà soddisfare anche l’equilibrio sul sistema R. Si procede con passi iterativi soddisfando entrambi i sisitemi algebrici.

Campagna di simulazioni

La campagna di simulazioni viene condotta in AceFEM. Si concentra su due differenti strutture. Si considera dapprima un patch test, per verificare la correttezza del codice e per avere un riscontro di come la variazione delle caratteristiche di hardening influenzi la risposta strutturale e l’andamento delle tensioni interne. Successivamente si sposta l’analisi su una piastra (L x H) con una fessura di dimensione variabile. L’analisi parte utilizzando un elemento finito Q1, quadrilatero con funzioni di forma lineari, con comportamento elastico plastico ed incrudimento isotropo. Questo verrà poi trasformato in un elemento capace di riprodurre un comportamento elasto-platico con hardening di tipo cinematico.

Q1EPS

Per implementare un comportamento i hardening cinematico è necessario considerare che la variabile α diventerà una variabile vettoriale e la legge evolutiva varierà tenendo conto della traslazione della funzione di snervamento nello spazio delle tensioni.

\underline q=-{d\psi \over d\underline \alpha}
\\\:\\
f(\underline \sigma,\underline q)=\sqrt{{3\over 2}\|\underline \sigma_{\operatorname{dev}}-\underline q}-\sigma_{yo}

Quindi sarà calcolato in modo diverso anche il sistema per l’algoritmo predittore correttore:

\bold Q_{\alpha}=\underline \alpha-\underline \alpha_n-(\lambda-\lambda_n){\partial f\over \partial \underline q}

Questo si rispecchia nell’introdurre nuove variabili di laboro, aggiornare la definizione dell’energia di deformazione portando in conto il prodotto interno e nel ridefinire la legge evolutiva.

SMSFreeze[\[Alpha], {{\[DoubleStruckH]g[[4]], \[DoubleStruckH]g[[
    6]]}, {\[DoubleStruckH]g[[6]], \[DoubleStruckH]g[[5]]}}, 
 "Symmetric" -> True]
 \[Lambda] = \[DoubleStruckH]g[[7]];

\[CapitalPsi]e \[DoubleRightTee] \[Lambda]e/ 2 (Tr[\[DoubleStruckCapitalD]e])^2 + \[Mu]e Tr[\\[DoubleStruckCapitalD]e . \[DoubleStruckCapitalD]e] +  1/2 H Tr[Transpose[\[Alpha]] . \[Alpha]];

fg \[DoubleRightTee] SMSSqrt[3/2 Tr[Transpose[\[Sigma]dev - q] . (\[Sigma]dev - q)]] - (\[Sigma]YO);

Inoltre viene aggiornato il sistema per l’algoritmo predittore correttore portando in conto i nuovi valori:

\[Alpha]n = {{\[DoubleStruckH]gnIO[[4]], \[DoubleStruckH]gnIO[[
   5]]}, {\[DoubleStruckH]gnIO[[5]], \[DoubleStruckH]gnIO[[6]]}};
\[Lambda]n = \[DoubleStruckH]gnIO[[7]];

SMSFreeze[q, -SMSD[\[CapitalPsi]e, \[Alpha], "Symmetric" -> True], 
  "Symmetric" -> True];

\[DoubleStruckCapitalQ] = {\[DoubleStruckCapitalQ]\[Epsilon][[1, 
    1]], \[DoubleStruckCapitalQ]\[Epsilon][[2, 
    2]], \[DoubleStruckCapitalQ]\[Epsilon][[1, 2]], Qq[[1, 1]], 
   Qq[[2, 2]], Qq[[1, 2]], Q\[Lambda]};

Questo permette di generare un nuovo elemento finito.

• Q1EPS, quadrilatero ordine 1 con comportamento elasto plastico e hardening isotropo
• Q1EPS2, quadrilatero ordine 1 con comportamento elasto plastico e hardening cinematico

Caratteristiche di hardening

È possibile considerare un acciaio inossidabile viste le ottime proprietà meccaniche e il grande campo di applicazione negli apparati medicali, strumenti chirurgici, protesi e molto altro. È proprio la loro composizione chimica con il cromo presente almeno al 10.5 % e un ridotto contenuto di carbonio sotto al 1.2 % a renderli fortemente resistente alla corrosione e a portare alla generazione di microstrutture che tipicamente non sono presenti negli acciai e portato a tali proprietà meccaniche. Chiaramente al variare della composizione dell’acciaio varieranno diverse proprietà meccaniche compresa la risposta strutturale. Possiamo quindi considerarne dei parametri medi e affrontare uno studio parametrico. È utile ricordare che possiamo definire il modulo di hardening a partire dal carico di snervamento e dal carico critico come:

H={\sigma_u-\sigma_y\over {\varepsilon_{u,\%}\over 100} - \varepsilon _y}

Gli acciai inossidabili sono tipicamente più duttili e presentano la possibilità di accogliere deformazioni plastiche molto grandi. Queste proprietà plastiche si riflettono nel hardening e in un incremento della tensione di snervamento, tensione ultime e della durezza. Viene quindi scelto un modulo di hardening variabile tra 600 e 1600 MPa e una tensione di snervamento tra 175 e 450 MPa.

Risposta strutturale e hardening isotropo

Patch test

Il patch test è un test che viene effettuato su un dominio molto semplice discretizzato in modo irregolare ma con una mesh molto lassa. Questo permette, sfruttando condizioni al bordo di cui è nota la soluzione analitica, di testare nuovi algoritmi. Inoltre, permette anche di percepire come le variazioni di alcuni parametri materiali influenzino il comportamento strutturale andando ad agire su risposte materiali di per se più semplici.

In questo caso viene analizzato un dominio quadrato discretizzato in 5 elementi. Si impone una trazione monoassiale per la quale otterremo una tensione costante. Andremo allora ad applicare una condizione di carico variabile nel tempo che alterna compressione e trazione tramite il moltiplicatore λ.

Si può allora analizzare l’andamento delle tensioni nel tempo osservando esattamente un comportamento con hardening isotropo.

Tale comportamento ci mostra come la tensione σ_yy cresce linearmente sia con l’aumentare di σ_y che del modulo di hardening.

Piastra

Considerando invece una piastra L x H con una fessura a metà altezza di lunghezza pari ad L/4 è possibile osservare come varia la sua risposta strutturale. Andando ad applicare un carico di trazione verticale si osserva un allargamento della cricca con una concentrazione delle tensioni proprio a ridosso dell’apice della cricca. Ovvero ci sarà una separazione.

Quindi è possibile analizzare la deformazione al variare delle caratteristiche di hardening.

Tensione di snervamento

Si osserva una riduzione dell’allargamento della cricca al crescere della tensione di snervamento. Un materiale con una tensione di snervamento più alta presenterà un minor spostamento. I risultati in figura considerano un coefficiente di Hardening medio.

Coefficiente di Hardening

Al crescere del coefficiente di hardening si osserva una riduzione dello spostamento verticale. Questo è dovuto al fatto che considerando le stesse condizioni di carico ci sarà un materiale più rigido che risponderà con una minore deformazione plastica.

Si può osservare l’effetto complessivo sia di una variazione della tensione di snervamento che del coefficiente di hardening con un contour plot. La figura seguente mostra come al crescere della tensione di snervamento o del modulo di hardening si ha una deformazione plastica residua sempre più piccola. In particolare, in figura viene plottato lo spostamento verticale del bordo superiore della cricca al variare dei parametri di hardening.

Variazioni

Considerando il caso con i parametri di hardening in condizioni medie è possibile andare a vedere come il moltiplicatore del carico influenza lo spostamento al variare del tempo.

Un’altra variazione importante è data dalle tensioni σ_yy andandole a leggere a metà piastra verticalmente con una coordinata che varia dalla fine della cricca fino al bordo destro.

Fattore di intensificazione degli sforzi

Considerando i parametri di hardening fissati ad un valore medio è possibile osservare cosa succede al variare della lunghezza della cricca. Un parametro molto interessante da osservare è il fattore di intensificazione degli sforzi. Con la presenza della cricca si crea una distribuzioni di strees all’interno della piastra differente da quella che si avrebbe considerando una piastra omogenea. Precisamente si osserva lo stress aumentare avvicinandosi alla cricca.

Vengono analizzati differenti casi in cui la cricca ha lunghezza:

• L/8
• L/4
• L/2
• 3L/4

Hardening cinematico

I risultati ottenuti finora considerano l’elemento finito Q1EPS, ovvero un comportamento elasto plastico con hardening di tipo isotropo. Utilizzando invece l’elemento Q1EPS2 è possibile analizzare un hardening di tipo cinematico. Il comportamento è subito evidente andando ad effettuare la simulazione su un patch test.

Questo comportamento rispecchia maggiormente ciò che accade realmente nei metalli dove la tensione di snervamento a trazione aumenta a spese della tensione di snervamento a compressione. Questo è noto come effetto Bauschinger. Ovvero la tensione di snervamento a compressione di un materiale portato a snervamento per trazione è inferiore a quella del materiale originale. Si ottiene lo stesso lo stesso effetto con una storia di carico inversa. Questo effetto è quanto testimoniato dell’asimmetria del ciclo di isteresi di un materiale duttile.

Risposta macroscopia

Evidenze di un hardening di tipo cinematico appaiono anche andando ad effettuare la simulazione con un carico variabile nel tempo secondo una legge triangolare sulla piastra con la cricca o su una membrana di Cook.

La variazione della tensione di snervamento a compressione si rispecchia in un maggior effetto di snervamento nella fase di compressione.

Nell’incrudimento isotropo il materiale ricorda la precedente fase di carico e lo snervamento successivo è rappresentato da un punto a tensione di snervamento maggiore ma aumentata in modo omotetico. Invece, nell’incrudimento isotropo, si ha la traslazione della superficie limite. Questo backstress si rispecchia in una riduzione della tensione di snervamento su una direzione differente di quella del primo ciclo di carico.

È possibile esportare anche l’animazione per vedere come varia lo spostamento e la deformazione della membrana di Cook lungo la storia di carico (passi temporali).

Codice

Il codice è disponibile alla repository GitHub: https://github.com/mastroalex/elasto-plastic

In particolare:

• Variazioni delle proprietà di hardening ../final_test_ripetuti_11_passi/
• Variazioni delle proprietà della cricca ../final_test_vari/
• Hardening cinematico ../final_kinematic/

Riferimenti

• von Mises Material Properties – Autodesk Knowledge
• Mechanical properties of structural stainless steels – Tylek, I and Kuchta, K
• Materiali per l’ingegneria civile ed industriale – Callisters, W and Rethwisc, D
• Scienza e tecnologia dei materiali – Smith, W and Hashemi, J
• Stainless Steels for Medical and Surgical Applications – Nutt, M