Modelování vícesložkové difúzní fázové transformace v pevných látkách

1. Programy a Algoritmy Numerické Matematiky 14 (2008) Dolní Maxov Modelování vícesložkové difúzní fázové transformace v pevných látkách Jiří Vala Brno University of Technology, Faculty of Civil Engineering, Czech Republic Vala.J@fce.vutbr.cz tradiční počasí na PANM

2. abstrakt přednášky pro PANM 14 Dolní Maxov 2008

3. Ignis Brunensis FAST VUT v Brně

4. KLÍČOVÁ SLOVA • difúzní a masivní fázová transformace • Onsagerův termodynamický princip • matematický model: stacionární a nestacionární případ • existence řešení a konvergence numerických algoritmů • praktické výsledky numerických simulací Ni-Al-Cr-Ta alloy superaustenitic iron NICROFER Ni superalloy CMSX4 Bi-Sn-Zn alloy

5. OBSAH • Historické poznámky • Difúzní a masivní fázová transformace – fyzikální podstata • Matematický model – stacionární a nestacionární případ • Onsagerův princip maximální rychlosti disipace • Numerické přístupy • Výpočty s podporou softwaru MATLAB a MAPLE • Existenční a konvergenční analýza • Příklad výpočtové simulace • Užitečná zobecnění

6. Historické poznámky: Onsager (1931) – difúzní tok pro každou složku ve vícesložkové soustavě = matice kinetických součinitelů x gradient chemických potenciálů jednotlivých složek, formulace termodynamického extremálního principu Callen (1960), Groot & Mazur (1962), …, Glicksmann (2000) – zobecnění původního Onsagerova přístupu: difúzní tok → “extensive quantity”, chemický potenciál → “intensive variable” Fan & al. (1999) – přenos rozpuštěných složek (“solute drag”)v polykrystalech, binární slitiny, substitučnísložky Schneider & Inden (2004) – ostré rozhraní, 3 typy okrajových podmínek: i) lokální podmínky rovnováhy s oddělením složek (“with partitioning”), ii) lokální podmínkyrovnováhy se zanedbatelným oddělováním (“with negligible partitioning”), iii) para-rovnovážnépodmínky, pohyblivá xpevná referenční soustava Svoboda & al. (20061) – rozhraní konečné tloušťky (ostré rozhraní jako limitní případ), bez umělých okrajových podmínek (uzavřená soustava), substituční slitiny, numerická analýza pro stacionární případ Svoboda & al. (20062) – substituční i intersticiální složky, difúze a creep s neideálními zdroji a místy zániku (“sources& sinks”)vakancí …

7. Difúzní a masivní fázová transformace – fyzikální podstata: trojí podstata difúzních procesů ve vícesložkových systémech: a) vakanční mechanismus pro “pomalu” difundující substituční složky b) “rychlý” pohyb atomů intersticiálních složek c) existence jistých zdrojů a míst zániku vakancí základní předpoklady: • materiály obsahující r + s složek • dvě fáze α a γ, oddělené rozhraním konečné tloušťky h, jehož vlastnosti odpovídají ideální kapalině β • Onsagerův thermodynamický extremální princip • evoluce molárních frakcí c1(x,t),…, cr(x,t) a rychlosti pohybu mezifázového rozhraní v(t) • jednorozměrný sdružený problém objemové difúze a migrace rozhraní

8. xL(t)xR(t) pohybující se rozhraní v průběhu fázové transformace

9. Matematický model – stacionární a nestacionární případ: molární frakce (podíl) c1(x,t) +…+ cr(x,t) = 1† c* = (c1, …, cr), c = (c1, …, cr−1) difúzní toky j1(x,t) +…+ jr(x,t) = 0† j = (j1, …, jr−1) zachování hmotnosti cۤ− vc’ + j’ = 0 cL = cR − CR/ v Cis theintegral ofcۤfromxLtox vyčíslení Gibbsovy energie: integrály z c*∙μ(c*),meze od xLdoxR chemické potenciály μ(c*) = wf μf(c*), součet přesf{α, β, γ} μf(c*) = μ0f + RT log c* + φf(c*) †jen pro substituční složky, pro intersticiální složky bez omezení

10. Onsagerův princip maximální rychlosti disipace: chemické potenciály váhové funkce Josiah Willard Gibbs 1839-1903 vazební podmínka pro (Gibbs-Duhem)

11. celková Gibbsova energie soustavy a její časová derivace zákon zachování hmotnosti pro difúzní toky (uzavřená soustava) rychlost disipace příslušná objemové difúzi a pohybu mezifázového rozhraní

12. Evoluční rovnice: Onsagerův termodynamický princip: kinetika soustavy odpovídá variaci po dosti dlouhých výpočtech vychází Lars Onsager 1903-1976

13. Stacionární řešení: rozklad chemických potenciálů podmínka stacionárního stavu integrace přes vektorově … (následuje mnoho formálních substitucí, podrobněji dále) vyjde r+s−1 rovnic pro r+s−1 & 1 rovnice pro konstantnírychlostv

14. Nestacionární řešení: evoluční rovnice pohybující se okrajové podmínky slabá formulace (testovací funkce  , malý kompaktní nosič) & 1 rovnice pro rychlostv(t)

15. Numerické přístupy: h rychlost rozhraní v = (M/Ω) ∫ c*∙μ(c*) dx 0 podrobněji: … evoluční rovnice Nj + B(c)c’+ K(c)c = 0 B(c)c’+ [K(c)+ (N/Ω)v]c = (N/Ω)(vcL + C) metoda přímek, Crankovo-Nicholsonovo schéma (zjednodušený tvar, ekvidistantní síť, s od 1 do m) Bs(cs–cs–1) + ½ δ[Ks+ (Ns/Ω)v] (cs+cs–1) = δ (Ns/Ω) (vc0 + Cs) & rychlost rozhraní v numerickou integrací (Simpsonovo pravidlo), obdobně též Cs & spektrální analýza (Fourierova metoda), zejména ve fázích α , γ

16. podrobněji: chemické potenciály stručné označení intersticiální složky substituční složky

17. Očekávatelná otázka:Jak se určí chemické potenciály  ? Pro čistě substituční slitinu (podle teoretického i experimentálního bádání kolegů z Leobenu) Fe-Cr-Ni s převahou Fe (tedy poněkud vylepšenou ocel) je lze vyjádřit (po předzpracování MAPLEm) např. pro speciální teplotu 1030 K takto: function muN=admui030(c);c1=c(1); c2=c(2); c3=c(3);p=1043; p1=-1354.5; p2=-468; p12=2373; p13=1100; p113=1100; p123=550; pX112=617;gp=p+p1*c1+p2*c2+p12*c1*c2+p13*c1*c3+p113*c1^2*c3+p123*c1*c2*c3+pX112*c1*c2*(c1-c2);gp1=p1+p12*c2+p13*c3+2*p113*c1*c3+p123*c2*c3+pX112*c2*(c1-c2)+pX112*c1*c2;gp2=p2+p12*c1+p123*c1*c3+pX112*c1*(c1-c2)-pX112*c1*c2;gp3=p13*c1+p113*c1^2+p123*c1*c2;gpa=gp1*c1+gp2*c2+gp3*c3;q=3.22; q1=-2.228; q2=-1.37; q12=4; q13=-0.85;gq=q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3;gq1=q1+q12*c2+q13*c3;gq2=q2+q12*c1;gq3=q13*c1;gqa=gq1*c1+gq2*c2+gq3*c3;a12=7241.75; a13=10529.6; a23=-2282.51; aM12=17670.5; aM32=463.152; aS5=-4.74062e-013; aS15=-1.11983e-044; aS25=-1.74984e-075;

18. m1=(-2*aM12*c2+a13)*c1^2+(2*aM32*c2*c3-aM32*c2^2-a12*c2+a23*c2+c2*a13+2*aM12*c2^2+2*aM12*c2-a13)*c1-aM32*c2^3-a23*c2+3*aM32*c2^2*c3+a23*c2^2-aM12*c2^2+a12*c2+a13*c3-aM32*c2*(c3-c2)-aM32*c2*c3;m2=(-2*aM12*c2+a13+aM12)*c1^2+(2*aM32*c2*c3-aM32*c2^2-a12*c2+a23*c2+c2*a13+2*aM12*c2^2+a12-2*aM12*c2-a13)*c1-aM32*c2^3-2*aM32*c2*c3+3*aM32*c2^2*c3+a23*c2^2+a23*c3-a23*c2+aM32*c3*(c3-c2)-aM32*c2*(c3-c2);m3=(-2*aM12*c2+a13)*c1^2+(2*aM32*c2*c3-aM32*c2^2-a12*c2+a23*c2+c2*a13+2*aM12*c2^2)*c1+3*aM32*c2^2*c3+a23*c2^2-aM32*c2^3;m1A=aS5*gp^5*log(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)+5*aS5*gp^4*(gp1-gpa)*log(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)+aS5*gp^5/(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)*(q1+q12*c2+q13*c3-c1*(q1+q12*c2+q13*c3)-c2*(q2+q12*c1)-q13*c1*c3);m2A=aS5*gp^5*log(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)+5*aS5*gp^4*(gp2-gpa)*log(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)+aS5*gp^5/(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)*(q2+q12*c1-c1*(q1+q12*c2+q13*c3)-c2*(q2+q12*c1)-…m1=(-2*aM12*c2+a13)*c1^2+(2*aM32*c2*c3-aM32*c2^2-a12*c2+a23*c2+c2*a13+2*aM12*c2^2+2*aM12*c2-a13)*c1-aM32*c2^3-a23*c2+3*aM32*c2^2*c3+a23*c2^2-aM12*c2^2+a12*c2+a13*c3-aM32*c2*(c3-c2)-aM32*c2*c3;m2=(-2*aM12*c2+a13+aM12)*c1^2+(2*aM32*c2*c3-aM32*c2^2-a12*c2+a23*c2+c2*a13+2*aM12*c2^2+a12-2*aM12*c2-a13)*c1-aM32*c2^3-2*aM32*c2*c3+3*aM32*c2^2*c3+a23*c2^2+a23*c3-a23*c2+aM32*c3*(c3-c2)-aM32*c2*(c3-c2);m3=(-2*aM12*c2+a13)*c1^2+(2*aM32*c2*c3-aM32*c2^2-a12*c2+a23*c2+c2*a13+2*aM12*c2^2)*c1+3*aM32*c2^2*c3+a23*c2^2-aM32*c2^3;m1A=aS5*gp^5*log(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)+5*aS5*gp^4*(gp1-gpa)*log(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)+aS5*gp^5/(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)*(q1+q12*c2+q13*c3-c1*(q1+q12*c2+q13*c3)-c2*(q2+q12*c1)-q13*c1*c3);m2A=aS5*gp^5*log(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)+5*aS5*gp^4*(gp2-gpa)*log(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)+aS5*gp^5/(q+q1*c1+q2*c2+q12*c1*c2+q13*c1*c3)*(q2+q12*c1-c1*(q1+q12*c2+q13*c3)-c2*(q2+q12*c1)-… … a následuje ještě asi 1500 obdobných programových řádků v MATLABu. Složitější případy ponecháme na fantazii laskavého čtenáře.

19. Výpočty s podporou softwaru MATLAB aMAPLE: 1.vypočtou se všechny molární frakce cs(pro celá s od 0 do m) v každém časovém kroku: nastaví se molární frakce cs×na základě počátečních podmínek (v prvním časovém kroku t = τ) nebo na základě hodnot získaných v předešlém časovém kroku t−τ, poté se odhadnou cs jako cs× 2.numericky se vypočte rychlost pohybu rozhraní v Simpsonovým pravidlem 3.nastaví se Bs, Ks a Nspomocí cs, dále se nastaví c0 jako sloupcový vektor r+s–1 neznámých parametrů (v programovém kódu lze využít komplexní aritmetiky) 4.Postupně se řeší cs (vždy jen ze soustavy r+s–1lineárních algebraických rovnic se stejným počtem neznámých), poté se vypočtou (stále jako jisté funkce c0) cs a ċs as (cs– cs*) / τ ,kde cs* pocházejí z předchozího časového kroku 5.určí se c0, přepočtou se všechna cs, použije se vhodný indikátor chyby (zachycující změny v hodnotách v i všech cs) 6. numericky se vypočtou integrályCs 7.rozhodne se o ukončení iterací, jinak se vrátí ke kroku 2

20. Realističtější výpočty: • stejnoměrná síť jen uvnitř rozhraní, nestejnoměrná síť v ostatních fázích • xL a xR aktualizovány v každém časovém kroku, nejsou tedy obecně identické s uzly sítě (je nutná interpolace) • Nedostatečná kvalita interpolačních schémat – kubické splajny používány namísto lineárních pro molární frakce, lineární splajny namísto konstant na intervalech (Crank-Nicholson) pro materiálové charakteristiky, alternativně i exponenciální funkce (blízké semianalytickému řešení Fourierovou metodou vně rozhraní) • materiálové charakteristiky získávány jako složité funkce od kolegů z Leobenu (polynomy až 25. řádu, vynásobené logaritmy dalších polynomů, apod.) – podpora softwaru MAPLE (resp. toolboxu MATLABu symbolic)

21. Existenční a konvergenční analýza: Fan & al. (1999) – algoritmus pro binární slitiny Mayer &.Simonett (1999) – existence of klasického řešení Schneider &. Inden (2004) – algoritmus pro vícesložkové slitiny, ostré rozhraní, umělé okrajové podmínky V. (2007) – předběžný výsledek pro stacionární případ (reflexivní Sobolevovy prostory, Gronwallovo-Bellmanov lemma);obecněji?časová diskretizace, Rotheho posloupnosti počáteční molární frakce Cr-Ni-Fe Příklad výpočtové simulace:

22. Vývoj molárních frakcípro Cr a Ni v čase

23. Molární frakce pro různé difúzní vlastnosti vynásobeno 1.1 (čárkovaně) či 0.9 (čerchovaně)

24. Molární frakce pro různé mobility rozhraní vynásobeno 1.15 (čárkovaně) či 0.85 (čerchovaně)

25. Chemické potenciály celkové potenciály, první aditivní část (čárkovaně), „linearizovaný”přístup (čerchovaně)

29. Skokyμi(h) −μi(0) for i {1,2,3} jako funkce tloušťky rozhraní Fe Cr Ni J. Svoboda, J. Vala, E. Gamsjäger, F. D. Fischer: A thick-interface model for diffusive andmassivephase transformation in substitutional alloys, Acta Materialia 54 (2006), 3953-3960, available atwww.sciencedirect.com

30. Rychlost pohybu rozhraní v jako funkce: i) tloušťky rozhraní h ii) teploty T

31. Užitečná zobecnění: • realistické 1D výpočty pro vícesložkové soustavy se substitučními i intersticiálními složkami, obecný nestacionární případ (dosud otevřené problémy: nejistá data, komplikované chemické potenciály, není garantováno žádné konečné „steady-state solution“…) • zahrnutí elastických a creepových přetvoření, následně rychlosti smršťování a roztahování (objemových změn), podle mechanismu b) • zahrnutí různých dalších toků, jmenovitě i) toku částinc vyvolaný gradientem teploty (Soretův jev) a ii) transportu tepla poplatného gradientu koncentrace (Dufourův jev) • obtíže 2D and 3D modelování: i) další hnací síly související se zakřivením rozhraní, ii) difúze podél (nejen kolmo na) rozhraní, iii) možná nerovnoměrná migrace rozhraní (nekonstantní tloušťka rozhraní ?) → FDM, FEM, FVM, …

32. Tento výzkum podporuje Grantová agentura AV ČR, reg.č. IAA200410601 (2006-08). Tuto cestu nesponzorovala GA AV ČR. DĚKUJI VÁM ZA POZORNOST. Ještě jedna velmi nekomerční reklama: 7. MATEMATICKÝ WORKSHOP S MEZINÁRODNÍ ÚČASTÍ FAST VUT v Brně – 16.října 2008 http://math.fce.vutbr.cz/konference