4. Symulacje numeryczne i weryfikacja doświadczalna
4.2. Parametry modelowania
4.2.2. Dobór parametrów siatek obliczeniowych
Metody numeryczne rozwiązywania równań różniczkowych obarczone są błędem [73, 74]. Błąd ten związany jest z zamianą przyrostów nieskończenie małych zmiennych całkowania na przyrosty skończone. Ze zmniejszeniem kroku dyskretyzacji zmiennych rozwiązanie numeryczne równania różniczkowego dąży do rozwiązania analitycznego [75
21
]. Wymaga to jednak większego zapotrzebowania na pamięć operacyjną oraz wydłuża czas obliczeń. Typowe wielkości komórek AK spotykane w literaturze [ , 27, 39, 40] mieszczą się w przedziale od 0.8 do 2 μm. W prezentowanej pracy zdecydowano się na stosowanie kroku a sieci gęstej (pole dyfuzji i AK) równego 1 μm.
Tab. 4.2. Termofizyczne parametry użyte w modelowaniu
Współczynnik przewodzenia ciepła: W/(m⋅K) Źródło
– cieczy λL 30 [41]
– austenitu λγ 20 [76
–grafitu
]
λgr 20 [76]
Współczynnik dyfuzji węgla w: m2 /s
– cieczy DL 1.25⋅10-9 [77
– austenicie
]
Dγ 5⋅10-10 [41]
Entalpia przemiany fazowej: J/m3
– cieczy w austenit ∆HL→γ 19.71⋅108 [76] – cieczy w grafit ∆HL→gr 16.16⋅105 – austenitu w grafit ∆Hγ→gr 8.8⋅105 Ciepło właściwe: J/(m3⋅K) – cieczy cv,L 5.6⋅106 [76] – austenitu cv,γ 5.84⋅106 [41] – grafitu cv,gr 17.84⋅105 [76] Gęstość: kg/m3 – cieczy ρL 7⋅103 [78 – austenitu ] ργ 7.3⋅103 [41] – grafitu ρgr 2.23⋅103 [76]
Współczynnik Gibbsa-Thomsona na granicy: m⋅K
– austenit-ciecz Γγ/L 1.9⋅10-7 [41]
– grafit-ciecz Γgr/L 7⋅10-6
– grafit-austenit Γgr/γ 10.5⋅10-6
Amplituda zmiany energii międzyfazowej na granicy:
– austenit-ciecz δΓ,γ/L 0.04 [40]
– grafit-ciecz δΓ,gr/L 0
– grafit-austenit δΓ,gr/γ 0
Średni kinetyczny współczynnik wzrostu: m/(s⋅K)
– austenitu w cieczy μ0,γ/L 10-3 [40]
– grafitu w cieczy μ0,gr/L 10-8
– grafitu w austenicie μ0,gr/γ 10-8
Amplituda zmiany kinetycznego współczynnika wzrostu:
– austenitu w cieczy δμ,γ/L 0.1 [40]
– grafitu w cieczy δμ,gr/L 0
– grafitu w austenicie δμ,gr/γ 0
Rząd symetrii kryształu:
– austenitu mS 4 [40]
Na etapie testowania prezentowanego modelu przeprowadzono obliczenia dla znacz-nie mznacz-niejszego kroku siatki tj. a = 0.15 μm. Opis parametrów modelowania jak i rezultaty obliczeń znajdują się w pracy [22]. Cześć tych wyników została przedstawiona na rys. 4.1 a-d. Stosowanie tak małego kroku ma jednak zaletę: dokładniej zachowywany jest założo-ny kąt orientacji krystalograficznej. W przedstawiozałożo-nym przykładzie założono go na
po-ziomie 22.5 stopnia względem pionowej osi układu współrzędnych. Minusem jest mniej-szy dopuszczalny krok czasowy oraz zwiększenie ilości komórek, niezbędnej dla modelo-wania obszaru fizycznego o tej samej wielkości, jak w przypadku stosowania siatki z dużym krokiem przestrzennym. Opis wpływu wielkości kroku siatki AK na kierunek wzrostu dendrytu można znaleźć np. w pracy [40].
a) b)
c) d)
Rys. 4.1. Wzrost dendrytu austenitu przy małym kroku siatki obliczeniowej a = 0.15 μm
Innym czynnikiem wpływającym na zachowanie się założonego kierunku wzrostu dendrytu austenitu jest szybkość odprowadzania ciepła z układu. Przykłady wzrostu den-drytu austenitu dla różnych szybkości chłodzenia pokazano na rys. 4.2. Obliczenia wyko-nano dla siatki AK o długości boku a = 1 μm. Dla każdego z czterech ziaren założono kąt θ0 orientacji krystalograficznej (kierunek najszybszego wzrostu) równy 22.5° względem pionowej osi.
I 18 s 12 s 7 s 5 s II 52 s 25 s 16 s 10 s III 86 s 38 s 25 s 15 s IV 120 s 50 s 33 s 20 s a) b) c) d)
Rys. 4.2. Wpływ szybkości chłodzenia na zachowanie się kierunku wzrostu głównych gałęzi dendrytu dla stężenia węgla 4.25%. Poszczególne kolumny oznaczają następującą szybkość
chłodzenia: a) 0.25 K/s; b) 0.5 K/s; c) 1 K/s; d) 2 K/s. Czasy symulacji umieszczono pod rysunkami. Długość boku komórki AK równa jest 1 μm
Dla szybkości chłodzenia równej 0.25 K/s kierunek wzrostu gałęzi dendrytu nie-znacznie odbiega założonego (kolumna a na rys. 4.2). Przy dwukrotnie większej prędkości chłodzenia gałęzie zaczynają zmieniać kierunek wzrostu przy pewnej wielkości dendrytu (rys. 4.2 II-b). Efekt ten jest jeszcze lepiej widoczny na rys. 4.2 II-c, gdzie prędkość chło-dzenia jest czterokrotnie większa. Dalsze zwiększanie szybkości chłodzenia do wartości 2 K/s powoduje, że główne ramiona dendrytu rosną wzdłuż przekątnej już w przypadku bardzo małego ziarna (rys. 4.2 I-d). Przy tej szybkości chłodzenia pojawiają się gałęzie
wtórne. Rosną one niestety tylko po jednej stronie gałęzi pierwotnych. Tego niekorzystne-go efektu niestety nie udało się dotychczas wyeliminować. Symetryczny wzrost wtórnych gałęzi udaje się uzyskać w przypadku, gdy założy się kierunek wzrostu po przekątnej siat-ki. Dlatego w symulacjach przeprowadzonych w tej pracy zakładano kąt równy 45°.
Kolejnym czynnikiem wpływającym na przyciąganie głównych gałęzi dendrytów do kierunku przekątnej komórek jest początkowe stężenie węgla. Na rys. 4.3 pokazano przy-kład symulacji wzrostu dendrytu austenitu dla stężenia węgla równego 0.6 %wag. Dla każ-dego z obliczeń założono różne szybkości chłodzenia tj. 0.5, 1.5, 2.5 i 3 K/s. Reszta para-metrów modelowania była taka sama (oprócz stężenia węgla) jak w obliczeniach z rys. 4.2.
I 9 s 4 s 2 s 2 s II 18 s 6 s 5 s 3 s III 29 s 9 s 7 s 5 s IV 40 s 13 s 8 s 6 s a) b) c) d)
Rys. 4.3. Wpływ szybkości chłodzenia na kierunek wzrostu dla stężenia węgla 0.6%. Poszczególne kolumny oznaczają następującą szybkość chłodzenia: a) 0.5 K/s; b) 1.5 K/s; c) 2.5 K/s; d) 3.0 K/s.
Z porównania rys. 4.2 i 4.3 widać, że w przypadku mniejszego stężenia węgla zało-żony początkowy kierunek wzrostu zachowuje się lepiej. Ze zwiększeniem szybkości chłodzenia od 0.5 do 3 K/s zmienia się kąt kierunku wzrostu głównych ramion dendrytu (od 22.5 do 45 stopni). Należy zwrócić uwagę, że zmiana kierunku wzrostu następuje wol-niej, niż w przypadku stężenia węgla równego 4.25% wag., i zaczyna się ona przy więk-szych szybkościach chłodzenia. Pokazuje to, że zdolność modelu na zachowanie się zało-żonego kierunku wzrostu gałęzi dendrytów zmniejsza się ze zwiększeniem stężenia węgla. Dopuszczalne wartości kroków czasowych dla pola temperatury i stężenia, wyliczo-ne z równań (3.12) i (3.16) na podstawie danych z tabeli 4.1 różnią się od siebie o kilka rzędów wielkości. Kroki te różnią się też między sobą w zależności od faz (różne wartości współczynników termofizycznych). W celu zapewnienia stabilności obliczeń należy wy-brać najmniejszą dopuszczalną wartość i stosować ją dla wszystkich faz występujących na siatce. Dla pola temperatury fazą, która wymaga zastosowania najmniejszego kroku cza-sowego, jest grafit z uwagi na najmniejszy stosunek ciepła właściwego do współczynnika przewodzenia ciepła. Natomiast dla pola stężenia fazą taką jest ciecz. Mając założony krok siatki dla pola stężenia (a = 1 μm), można tak dobrać wielokrotność nk kroku siatki tempe-raturowej, aby dopuszczalne kroki czasowe dla obu siatek były bliskie. W tym celu należy porównać ze sobą stronami zależności (3.12) i (3.16) i wyliczyć stosunek kroku czasowego dla pola temperatury do pola stężenia. W ten sposób otrzymuje się następującą zależność
L gr v gr L C gr T k D c n , , max, , max, λ = τ ∆ τ ∆ = (4.2)
Wartość nk, wyliczona z powyższej zależności, wynosi w przybliżeniu 94.7. W obliczeniach przyjęto nieco mniejszą wartość nk = 80.
Na rys. 4.4 przedstawiono rozkład temperatury interpolowanej z sieci rzadkiej na sieć gęstą dla czterech wartości nk tj. 2, 5, 20 i 80. Wszystkie obliczenia wykonano dla tych samych warunków. Zarodek austenitu umieszczono w środku siatki obliczeniowej o rozmiarach 240×240 komórek AK i długości boku sieci gęstej a = 1 μm. Obszar chło-dzono zewnętrznym strumieniem, który w przypadku braku przemiany zapewnia obniżanie się temperatury z szybkością 10 K/s. Wyniki przedstawiono dla symulowanego czasu kry-stalizacji ok. 2.5 s.
a) b)
c) d)
Rys. 4.4. Rozkład temperatury dla różnych wartości współczynnika nk: a) 2; b) 5; c) 20; d) 80