Microsoft Word - numero 11 art 4 M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 36 Modello coesivo per l’avanzamento di fratture mediante rilascio nodale di strutture discretizzate con elementi finiti M. Minotti, P. Salvini Università di Roma “Tor Vergata”, Dipartimento di Ingegneria Meccanica,Via del Politecnico, 1 – 00133 Roma minotti@ing.uniroma2.it, salvini@uniroma2.it G. Mannucci Centro Sviluppo Materiali S.p.A. di Roma,Via di Castel Romano, 100/102 – 00100 Roma g.mannucci@c-s-m.it RIASSUNTO. La simulazione numerica della propagazione di una frattura in MODO I, viaggiante ad elevata velocità in un acciaio a comportamento duttile è realizzata attraverso un modello coesivo che governa la distribuzione delle forze di rilascio nodale. Come noto, la ricerca di un valore di tensione all’apice non ha alcun senso nel caso elastico; infatti, la tensione può essere valutata solo mediante fattori di campo. Nel caso elasto- plastico, incrudimento e softening di origine geometrica o legato al progressivo danneggiamento influiscono decisamente sull’andamento esponenziale del campo tensionale. È possibile comunque individuare un valore di riferimento, di entità finita, mediante estrapolazione delle tensioni elasto-plastiche nella zona di inizio softening della frattura. Tale grandezza può essere presa come fattore di riferimento per il calcolo delle forze di rilascio coesive e quindi dell’energia dissipata. Nel lavoro viene discusso come determinare, dal campo di tensione elasto-plastico locale, il valore che governa la zona coesiva al variare del T-stress. ABSTRACT. High speed MODE I crack growth in elastic-plastic materials, involving large scale plasticity and dynamic effects connected to rapid propagation, is faced through a cohesive model to compute force nodal release. As it is well known, the search of a maximum stress in the crack tip is not worthwhile in the elastic case, only field factors may be used instead. In the elastic-plastic case, hardening and softening due to geometric reasons or due to local damaging have a deep influence on the exponential behaviour of the measured stress field. However, it is possible to introduce a reference value (finite), which is the stress at the beginning of the softening zone. This entity is taken as a scale factor for cohesive force evaluation that is necessary for the accounting of the dissipated energy. In the present work a procedure to compute this reference value by stress- strain field is presented, for a variety of T-stress conditions. PAROLE CHIAVE. Elementi finiti, meccanica della frattura, modello coesivo, forze di rilascio nodale. INTRODUZIONE li strumenti messi a disposizione dalla meccanica della frattura hanno subito importanti sviluppi negli anni, per cui si dispone oggi di tecniche di progettazione molto affidabili per l’individuazione della criticità di una ampia classe di difetti [1]. In alcuni casi tuttavia, si considerino ad esempio strutture soggette a pressione interna, non basta conoscere le sole condizioni di innesco, ma occorre anche prevedere il comportamento evolutivo e quindi dinamico del difetto in modo da controllare la combinazione svuotamento-propagazione della frattura che, in buona sintesi, governa G http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 37 la pericolosità dell’evento. Questo genere di previsioni è molto più complesso e non gode ancora di tecniche di valutazione affermate e applicabili ad una vasta tipologia di casi [2]. La propagazione può essere controllata, estremizzando, da un regime di comportamento duttile o fragile del materiale in fratturazione; essa può estendersi anche per tratti piuttosto consistenti e quindi in condizioni di quasi stazionarietà [3]. Nel caso di acciai tenaci, se la driving force si mantiene sufficientemente elevata durante la propagazione, si assiste ad un comportamento fortemente duttile della frattura, caratterizzata da una velocità di avanzamento molto inferiore rispetto a quella di Raleigh [4]. La propagazione è associata ad un consistente fronte plastico che precede l’apice della cricca. In tali particolari condizioni i parametri di controllo usuali per lo studio della meccanica della frattura mostrano tutti i loro limiti di utilizzo. La frattura presenta fenomeni locali molto complessi (tunneling, necking, continue modificazioni della direzione del piano di fratturazione) per cui risulta decisamente difficile pretendere una rappresentazione locale dello stato tensionale. Una metodologia di analisi alternativa, per di più estremamente vantaggiosa dal punto di vista computazionale, è resa disponibile dall’utilizzo di modelli coesivi o di progressivo rilascio dei lembi di frattura [5]. Limitatamente alla frattura duttile, è fondamentale il modo in cui si simula la dissipazione di energia in prossimità dell’apice della cricca. Una prima possibilità consiste nell’introdurre i complessi legami costitutivi che considerano il danneggiamento progressivo del materiale [6]; tale approccio, sebbene utile a prevedere molti aspetti della frattura duttile, non è adatto allo studio della propagazione della frattura ove si presentino deformazioni plastiche estese e di notevole entità associati ad altrettanto elevati strain-rate [7]. D’altra parte il limite imposto alla dimensione della cella elementare condiziona fortemente l’efficienza del calcolo nel caso di strutture di considerevoli dimensioni. A questi svantaggi si aggiungono le difficoltà numeriche derivanti dall’impiego di tali modelli nei codici di calcolo con risoluzione esplicita, opportunamente utilizzati per l’analisi di fenomeni in cui sono presenti elevate forze dinamiche associate a grandissimi spostamenti. Per tali motivi, i legami costitutivi che affrontano il danneggiamento localizzato sono consigliabili solo per l’analisi di strutture di dimensioni modeste e per atti di propagazione contenuti (i.e. nella fase di innesco e di formazione stabile della frattura). L’alternativa è rappresentata dall’uso di metodi locali di computazione basati sulla tecnica del rilascio controllato dei nodi; essi sono in grado di affrontare il fenomeno della propagazione delle fratture, in condizioni stabili o instabili, anche per lunghe distanze, pur richiedendo osservazioni sperimentali di fratturazione simili a quelle in studio. Questo approccio è di particolare interesse nella modellazione della frattura di strutture a guscio ove le sollecitazioni agenti sui nodi di frattura possono essere assimilate semplicemente a forze. In tal caso si utilizzano convenientemente gli elementi di tipo shell, dove l’effetto dello spessore viene portato in conto mediante alcuni accorgimenti. Il comportamento della cricca durante la propagazione può essere valutato mediante l’impiego di opportuni criteri, basati su parametri caratteristici della frattura duttile/fragile, estesi alle condizioni dinamiche. Da questo punto di vista il Jintegral può essere adatto allo scopo, ma presenta alcuni svantaggi dal punto di vista computazionale che ne rendono complessa l’applicazione, soprattutto nei codici di calcolo di tipo esplicito [8]. Secondo recenti studi il CTOA (Crack Tip Opening Angle) è un parametro che può descrivere correttamente le condizioni di propagazione della frattura, soprattutto in campo duttile, in quanto si mantiene costante in caso di propagazione stazionaria e non richiede l’adozione di elementi finiti particolari, o di prefissate dimensioni. L’energia specifica dissipata durante la propagazione rappresenta un altro parametro di estremo interesse nell’ambito dello studio della frattura, la cui efficacia è indipendente dalla natura duttile o fragile della propagazione. In particolare l’Essential Work of Fracture (Wf) [9] costituisce un termine direttamente legato alla capacità del materiale di deformarsi plasticamente e fornisce l’entità dell’energia dissipata, nella regione di materiale prospiciente l’apice della cricca, nota come Fracture Process Zone (FPZ). In essa si sviluppa il cedimento del materiale con formazione di nuove superfici di frattura [1,3,4]. È nella zona di processo che i modelli coesivi simulano il progressivo indebolimento strutturale attraverso l’applicazione, sui nodi contenuti nella FPZ, di forze di rilascio la cui entità è funzione dello stato tensionale agente e della distanza dall’apice della frattura. Nei codici di calcolo agli elementi finiti con formulazione esplicita, risulta assai difficoltoso valutare l’istantaneo campo di tensioni agenti, soprattutto a causa della dinamica di alta frequenza che caratterizza tale formulazione. È necessario quindi definire una procedura che consenta di determinare con accuratezza l’entità della tensione di apertura (MODO I) in prossimità dell’apice, al fine di tarare correttamente il modello coesivo ed ottenere una simulazione corretta del fenomeno della propagazione. STATO TENSIONALE DEL FRONTE DI CRICCA on l’ausilio degli elementi finiti è stato condotto uno studio approfondito sul campo di tensioni in prossimità dell’apice di una cricca sollecitata prevalentemente in MODO I ed in condizioni di elasto-plasticità. In particolare si è indagata l’influenza sul valore delle tensioni ortogonali al piano di frattura, esercitata da sollecitazioni imposte parallelamente all’asse della cricca che si vanno a sovrapporre a quella autoindotta dalla geometria del pezzo (T-stress), C http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 38 utilizzando il criterio di Von Mises per il calcolo delle tensioni equivalenti. Il modello adoperato per l’analisi consta di una piastra quadrata di notevoli dimensioni, lato 2 m e spessore 1 mm, che presenta una cricca centrale di estensione 20 mm (Fig. 1), cosicché lo stato di tensione all’apice della frattura non risulti influenzato da significativi effetti di bordo. Figura 1: Geometria e disposizione dei cairchi utilizzati per il modello agli elementi finiti. La piastra è sollecitata simmetricamente con un carico distribuito S sui bordi paralleli ai lembi della cricca e T perpendicolarmente. Le entità dei suddetti carichi sono espresse in rapporto alla tensione di snervamento 0 . L’elemento finito utilizzato è di tipo piastra a 4 nodi adatta ad analisi strutturali in campo elasto-plastico e grandi spostamenti. Il materiale utilizzato è un acciaio di grado API X65 con valori di Ramsber-Osgood equivalenti pari a 13n  e 0.56  . Figura 2: Variabilità di σyy (MPa) in stato di tensione piana e deformazione piana a parità di T agente. Le simulazioni sono state effettuate in regime di grandi spostamenti sia in stato di tensione piana che di deformazione piana, per diversi valori di tensione di apertura S, ed a parità di S per diversi valori del rapporto T. In Fig. 2 si possono apprezzare i campi di isotensione nelle due condizioni di tensione piana e deformazione piana per due valori opposti del campo aggiuntivo T. Nelle Figg. 3a e 3b sono riportati gli andamenti yy relativi al caso di tensione piana e deformazione piana della porzione di materiale plasticizzata posta di fronte l’apice di frattura, avendo a parametro diversi valori della sollecitazione T. Dal confronto tra le due è possibile rilevare come nel caso di tensione piana, carichi paralleli di trazione influenzino marginalmente il valore della tensione in direzione ortogonale; al contrario se la sollecitazione parallela è di compressione lo stress yy aumenta sensibilmente. Effetti analoghi ma invertiti si possono osservare in regime di deformazione piana, dove solo carico di trazione parallelo alla cricca determina un innalzamento consistente della componente ortogonale della http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 39 tensione. Per quanto attiene l’andamento delle curve di intensificazione, esse presentano un trend di tipo esponenziale che è governato sia dai parametri di Ramsber-Osgood che dallo stato tensionale nominale presente attorno all’apice della frattura [10]. In ogni caso il valore ultimo risulta limitato, ossia finito. Esso assume un’importanza notevole in quanto a partire dal suo valore si sviluppano le forze coesive di rilascio progressivo nella zona di processo. Come si vede dalle Figg. 3a e 3b, i valori della tensione riportati non coinvolgono anche il vero e proprio apice, ciò è dovuto al mesh refinement che impone l’esclusione del primo nodo vincolato. Nel nostro caso ciò non costituisce una limitazione significativa in quanto il valore di riferimento viene comunque determinato mediante estrapolazione e non da lettura diretta. (a) (b) Figura 3: Valori di σyy in tensione piana (a) e in deformazione (b) per diverse tipologie di carico. IL MODELLO COESIVO ’analisi di propagazione della frattura nei materiali duttili per mezzo degli elementi finiti, deve necessariamente tener conto che davanti all’apice di frattura si verifica una forte crescita delle tensioni, accompagnata da estese deformazioni in campo plastico. In queste condizioni il distacco dei lembi non avviene in modo istantaneo ma segue i seguenti stadi: formazione di una superficie libera intorno un difetto microscopico, crescita della cavità attorno al difetto favorita dalla deformazione plastica, coalescenza tra le cavità contigue, frattura. Per simulare correttamente su scala macroscopica l’evolversi dei fenomeni succitati si considera una zona di processo (Fracture Process Zone − FPZ), dove hanno luogo i fenomeni coesivi che si accompagnano alla formazione delle nuove superficie di frattura. Per lo studio di frattura in cui lo spessore della stessa è molto inferiore rispetto alle altre dimensioni dell’oggetto fratturato, è lecito considerare una FPZ monodimensionale e quindi caratterizzabile con la sola lunghezza Δ chiamata distanza di estinzione. La valutazione di Δ non segue considerazioni particolari e spesso è fissata in relazione alla dimensione degli elementi di frattura impiegati. Questa assunzione, sebbene conduca ad una dipendenza dalla mesh utilizzata, può essere accettabile nel caso di frattura fragile ove la zona di danno, localizzata sull’apice, è molto piccola e la dimensione degli elementi è la minima utilizzabile per rilasciare le forze. Al contrario la scelta di Δ diventa fondamentale nel caso di frattura duttile quando il cedimento del materiale può interessare un’ampia regione in prossimità dell’apice. Nel presente lavoro è stata simulata la frattura duttile considerando un modello coesivo monodimensionale la cui lunghezza Δ è determinabile ricostruendo una prova sperimentale DWTT mediante l’analisi agli elementi finiti. Le forze di rilascio sui nodi sono calcolate utilizzando la formulazione proposta da Rydholm et. al. [11] che, durante la propagazione tra due nodi successivi, impone la costanza dell’energia specifica dissipata dalla frattura (EFD) a meno di oscillazioni derivanti dalla natura dinamica del fenomeno. Questa condizione deve certamente essere soddisfatta durante la propagazione stazionaria e può essere giudicata accettabile anche in condizioni di propagazione instabile quando, per la dimensione ridotta degli elementi di frattura, si possano trascurare le variazioni di energia specifica dissipata nel passaggio tra due nodi successivi. Supponendo per semplicità che un solo nodo si trovi all’interno della distanza Δ dall’apice di frattura, l’energia specifica EFD può dunque essere espressa in relazione alla forza di rilascio mediante la (1): 0 10 20 30 40 Distance from crack tip [mm] 400 500 600 700 800 900 [M P a] Legend T = 0.8 T = 0.5 T = 0.2 T = 0.0 T = -0.2 T = -0.5 syy PLANE STRESS 0 10 20 30 40 Distance from crack tip [mm] 400 500 600 700 800 900 [M P a] Legend T = 0.8 T = 0.2 T = 0.0 T = -0.2 T = -0.5 syy PLANE STRAIN L http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 40       0 2 dx dx dV F B EFD (1) a) b) Figura 4: a) Schema teorico del modello coesivo; b) definizione del CTOA. Essendo B lo spessore del guscio e F, δ e V rispettivamente la forza di chiusura, la distanza dall’apice e lo spostamento dalla linea di frattura del nodo all’interno della zona di coesione. Il fattore 2 nella formula precedente è introdotto per considerare l’energia dissipata su entrambi i fianchi della cricca. Imponendo la costanza di EFD si ottiene che per ogni istante di tempo deve valere l’uguaglianza seguente:           I B VF EEI B VF f FDFD f 00 2 0 2   (2) dove F0 è la forza agente sul nodo prima della fase di rilascio, Vf lo spostamento finale del nodo quando si trova alla distanza Δ dall’apice mentre per I(α) vale la seguente:             1 1 d du d I (3)     release nodeBegin 1 release node of End 0  (4) La funzione 0 F F   esprime la legge di riduzione della forza sui nodi ed è l’incognita da determinare in relazione alla posizione δ mentre f V u V  è lo spostamento adimensionale. Assumendo una relazione lineare tra la forza agente su un nodo nella zona di coesione ed il relativo spostamento, ne consegue che  1 u   . Sostituendo nelle equazioni (2) e (3) si deduce la funzione di rilascio delle forze da applicare nella zona di coesione: 1 k      (5) Quest’ultima è stata impiegata come legge di forza nello strato coesivo anche in presenza di più nodi in fase di rilascio (Δ maggiore della dimensione degli elementi di frattura) ed utilizzando i dati sperimentali della prova DWTT per determinare il valore di Δ. L’espressione della EFD (1) è in accordo con la definizione di Essential Work Of Fracture proposta da Cotterell [9]. Usando le precedenti relazioni ed introducendo la tensione di rottura del materiale σu (in relazione con il valore F0) si può stimare l’assorbimento locale di energia (6):  CTOAE uDL tan (6) http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 41 Di fatto durante la propagazione stazionaria della frattura esiste un legame univoco tra l’angolo di frattura, la dimensione della zona di coesione e l’energia specifica dissipata per l’avanzamento. La scelta di imporre CTOAC e distanza Δ equivale, ai fini della meccanica della frattura duttile, all’aver fissato il livello di dissipazione energetica sull’apice. Al contempo, potendo assumere C CCOD CTOA   ne resta determinata anche la separazione tra le superfici di frattura. Nel caso comune di più nodi all’interno della zona coesiva, l’evoluzione della forza nodale può essere considerata indipendente per ciascun nodo e pertanto si può continuare ad usare la (1) per ciascuno dei nodi, mentre l’energia specifica di frattura si dovrà calcolare come somma di tutti i contributi:   NODES i i i i FD dx dx dV F B E   0 12 (7) TECNICA DEL RILASCIO CONTROLLATO DEI NODI a simulazione della propagazione della frattura duttile mediante l’utilizzo degli elementi finiti è basato sulla tecnica del rilascio progressivo dei nodi, e consiste nell’implementazione di un layer coesivo posizionato a partire dall’apice virtuale della cricca. All’interno di esso si ha un graduale allontanamento delle superfici di frattura fino all’effettiva separazione, al termine della FPZ. In pratica, come illustrato nella Fig.5a, si possono distinguere un apice virtuale, dove ha inizio il softening ed uno finale, ove si ha il definitivo allontanamento dei lembi della frattura. Tale modello è stato implementato nell’algoritmo di risoluzione di un codice esplicito agli elementi finiti, opportunamente realizzato, in grado di simulare la propagazione della frattura duttile su geometrie generali, con la possibilità di computare anche impatti con corpi rigidi. Secondo l’algoritmo, sui nodi di frattura che si trovano all’interno della zona di processo agiscono delle forze che si oppongono al rilascio dei nodi stessi. Tali forze nodali sono massime in corrispondenza di un nodo sito all’apice virtuale della frattura (vedi Fig.5a) e si annullano quando il nodo esce dalla zona di processo, ovvero ad una distanza dall’apice maggiore o al più uguale alla dimensione della zona coesiva . Conseguentemente, un nodo potrà considerarsi completamente rilasciato solo nel momento in cui si troverà al di fuori di una dimensione caratteristica; essa non dipende in realtà solo dal materiale, ma anche dalle condizioni generali di constraint; quella geometrica più immediata è lo spessore. L’andamento della forza nodale ipotizzato, in funzione della distanza δ dall’apice stesso, segue la legge di distribuzione data dall’equazione (5):   0 1 k F F          (8) avendo indicato con δ la distanza di un nodo in fase di rilascio dall’apice della frattura e con Δ la dimensione della zona di processo ovvero la distanza ove le forze si estinguono. Tale modello, di natura monodimensionale, permette di computare direttamente l’energia dissipata per la formazione di nuove superficie di frattura, attraverso semplice integrazione della (1), in tutto il dominio della FPZ. (a) (b) Figura 5: Schema del modello coesivo monodimensionale (a) e distribuzione delle forze nodali (b). L http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 42 Il valore dell’esponente k della (8), per strutture caratterizzate da spessore sottile e realizzate con materiali a comportamento duttile, viene solitamente assunto pari a 0.5 [11]. Tuttavia, osservazioni sperimentali su provini di consistente spessore hanno evidenziato un considerevole rapporto di riduzione dello spessore in vicinanza dell’apice della cricca, con un effetto importante nel passaggio da forze a tensioni coesive. Questa considerazione ha un effetto notevole sulle sollecitazioni coesive stesse (computate per unità di spessore) che non era stato evidenziato in [11] ove gli spessori erano di fatto molto sottili. Per tenere in considerazione l’influenza di tale necking sulla distribuzione di tensioni coesive nella FPZ, è sufficiente imporre 1.5 come valore di k, che deriva dal considerare una variazione di spessore della cricca dal suo valore massimo al Virtual Tip a quello nullo al Real Crack Tip secondo una modulazione lineare (vedi Fig.5a). FORMULAZIONE BI-PARAMETRICA E TENSIONE DI RIFERIMENTO C a corretta valutazione della tensione di riferimento σC è fondamentale per l’applicazione del modello coesivo ad una struttura discretizzata con gli elementi finiti: ogni forza coesiva applicata ai nodi di frattura è funzione di questa tensione. In aggiunta, il valore dell’Essential Work of Fracture dipende sostanzialmente dall’intensità della σC. Una corretta valutazione del Wf determina il vantaggio di poter utilizzare quest’ultimo come parametro di controllo alternativo per la propagazione della frattura. Il Wf è una caratteristica fortemente locale, quindi si presta ad essere usato in alternativa a quei parametri consolidati della meccanica della frattura (CTOA, COD, etc.), che in determinate circostanze risultano inadatti a valutare le condizioni di propagazione della cricca: un esempio sono i gasdotti interrati allorché si abbia una frattura di elevata velocità che attraversa un Crack-Stopper di contenimento. Per determinare la tensione di riferimento σC, e prendendo spunto dalla teoria HRR [1], è stata definita una funzione analitica a formulazione bi-parametrica che è in grado di caratterizzare compiutamente l’andamento della sola componente di tensione agente in direzione perpendicolare alla cricca: 1 1 a yy x f b xs s               (9) Dove :  σyy : tensione perpendicolare alla cricca  a, b: parametri caratteristici da definire caso per caso  x: distanza nodale dall’apice  s: spessore. Se la (9) risulta valida, in una simulazione agli elementi finiti allorché siano noti i valori di tensione e distanze nodali dall’apice della frattura di almeno una coppia di nodi, è possibile determinare il valore della σyy in un qualsivoglia punto della regione di materiale plasticizzata e prospiciente la cricca. La (9) è evidentemente non lineare nei parametri a e b, ma presenta la peculiarità di essere facilmente linearizzabile mediante il passaggio ai logaritmi, consentendo una valutazione pressoché immediata dei suoi parametri oltre che facilitandone la ricerca ottimale se si dispone di più valori:   1 ln ln 1 a yy b x s              Ponendo la condizione  lnB b :    ln ln 1yy xB a s     (10) Dalla (10) si evince come B rappresenti proprio il valore della tensione in corrispondenza dell’apice, e quindi il valore di riferimento cercato per la taratura del modello coesivo. Il valore di riferimento è dipendente da fattori quali la natura del materiale, le condizioni di carico e i fattori di natura geometrica, i cui molteplici effetti risultano difficilmente disaccoppiabili. L http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 43 Al fine di verificare l’affidabilità della (9) nella caratterizzazione dello stato tensionale in una piastra criccata costituita da un materiale a comportamento elasto-plastico, sono state condotte ad hoc delle simulazioni agli elementi finiti, applicate al modello descritto nel paragrafo precedente, e caratterizzate dalle seguenti configurazioni di carico: Figura 6: Rappresentazione dell’area di campionamento mobile. Caso 1: Presenza di un solo carico distribuito in direzione perpendicolare alla frattura (MODO I). Caso 2: Carico distribuito perpendicolarmente alla frattura e presenza di ulteriore momento flettente. Caso 3: Carico distribuito perpendicolarmente alla frattura e aliquota aggiuntiva di T-stress. Caso 4: Carico distribuito in direzione perpendicolare alla frattura (MODO I) con una distribuzione di forze coesive davanti l’apice della cricca. In Fig.7 sono riportati gli andamenti della yy per le diverse tipologie di carico presentate, relativamente alla regione di materiale posta di fronte l’apice e già interessata da plasticità. La sovrapposizione dei trend ottenuti dalle simulazioni numeriche e dall’applicazione della (9) evidenziano come la funzione bi-parametrica (curve a tratto continuo) sia in grado di interpolare con estrema accuratezza i dati ottenuti dalle simulazioni agli elementi finiti. Figura 7: Confronto tra i dati delle simulazioni FEM e fitting della funzione bi-parametrica. La (9), nella sua forma linearizzata (10) risulta particolarmente indicata per valutare la tensione di riferimento all’interno di un codice agli elementi finiti con valutazione attendibile degli stress agenti. Le elevate variazioni che invece si presentano nella formulazione esplicita richiedono necessariamente l’utilizzo delle medie mobili, nel tempo e nello spazio. Il modo migliore per computarle consiste nell’introduzione di una finestra dinamica di campionamento (vedi Fig. 6). Quest’ultima si muove solidalmente all’apice della frattura, e presenta un’estensione pari a circa 10 volte lo spessore. Per i nodi che ricadono nell’area di campionamento e che appartengono ad elementi già plasticizzati, vengono memorizzati i valori istantanei di stress e delle distanze nodali dall’apice virtuale (Fig. 5b). Per una più rapida ed efficiente valutazione dei parametri, la finestra di campionamento è divisa in due sottoaree: la prima, più significativa in quanto più vicina all’apice della cricca, risulta però meno estesa della seconda. Particolareggiando la (10) con le coppie di valori campionati per ogni nodo, si perviene ad un sistema sovra determinato di equazioni algebrico lineari, la cui risoluzione prevedrebbe il ricorso a tecniche di regressione che risultano sconvenienti in termini computazionali, vista la necessità di operare su una base di numerosissimi time step negli elementi finiti a http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 44 formulazione esplicita. Un approccio alternativo, assai più vantaggioso, consta nel raggruppare le equazioni del sistema in funzione delle sottoaree di campionamento precedentemente definite (Fig. 6). Il risultato è un semplice sistema di due equazioni lineari nelle due incognite a e B.     1 1 1 2 2 2 ln ln 1 ln ln 1 i y SubAi SubA SubA i y SubAi SubA SubA x n B a s x n B a s                            →   1 2 2 1 1 1 SubA SubA SubA SubA SubA n D A n a n C E n B A a C n                (11) con   1 ln y SubA i A   1 ln 1 i SubA x C s        2 ln 1 i SubA x E s          2 ln y SubA i D   (12) APPLICAZIONE DEL MODELLO COESIVO ll’interno del codice di calcolo PiCPro®, solutore agli elementi finiti a formulazione esplicita [12,13,14,15] di proprietà del Centro Sviluppo Material S.p.A. di Roma, è stato implementato il modello coesivo che utilizza le (11) per la determinazione della tensione di riferimento σC. Il PiCPro® è uno strumento molto potente, utilizzato prevalentemente per lo studio e l’analisi della frattura longitudinale di gasdotti interrati. Il codice dispone di features specifiche che consentono di simulare i fenomeni complessi d’interazione che coinvolgono la condotta durante il suo danneggiamento, quali ad esempio il contatto e l’eventuale espulsione del terreno circostante e la decompressione del gas in fuoriuscita dai lembi fratturati del tubo. L’applicazione di un modello coesivo ad una struttura discretizzata con gli elementi finiti prevede che ai nodi di frattura, che risiedono all’interno della zona di processo, siano applicate delle forze vincolari proporzionali alla distanza nodale dall’apice della cricca. È evidente quindi come la presenza di almeno un nodo all’interno della zona di processo sia una condizione necessaria per il funzionamento del modello coesivo, comportando conseguentemente una correlazione vincolante tra estensione della FPZ ed il grado di risoluzione della mesh utilizzata. Una modalità d’implementazione che permette di eludere tale dipendenza, mantenendo inalterata l’efficacia del modello, si ottiene determinando sistemi di forze nodali equivalenti; essi sono ottenuti imponendo la distribuzione di tensioni coesive (Fig. 4a) equivalenti al carico applicato, e garantendo la non violazione delle condizioni di equilibrio e congruenza degli elementi sulla frattura. In questo modo è possibile estendere l’applicabilità del modello coesivo anche a strutture modellate mediante elementi finiti di dimensioni rilevanti od anche, all’opposto, a materiali caratterizzati da scarsa duttilità (Δ molto bassi); il tutto senza ricorrere a mesh refinement troppo spinti a tutto vantaggio dei tempi di simulazione. Di seguito vengono presentati i risultati ottenuti dall’applicazione del codice PiCPro® su una tipologia di gasdotto interrato le cui caratteristiche sono riportare in Tab. 1: Materiale E [MPa] Rs [MPa] Rm [MPa] Δ [mm] Lunghezza [m] Raggio [inch] Spessore [mm] Gas Pressione [bar] X100 206000 775 833 35 12 48 18.4 metano 187 Tabella 1: Caratteristiche del gasdotto Per valutare l’affidabilità del modello e testarne le proprietà di mesh-indipendenza, sono state preparate due diverse discretizzazioni agli elementi finiti della condotta (Fig. 8). La mesh di Fig. 8a è composta da 10320 elementi e 10951 nodi, e presenta un infittimento progressivo delle maglie in direzione circonferenziale, dalla generatrice inferiore del tubo a quella superiore. In Fig. 8b è rappresentata la seconda tipologia di mesh adottata, composta da 17984 elementi e 18395 nodi, in cui è visibile un graduale e continuo mesh rifinement in direzione contraria all’asse z, caratterizzato dall’assenza di elementi per la transizione tra regioni della mesh a diverso infittimento. La generatrice superiore del tubo è nella mesh TIPO 1 costituita da elementi di egual lunghezza, pari A http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 45 a 11 mm; viceversa, nella seconda tipologia si osservano elementi ad estensione variabile progressivamente dai 24 mm ai 46 mm (per leggibilità grafica solo una parte assiale del tubo è evidenziate in Fig. 8). (a) (b) Figura 8: Mesh ad infittemento circonferenziale (a); Mesh ad infittemento assiale progressivo (b). La mesh di Fig. 8a è composta da 10320 elementi e 10951 nodi, e presenta un infittimento progressivo delle maglie in direzione circonferenziale, dalla generatrice inferiore del tubo a quella superiore. In Fig. 8b è rappresentata la seconda tipologia di mesh adottata, composta da 17984 elementi e 18395 nodi, in cui è visibile un graduale e continuo mesh rifinement in direzione contraria all’asse z, caratterizzato dall’assenza di elementi per la transizione tra regioni della mesh a diverso infittimento. La generatrice superiore del tubo è nella mesh TIPO 1 costituita da elementi di egual lunghezza, pari a 11 mm; viceversa, nella seconda tipologia si osservano elementi ad estensione variabile progressivamente dai 24 mm ai 46 mm (per leggibilità grafica solo una parte assiale del tubo è evidenziate in Fig. 8). (a) (b) Figura 9: Andamento del CTOA (a) ed Essential Work of Fracture (b) di una frattura viaggiante a velocità costante (230000 mm/s) su un modello con mesh di TIPO 1. (a) (b) Figura 10: Andamento del CTOA (a) ed Essential Work of Fracture (b) di una frattura propagante a velocità costante (230000 mm/s) su un modello con mesh di TIPO 2. http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 46 Nelle Figg. 9 e 10 sono riportati gli andamenti del CTOA e dell’Essential Work of Fracture relative a simulazioni condotte con le due tipologie di mesh (Figg.8a e 8b) avendo imposto alla cricca una velocità di propagazione costante pari a 230000 mm/s. Il confronto diretto tra i trend di Figg. 9a e 10a evidenzia come la variazione di assetto dei lembi di frattura durante la propagazione si mantenga pressoché circoscritta lungo tutta l’estensione della condotta, in un range compreso tra i 15° ed i 20°, al netto delle perturbazioni di natura dinamica, amplificate anche dalla particolare modalità di propagazione utilizzata. (a) (b) Figura 11: Variazione della velocità della cricca (a) ed Essential Work of Fracture (b) di una frattura in controllo di CTOA (CTOAc = 16°) su un modello con mesh di TIPO 1. (a) (b) Figura 12: Variazione della velocità della cricca (a) ed Essential Work of Fracture (b) di una frattura in controllo di CTOA (CTOAc = 16°) su un modello con mesh di TIPO 2. I picchi presenti negli andamenti del CTOA e dell’Essential Work of Fracture sono causati dalla trasmissione di onde elastiche a bassa frequenza, che “rimbalzano” da un estremo all’altro del tubo, determinando una più o meno accentuata modificazione del profilo di frattura (vedi Fig. 13). A tal proposito è interessante osservare come la diversa configurazione delle due mesh determini una sostanziale differenza della trasmissibilità delle perturbazioni dinamiche. In particolare la mesh di TIPO 1, essendo caratterizzata da una dimensione media degli elementi decisamente superiore rispetto alla mesh di TIPO 2, riesce a filtrare le onde elastiche a più alta frequenza, che invece sono apprezzabili nel modello ad infittimento assiale progressivo (si notino in Fig. 10a i tre picchi in sequenza di magnitudo 30°). I trend relativi all’Essential Work of Fracture (Figg. 9b e 10b) mostrano una variabilità accettabile, analoga a quella evidenziata dagli andamenti del CTOA, presentando i caratteristici picchi determinati dalla propagazione delle onde elastiche. Questi ultimi, comunque, presentano valori proporzionalmente meno elevati rispetto a quelli di Figg. 9a e 9b, per effetto dell’introduzione della media mobile utilizzata non solo per la computazione delle forze coesive ma anche dell’Essential Work of Fracture. Gli andamenti di Figg. 11 e 12 sono stati ottenuti utilizzando una modalità di calcolo diversa che non prevede l’imposizione della velocità di avanzamento, ma lo studio dell’evoluzione della propagazione governata dai parametri caratteristici della meccanica della frattura, nel caso specifico il CTOA (CTOAc). Tale modalità è definita a “propagazione libera”, ovvero l’avanzamento della frattura è dipendente dal confronto tra il valore attuale del parametro descrittivo scelto ed il valore critico del medesimo, che è una proprietà intrinseca del materiale del tubo: la velocità con cui la cricca viaggia è ottenuta mediando i valori relativi agli atti elementari di avanzamento/arresto. Il valore di 16° per il CTOAc è stato http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 47 determinato dall’osservazione dei grafici di Figg. 9a e 10a, in modo tale da poter prevedere una velocità di propagazione vicina a quella imposta (230000 mm/s) nelle precedenti simulazioni. (a) (b) Figura 13: Effetto delle onde elastiche sul profilo di frattura per la mesh TIPO 1(a) e TIPO 2(b). Figura 14: Andamento della tensione di Von Mises [MPa ] nel tubo in due step di calcolo successivi. L’osservazione delle Figg. 11a e 12a mostra effettivamente una velocità di avanzamento che si attesta su un valore medio circa pari a 230000 mm/s, con una variabilità dei trend assai contenuta e visibilmente indipendente dalla particolare discretizzazione utilizzata (soprattutto se confrontati con gli andamenti di Figg. 9a e 10a). I trend rappresentativi dell’Essential Work of Fracture nelle Figg. 11b e 12b mostrano come l’energia specifica per la formazione di nuove superfici di frattura sia caratterizzata da una sorprendente stabilità in valore, determinata molto probabilmente dall’assetto della frattura condizionato dalla modalità a “propagazione libera” con controllo del CTOA. Questo risultato è paritcolarmente importante in virtù della considerazione che nel secondo modello la dimensione degli elementi è addirittura inferiore a quella della distanza di estinzione. Figura 15: Andamento della deformazione totale nel tubo in due step di calcolo successivi. Variazione profilo Variazione profilo http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it M. Minotti et alii, Frattura ed Integrità Strutturale, 11 (2009) 36-48; DOI: 10.3221/IGF-ESIS.11.04 48 Tutto ciò suggerisce quindi come l’Essential Work of Fracture possa essere utilizzato in alternativa al CTOA come parametro valutativo della propagazione di una cricca in tutte quelle circostanze in cui la configurazione geometrica posseduta dai lembi di rottura risulta inadatta od insufficiente per la valutazione delle effettive condizioni di propagazione. La necessità di effettuare medie spazio-temporali rispetto alla posizione attuale della frattura è infine evidenziata con chiarezza dalle Figg. 14 e 15. Esse si riferiscono a due step temporali di poco successivi per i quali si può assumere che la cricca abbia subito un avanzamento di entità molto modesta. Il confronto sinistro-destro mostra che le tensioni equivalenti variano con grandissima velocità, essendo in gran parte dominate da sovrapposizione di onde elastiche di varia natura. Tuttavia la ricerca di valori medi ha senz’altro un significato importante in quanto le deformazioni totali, meno sensibili agli effetti d’onda, risultano molto più costanti nel tempo come visibile in Fig. 15. Quindi il fenomeno, pur se mascherato da importanti effetti dinamico-elastici che rendono impossibile la stima tensionale tout court, consente purtuttavia di estrapolare una stato di tensione persistente. È necessario considerare a riferimento spaziale non la posizione fisica rispetto al tubo, ma quella relativa all’apice della frattura. CONCLUSIONI i è proposta una formulazione analitica bi-parametrica per la determinazione del valore di tensione, agente in direzione ortogonale al piano di frattura, di una struttura criccata sollecitata prevalentemente in MODO I e costituita da un materiale a comportamento elasto-plastico. Il valore di tensione calcolato in concomitanza dell’apice è utilizzato come valore di riferimento per un modello coesivo (utilizzato per la simulazione degli effetti di softening correlati alla propagazione della frattura) implementato all’interno di un codice agli elementi finiti a formulazione esplicita per l’analisi della frattura dei gasdotti interrati. Sono stati inoltre proposti i risultati ottenuti dall’applicazione del suddetto codice per differenti tipologie di mesh e modalità di propagazione su una condotta di 12 m esercita con metano a 187 bar. I trend relativi alla variazione di velocità di propagazione e di assetto (CTOA) della frattura hanno mostrato un’ottima efficienza del modello coesivo nella simulazione del danneggiamento del materiale prospiciente l’apice della cricca; in particolare si è potuta apprezzare una notevole stabilità nella computazione dell’Essential Work of Fracture, consentendo di adoperare quest’ultimo come un ulteriore parametro di valutazione per la propagazione della frattura. BIBLIOGRAFIA [1] T.L. Anderson, Fracture Mechanics, Fundamentals and Applications, CRC Press, (1995). [2] G. Berardo, P. Salvini, G. Mannucci, G. Demofonti, in: International Pipeline Conference, Calgary, Canada, (2000). [3] W.A. Maxey, in: 5th Symposium on Line Pipe Research, AGA, (1974). [4] M.F. Kanninen, C.H. Popelar, Advanced Fracture Mechanics, Oxford University Press, (1985). [5] G. Berardo, P.Salvini, M. Di Biagio, in: XIII European Conference on Fracture, Elsevier, San Sebastian, (2000). [6] J. Lemaitre, J.L. Chaboche, Mechanics of Materials, Cambridge University Press, Cambridge (1994). [7] N. Bonora, P. Salvini, A. Fonzo, in: Atti XXX Convegno AIAS, Alghero (SS), (2001). [8] G. Demofonti, M. Minotti, P. Salvini, Frattura ed Integrità Strutturale, 5 (2008) 39. [9] B. Cotterel, J.K. Reddel, Int. Journal of Fracture, 13 (1977) 267. [10] C. Betegon, J.W. Hancock, Journal of Applied Mechanics, 58 (1991) 104. [11] G. Rydholm, B. Fredriksson, F. Nilsson, in: Proc. of 1st Int. Conf. On Numerical methods in fracture mechanics, Swansea, (1978). [12] J.R. Rice and G.F. Rosengren, J. of the Mechanics and Physics of Solids, 16 (1968) 1. [13] C.E. Turner, O. Kolednik, Fatigue & Fracture of Engineering Materials & Structures, 17 (1994) 1089. [14] J.O. Hallquist, LS-DYNA3D Theoretical Manual, Livermore Software Technology Corporation, (1991). [15] T. Belytschko, J.I. Lin, C.S. Tsay, Computer methods in applied mechanics and engineering, 42 (1984) 225. S http://dx.medra.org/10.3221/IGF-ESIS.11.04&auth=true http://www.gruppofrattura.it