Numeričko modeliranje oštećenja u heterogenim materijalima primjenom teorije faznih polja

Material failure is still one of the central issues in modern engineering. Its prediction and prevention are being tackled early in the product design stage. Materials used in modern engineering often exhibit significant microstructural heterogeneity. The size, shape, distribution and properties of microconstituents considerably influence the heterogeneous material properties. Nowadays, numerical simulations play an essential role in component design and material development, gradually supplanting and replacing expensive and time-consuming experiments. However, it is worth noting that many complex fracture processes occur at microstructural scale making the fracture analysis an especially challenging and interesting problem. One of the methods for the numerical modelling of fracture capable of efficiently recovering these complex fracture processes is the recently emerged phase-field method. It approximates the sharp crack discontinuity with a diffusive band regulated by a length-scale parameter, thus separating the broken and intact material states. Although extensive research has been carried out on the development of phase-field fracture theory over the past decade, certain challenges still exist in the computational implementation of the method. Within the finite element framework, a fine spatial discretization is often required to resolve the smooth phase-field distribution regulated by a small length scale parameter. Thus, coupling phase-field method with an inefficient solution scheme can be computationally rather expensive. In this work, a novel generalized phase-field framework capable of simultaneously recovering brittle, ductile, and fatigue fracture in three-dimensional settings is developed. The robustness and accuracy of the results is ensured by the development of an efficient residual control algorithm and its implementation in the commercial finite element software ABAQUS. Major advantage of such implementation is the high usability of different underlying solvers, convergence criteria and other additional options including automatic incrementation, element deletion, coupled contact analysis, thread parallelization or restart analyses, used and thoroughly discussed in this work. The full source code together with the presented examples and explanations is made publicly available, thus promoting the phase-field fracture methodology. The proposed implementation is exhaustively tested on a large number of different benchmark examples, verified and validated in comparison with numerical and experimental data from available literature. The importance of a stopping criterion within a staggered scheme is emphasised in the illustrative examples. The detailed discussions regarding the proposed implementation’s accuracy and CPU time usage are provided. Special attention is given to the verification of fatigue fracture examples through the parametric study. Main features of fatigue, including Wöhler and Paris law curves in low- and high-cycle regimes, are easily recovered without any additional criteria. The cycle skipping technique is implemented to allow for the calculation of a very high number of cycles on moderate size examples. Finally, the potential of the implementation is demonstrated on the different-sized samples of nodular cast iron with their actual heterogeneous microstructure obtained from microtomography. Specimens are tested in both monotonic and cyclic loading regimes. The results clearly show the size-effect behaviour as well as the influence of microstructural topology on fracture patterns.UVOD Neke od glavnih karakteristika modernih konstrukcija i strojnih komponenti su njihova sve veća složenost te učestalo korištenje novo razvijenih materijala s poboljšanim svojstvima. Time se zadovoljavaju potrebe za iznimnom sigurnošću, pouzdanošću i trajnošću, istovremeno zadovoljavajući visoke zahtjeve učinkovitosti i niskih cijena izrade, odražavanja i upravljanja. Povećanje životnog vijeka komponenti također je važan čimbenik u rješavanju gorućih pitanja kao što su klimatske promjene i očuvanje okoliša. Degradacija i lom materijala još su uvijek problemi koji često dovode do materijalnih i financijskih gubitaka uzrokovanih oštećenjima proizvoda i neisporučenim uslugama, ili u ekstremnijim slučajevima, dovode do gubitaka života. Prema tome, predviđanje i prevencija pojave oštećenja i loma materijala u ranoj fazi konstruiranja nedvojbeno su iznimno važni, i kao takvi, od velikog su interesa mnogim inženjerima i istraživačima. Numeričke simulacije nezaobilazan su dio konstruiranja proizvoda i razvoja novih materijala. Nakon odgovarajuće validacije, numeričke simulacije nadopunjuju te ponekad i u potpunosti zamjenjuju često skupe i dugotrajne eksperimente, smanjujući broj potrebnih prototipa i olakšavajući plansko održavanje. Većina materijala koji se koriste u modernim konstrukcijama i strojnim komponentama sadrži određenu razinu mikrostrukturne heterogenosti. Razni geomaterijali poput betona i stijena, čelici visoke čvrstoće, kompoziti i polimeri, ili materijali dobiveni postupcima poput sinteriranja ili aditivne proizvodnje (3D printanja) samo su neki od primjera takvih mikrostrukturno heterogenih materijala. Na njihova materijalna stvojstva u velikoj mjeri utječu pojedinačna svojstva mikrokonstituenata te njihova veličina, oblik i raspodjela. Upravo se na mikrostrukturnoj razini pojavljuje većina složenih procesa loma kao što su lokalizacija oštećenja, nastanak, propagacija, grananje i srastanje pukotina, što analizu loma heterogenih mikrostruktura čini izrazito izazovnom i zanimljivom. Stoga su numeričke analize oštećenja i loma heterogenih materijala koje u obzir uzimaju mikrostrukturnu topologiju i svojstva mikrokonstituenata važan čimbenik u procjeni pouzdanosti i strukturalne cjelovitosti, kao i predviđanju radnog vijeka modernih konstrukcija i strojnih komponenata. Također bi mogle dovesti do realnijeg opisa ponašanja materijala čime bi se olakšao razvoj novih materijala s poboljšanim svojstvima koja proizlaze iz manipulacije mikrostrukturom. Time bi se omogućila proizvodnja učinkovitijih, sigurnijih i jeftinijih konstrukcija izrađenih od heterogenih materijala. MODELIRANJE OŠTEĆENJA I LOMA Numeričke metode modeliranja oštećenja i loma najčešće su razvijene u sklopu metode konačnih elemenata (MKE). MKE se temelji na konceptu podjele fizikalnog modela na manje segmente (elemente) jednostavne geometrije i konačnog broja stupnjeva slobode. Takav se diskretni model, koji se sastoji od mreže konačnih elemenata, može opisati sustavom algebarskih jednadžbi i jednostavno riješiti. Numeričke metode modeliranja oštećenja i loma općenito se mogu podijeliti na diskretne i difuzne (tj. kontinuumske) prema načinu njihovog opisa oštećenja, odnosno loma. Diskretne metode pukotinu opisuju kao oštar geometrijski diskontinuitet u polju pomaka. Najpoznatije i najčešće korištene diskretne metode su metoda modeliranja kohezivnih zona (eng. cohezive zone modelling, CZM), razni algoritmi lokalne promjene mreže konačnih elemenata (eng. remeshing) te tehnike obogaćivanja konačnih elemenata čiji je najpoznatiji predstavnik proširena metoda konačnih elemenata, XFEM (eng. extended finite element method). Iako su se diskretne metode do sada pokazale točnim i vrlo korisnim u predviđanju oštećenja i loma, uz neizbježno korištenje dodatnih kriterija nastajanja i rasta pukotina, njihov najveći problem zapravo je numeričko praćenje diskontuiteta, odnosno pukotine. Tako kod pristupa koji se temelje na metodi konačnih elemenata pukotina raste samo duž rubova elementa što uzrokuje problem ovisnosti rezultata o veličini mreže (eng. mesh dependency) i problem ovisnosti rasta pukotine o usmjerenosti konačnih elemenata (eng. bias dependency). Taj problem najviše dolazi do izražaja pri složenim topologijama loma kakve se očekuju na miktrostrukturnim geometrijama ili složenim trodimenzijskim problemima. Posljednjih godina razvijene su i brojne bezmrežne metode koje omogućuju modeliranje razvoja pukotina bez remeshinga, kao npr. bezmrežna Galerkinova metoda (eng. Element-free Galerkin, EFG). Iako se često koriste za modeliranje linearnih i nelinearnih problema rasta pukotina, u usporedbi s MKE, takve metode imaju nekoliko još neriješenih problema među kojima se najviše ističe numerička nestabilnost i dugo vrijeme računanja. S druge strane, kontinuumske metode modeliranja oštećenja i loma, umjesto modeliranja diskretne pukotine i stvaranja novih površina, uvode parametar oštećenja na razini integracijske točke. Njime se kontrolira krutost materijala što omogućuje smanjenje naprezanja uzrokovano stvaranjem pukotine. U modelima utemeljenim na pretpostavci standardnog ili tzv. lokalnog kontinuuma pojavljuje se lokalizacija oštećenja koja uzrokuje lokalni gubitak eliptičnosti sustava diferencijalnih jednadžbi, zbog čega numerička rješenja često ne konvergiraju prema fizikalno smislenom rješenju. Kao što je ranije spomenuto, to se u okviru MKE-a očituje rezultatima ovisnim o gustoći i usmjerenosti mreže konačnih elemenata. Navedeni problemi riješeni su uvođenjem nelokalnih i gradijentno-poboljšanih kontinuumskih pristupa. Zajedničko svojstvo takvih modela je dodavanje parametra duljinske skale čime se omogućuje i računanje efekta veličine uzorka (eng. size-effect), što nije moguće s modelima temeljenim na pretpostavkama lokalne mehanike kontinuuma. Poznati nedostatci takvih, sad već standardnih, ne-lokalnih pristupa jesu rast zone oštećenja okomito na smjer rasta pukotine. Iako se taj nedostatak uspješno rješava uvođenjem parametra duljinske skale kao funkcije naprezanja, fizikalna utemeljenost takve pretpostavke je upitna. Drugi način rješavanja navedenog problema je primjena kontinuuma višeg reda (gradijentne, nelokalne i mikrokontinuumske teorije), međutim, to neizbježno uvodi nove komplikacije tokom numeričke ugradnje navedenih modela. MODELIRANJE HETEROGENIH MATERIJALA Iako je pretpostavka homogenosti materijala na materijalnoj točki makrorazine još uvijek valjana u mnogim primjerima iz inženjerske prakse, razvoj naprednih materijala korištenih u modernim konstrukcijama zahtijeva razmatranje materijalne mikroheterogenosti i njezinog utjecaja na konstitutivno ponašanje materijala na makrorazini. Eksperimentalne tehnike također su znatno napredovale u procjeni mikrostrukturne heterogenosti. Najbolji primjer je 3D rendgenska mikrotomografija koja je sada u stanju rutinski generirati realistične geometrijske modele mikrostruktura mnogih materijala na raznim skalama. Jedan od načina numeričkog modeliranja heterogenosti upotreba je izravnih numeričkih simulacija (eng. direct numerical simulation, DNS) pri čemu se heterogena mikrostruktura detaljno modelira izravno na makrorazini. To je ujedno i najpreciznija metoda za numeričku predikciju loma u heterogenim materijalima. Međutim, računalno je vrlo zahtjevna, posebno kad se radi o velikim razlikama između skala što zahtjeva vrlo gustu diskretizaciju konačnim elementima. S druge strane, razvijene su i višerazinske metode koje na odgovarajući način razdvajaju i prenose rješenja između različitih skala. Najpopularnija takva metoda zasigurno je računalna homogenizacija (eng. computational homogenization) koja uprosječuje rezultate dobivene analizom problema na mikrostrukturnoj razini, te ih prosljeđuje numeričkom modelu na makrorazini. Takva analiza provodi se na tzv. reprezentativnim volumenskim elementima (RVE). Računalna homogenizacija standardno se koristi u problemima elastičnog i elastoplastičnog ponašanja materijala te spregnutih termomehaničkih problema. Također, razvijene su i metode računalne homogenizacije drugog reda koje u obzir uzimaju i gradijente deformacije na makrorazini koristeći prethodno spomenute modele kontinuuma višeg reda. Međutim, kad je u pitanju modeliranje oštećenja i loma, još uvijek postoje veliki problemi i otvorena pitanja u računalnim metodama homogenizacije. Možda i najvažnije, odnosi se na upitno postojanje RVE-a prilikom pojave oštećenja i loma. Naime, po svojoj definiciji, RVE gubi svoju statističku reprezentativnost prilikom lokalizacije oštećenja i pojave loma. Iako alternativne homogenizacijske sheme ukazuju na mogućnost rješavanja ovog problema za kvazi-krhke materijale, čini se da je dosljedna i pouzdana primjena višerazinskih metoda predikcije oštećenja i loma na realnim konstrukcijama još uvijek daleko. Prema tome, razvoj pouzdanih modela predikcije oštećenja i loma na vrlo složenim geometrijama kakve se javljaju na mikrorazini heterogenih materijala područje je od velikog interesa. Razvoj takvih modela uz razmatranje stvarne topologije mikrostrukture nudi način za preciznije modeliranje složenih procesa loma u heterogenim materijalima te je detaljnije obrađen u ovom radu. TEORIJA FAZNOG POLJA Teorija faznog polja (eng. phase-field modelling) postala je iznimno popularna u posljednjem desetljeću. Primjena metode u području modeliranja loma samo je jedna od mogućnosti ove teorije. Općenito, ona se uspješno koristi za modeliranje sustava s oštrim prijelazima između faza, kao što je npr. interakcija između fluida i čvrstih tijela (eng. fluid-solid interaction). U teoriji faznog polja, uvodi se kontinuirana varijabla čijim se glatkim prijelazom razlikuju fizičke faze unutar određenog sustava. U slučaju loma, ova varijabla razdvaja oštećeni i neoštećeni dio materijala aproksimirajući, na taj način nagli diskontinuitet pukotine difuznim pojasom. Širina takvog pojasa regulirana je parametrom duljinske skale. Prema tome, teorija faznog polja može se svrstati u kontinuumske metode modeliranja oštećenja i loma. Zanimljivo je da je metoda neovisno razvijena od strane dvije znanstvene zajednice, fizičara i inženjera, čiji se pristupi i polazišne točke za dobivanje osnovnih jednadžbi razvoja loma pomoću faznog polja znatno razlikuju. S jedne strane, zajednica fizičara razvila je modele predviđanja dinamičkih lomova koristeći Ginzburg-Landau teoriju, izvorno izvedenu za elektromagnetske fenomene faznih prijelaza. S druge strane, modeliranje loma metodom faznog polja na kojoj se temelji ovaj rad, potječe od varijacijskog pristupa krhkom lomu gdje je osnovna Griffithova teorija loma pretvorena u problem minimizacije energije. Takav princip nalaže da pukotina može nastati ili rasti u čvrstom tijelu ukoliko tako rezultirajuća konfiguracija ima nižu ukupnu energiju u usporedbi s ostalim konfiguracijama u kojima nema nastajanja ili rasta pukotine. Smjerovi rasta pukotine prirodno se definiraju kao oni koji vode u konfiguraciju s minimalnom ukupnom energijom. Pukotine se mogu granati ili spojiti ukoliko to dovodi do konfiguracije s nižom energijom od jednostavnog rasta postojećih pukotina, bez potrebe za uvođenjem dodatnih kriterija rasta oštećenja i loma. Regularizacijom tako izvedenog energijskog funkcionala dobiva se sustav parcijalnih diferencijalnih jednadžbi koje u potpunosti određuju nastanak i rast pukotine. Upravo su prethodno navedene činjenice razlog popularnosti ove metode, koja se pokazala izrazito uspješnom u rješavanju složenih procesa loma, uključujući nastanak, rast, grananje i srastanje pukotina. Također, aproksimacijom diskretne topologije pukotine na fiksnoj mreži konačnih elemenata zaobilazi se složeni problem praćenja površina pukotina i značajno se pojednostavljuje numerička implementacija, posebno u trodimenzijskim postavkama. U zadnje vrijeme razvijen je velik broj takvih modela s različitom primjenom kao što su modeli faznog polja za predviđanje krhkog i duktilnog loma ili multifizikalni problemi kao što su lom uzrokovan termomehaničkim, elektromehaničkim ili hidrauličkim djelovanjem, što također pokazuje veliki potencijal ove metode. Međutim, velika većina ovih modela uglavnom je verificirana kroz usporedbe sa standardnim numeričkim primjerima (benchmark) i kvalitativnom usporedbom s eksperimentalno dobivenom slikom putanje loma. Manji je broj modela kvantitativno validiran u usporedbi s eksperimentom, što je preduvjet kako bi modeli predikcije loma primjenom faznog polja uistinu mogli postati pouzdana alternativa numeričkih simulacija rasta pukotina u inženjerskoj praksi. Dio modela testiran je i na heterogenim mikrostrukturnim geometrijama. Neki od primjera su radovi na geometriji dobivenoj mikrotomografijom betona ili nodularnog lijeva u jednom od autorovih radova, proizvoljnoj (randomiziranoj) mikrostrukturi grafita H-451 ili uranovog dioksida (UO2) gdje su svojstva materijala kalibrirana prema simulacijama molekularne dinamike. Tek je nedavno formulacija faznog polja proširena i na probleme zamornog loma. Predstavljeni modeli jasno pokazuju potencijal metode faznog polja za repliciranje glavnih značajki zamornog loma, uključujući nastanak pukotine, stabilne i nestabilne faze rasta, rezultat čega su poznata Wöhlerova krivulja i Parisov zakon. Navedeni modeli razvijeni su pod pretpostavkom elastičnog ponašanja materijala što odgovara takozvanom visokocikličkom zamornom režimu. Budući da je zamorni lom jedan od glavnih uzroka loma inženjerskih konstrukcija, a njegovo numeričko predviđanje i dalje predstavlja veliki izazov, za očekivati je da će se uskoro pojaviti i još napredniji modeli zamora s temeljitom eksperimentalnom validacijom. Međutim, u ovom trenutku još uvijek nedostaje generalni model koji je sposoban replicirati lom pri monotonom opterećenju, kao i značajke nisko- i visokcikličkog zamora. Ovaj rad se također bavi i s razvojem takvog modela. Određeni problemi u numeričkoj implementaciji modela unutar okvira metode konačnih elemenata i dalje postoje te su predmet intenzivnog istraživanja u znanstvenoj zajednici. Neki od glavnih problema su nekonveksnost temeljnog funkcionala slobodne energije s obzirom na stupnjeve slobode polja pomaka i faznog polja, način modeliranja ireverzibilnosti loma te dugo vrijeme računanja zbog često potrebnih, vrlo gustih mreža konačnih elemenata. Konačno, većina implementacija modela loma metodom faznog polja razvijena je u nekomercijalnim softverima otvorenog koda baziranim na MKE ili kao samostalni algoritmi. Na taj je način značajno ograničen potencijal primjene ove metode u praktičnim inženjerskim problemima izvan znanstvene zajednice. S druge strane, primjena modela loma temeljenih na metodi faznog polja unutar najčešće korištenih komercijalnih računalnih MKE paketa ključna je za promicanje metodologije kod ostalih inženjera, istraživača i studenata. Zbog toga su modeli razvijeni u ovom radu implementirani u komercijalni računalni MKE paket ABAQUS te objavljeni i slobodno dostupni online. HIPOTEZE, CILJEVI, ZNANSTVENI DOPRINOS I ZAKLJUČAK Cilj ovog rada je razviti i implementirati novu unificiranu i robusnu metodu za numeričko modeliranje krhkih, duktilnih i zamornih procesa oštećenja i loma, ovisno o materijalnim svojstvima i vrsti opterećenja. Metoda bi trebala biti općenita te omogućiti modeliranje složenih procesa loma uključujući iniciranja oštećenja, propagaciju, grananje i srastanje pukotina bez uvođenja dodatnih specijalnih uvjeta. Takva metoda bi trebala omogućiti učinkovito rješavanje složenih procesa loma prisutnih na mikrostrukturnoj razini heterogenih materijala. Hipoteza ovog rada je mogućnost razvoja takve numeričke metode primjenom teorije faznog polja za modeliranje oštećenja i loma koja pokazuje prednosti u odnosu na postojeće postupke, kod kojih se pukotina modelira kao geometrijski diskontinuitet. Dobiveni rezultati trebali bi jasno prikazati utjecaj mikrostrukturne topologije te utjecaj veličine uzoraka čime bi numeričke metode mogle pridonijeti boljoj procjeni strukturne pouzdanosti i sigurnosti te razvoju materijala s poboljšanim mehaničkim svojstvima. U sklopu ovog rada razvijena je poopćena numerička metoda rješavanja krhkog, duktilnog te zamornog oštećenja i loma temeljena na teoriji faznog polja. Metoda je implementirana u komercijalni MKE softver ABAQUS. Razvijen je „staggered“ algoritam s novim kriterijem konvergencije, temeljenom na provjeri norme reziduala, čime se poboljšavaja stopa konvergencije. Implementacija navedenog algoritma prva je takva implementacija iterativnog algoritma teorije faznog polja za rješavanje loma u paket ABAQUS. Time je također riješen problem ne-konveksnosti osnovnog funkcionala slobodne energije, koji često dovodi do numeričkih nestabilnosti. Sam algoritam je temeljito testiran čime je potvrđeno da je robustan i učinkovitiji od uobičajenih algoritama loma primjenom teorije faznog polja. Glavne prednosti implementacije unutar komercijalnog softvera su učinkovite iskorištene, uključujući računalnu paralelizaciju i opciju brisanja elemenata kako bi se smanjilo vrijeme računanja; korištenje različitih rješavača koji omogućuju računanje primjera s opterećenjem zadanim putem pomaka, ali i sile; korištenje naprednih kriterija konvergencije; korištenje automatske inkrementacije, opcije ponovnog pokretanja, odnosno nastavka analiza; opcije računanja problema konktakta. Razmatrane su različite 1D, 2D i 3D formulacije elemenata. Kompletan izvorni kod zajedno s primjerima predstavljenim u ovom radu, uputama i objašnjenjima javno je dostupan drugim istraživačima, studentima i inženjerima, promovirajući metodologiju loma primjenom teorije faznog polja. Kod je otvoreno dostupan na: https://data.mendeley.com/datasets/p77tsyrbx2/4 Kroz veliki broj primjera dokazano je da razvijena poopćena metoda može točno reproducirati krhki/duktilni ili zamorni lom, u skladu s ponašanjem materijala i uvjetima opterećenja. U skladu s tim, prikazana je i detaljna rasprava o točnosti rezultata i trajanju analize. Metoda je temeljito validirana i verificirana na primjerima krhkog i duktilnog ponašanja materijala, pod pretpostavkom materijalne homogenosti, usporedbom s numeričkim i eksperimentalnim rezultatima iz dostupne literature.