Index of /rozprawy2/10376

im. Stanisława STASZICA

w Krakowie

WYDZIAŁ ODLEWNICTWA

Rozprawa doktorska

Model krystalizacji

w układzie dwuskładnikowym

z

udziałem trzech faz z wykorzystaniem automatu

komórkowego

mgr inż. Daniel Gurgul

Promotor:

dr hab. inż. Andriy Burbelko, prof. AGH

promotorskiego N507 422536 finansowanego przez

MNiSW

Spis treści

Wykaz ważniejszych oznaczeń ... 6

Znaczenie ważniejszych indeksów ... 8

Streszczenie ... 9

Wstęp ... 10

1. Przegląd metod modelowania krystalizacji żeliwa ... 12

1.1. Modele analityczne ... 12

1.2. Komputerowe modele obliczeniowe ... 13

1.2.1. Modele numeryczne z analitycznym opisem przepływu ciepła ... 13

1.2.2. Modele numeryczne ... 14

2. Cel i zakres pracy ... 17

3. Opis modelu ... 18

3.1. Definicja automatu komórkowego ... 18

3.2. Typy granic międzyfazowych w stosowanym AK ... 19

3.3. Transport masy... 22

3.4. Przepływ ciepła ... 26

3.4.1. Zmiana pola temperatury ... 26

3.4.2. Rozwiązanie numeryczne ... 26

3.4.3. Dopuszczalny krok czasowy rozwiązania numerycznego ... 28

3.4.4. Interpolacja temperatury ... 30

3.5.1. Znane rozwiązania ... 31

3.5.2. Testowanie równań ... 33

3.6. Kierunek normalny do frontu ... 36

3.6.1. Znane metody ... 36

3.6.2. Testowanie metod ... 38

3.7. Krzywizna frontu ... 41

3.7.1. Znane metody ... 41

3.8. Prędkość przemieszczania się frontu ... 49

3.9. Zarodkowanie ... 52

3.9.1. Modelowanie zarodkowania ciągłego ... 52

3.9.2. Algorytm wyboru lokalizacji podkładek ... 54

3.9.3. Algorytm aktywacji podkładek ... 56

3.9.4. Algorytm równoczesnego wyboru lokalizacji i warunków aktywacji podkładek ... 57

3.9.5. Prawo zarodkowania ... 58

3.9.6. Porównanie działania algorytmów zarodkowania ... 58

3.10. Reguły zmiany stanów komórek AK ... 61

3.11. Warunki brzegowe i początkowe ... 63

3.12. Algorytm obliczeniowy ... 64

4. Symulacje numeryczne i weryfikacja doświadczalna ... 66

4.1. Opis wytopów ... 66

4.1.1. Analiza metalograficzna ... 69

4.2. Parametry modelowania ... 73

4.2.2. Dobór parametrów siatek obliczeniowych ... 73

4.3. Obliczenie dla stopu o składzie podeutektycznym ... 79

4.3.1. Modelowanie komputerowe ... 79

4.3.2. Jakościowa porównanie wyników ... 83

4.4. Obliczenie dla stopu o składzie eutektycznym ... 86

4.5. Obliczenia dla stopu o składzie nadeutektycznym ... 89

4.6. Jakościowa weryfikacja na podstawie modelowych substancji krystalizacji ... 92

5. Podsumowanie i wnioski ... 94

Wykaz ważniejszych oznaczeń

a – długość boku kwadratowej komórki automatu komórkowego, m

T

a – średni współczynnik wyrównywania temperatury, m2/s b – współczynnik zarodkowania, K

C – stężenie, uł. wag., kg/m3, % wag.

cv – ciepło właściwe przypadające na jednostkę objętości, J/(m3⋅K)

v

c – średnie ciepło właściwe przypadające na jednostkę objętości, J/(m3⋅K) D – współczynnik dyfuzji domieszki, m2/s

f – udział fazy stałej w komórce KInt – kod konfiguracji interfejsu

k – równowagowy współczynnik rozdziału domieszki pomiędzy fazą stałą a ciekłą

lp, ln – wektory poprzedni i następny wykorzystywane w metodzie geometryczno-wektorowej wyznaczania krzywizny, m

ms – rząd symetrii kryształu

Na – średnia ilość aktywnych podkładek zarodkowania w jednostkowej objętości, 1/m3

Na,V – średnia ilość aktywnych podkładek zarodkowania w objętości V stopu NF – średnia ilość aktywnych podkładek zarodkowania przypadająca na

jednost-kę powierzchni, 1/m2

NF,n – średnia ilość kulek grafitu przypadająca na jednostkę powierzchni, 1/m2 NV – średnia ilość podkładek zarodkowania w jednostkowej objętości, 1/m3

NS – średnia ilość podkładek zarodkowania na jednostkowej powierzchni, 1/m2

n – jednostkowy wektor normalny do interfejsu międzyfazowego

nk – wielokrotność kroku sieci rzadkiej w stosunku do kroku sieci gęstej nx, ny – współrzędne wektora n

p – liczba losowa o rozkładzie jednostajnym z zakresu [0...1) qT – wydajność źródła ciepła, W/m3

T – temperatura, K

TEq – temperatura równowagi fazowej, K, °C

TS – średnia temperatura w obszarze obliczeniowym lub temperatura na krzywej stygnięcia, K, °C

Ṫ – szybkość zmiany temperatury (stygnięcia), K/s t – numer kroku czasowego

u – wartość prędkości przemieszczania się frontu, m/s

u – wektor prędkość frontu, m/s

v – objętość komórki AK, m3 Δf – przyrost udziału fazy stałej ∆H – entalpia przemiany fazowej, J/m3

ΔT – lokalne przechłodzenie, K

ΔTźr – zmiana temperatury spowodowana wydzielaniem się ciepła przemiany fa-zowej, K

ΔTκ – przechłodzenie kapilarne, K

ΔTμ – przechłodzenie kinetyczne na froncie, K Δτ – krok czasowy, s

Δτmax,C – dopuszczalny krok czasowy w modelowaniu pola stężenia, s Δτmax,T – dopuszczalny krok czasowy w modelowaniu pola temperatury, s

Γ – współczynnik Gibbsa-Thomsona, K⋅m δΓ – amplituda zmiany energii międzyfazowej

δμ – amplituda zmiany kinetycznego współczynnika wzrostu θ – kąt pomiędzy wektorem n a osią X, rad

θ0 – kąt orientacji krystalograficznej, rad, stopnie ξ – bazowe funkcje interpolacji

κ – krzywizna frontu, 1/m

λ – współczynnik przewodzenia ciepła, W/(m⋅K)

λ – średni współczynnik przewodzenia ciepła, W/(m⋅K) μ – kinetyczny współczynnik wzrostu, m/(s⋅K)

μ0 – wartość średnia kinetycznego współczynnika wzrostu, m/(s⋅K) ρ – gęstość, kg/m3

σ – odchylenie standardowe τ – czas, s

∇ – operator nabla ˄ – i (koniunkcja)

Znaczenie ważniejszych indeksów

gr, L, γ – fazy, odpowiednio grafit, ciecz i austenit

i, j – współrzędne środka komórki na sieci automatu komórkowego lub pola stę-żenia

ir, jr – współrzędne węzła na sieci rzadkiej

N, E, S, W – oznaczenia współrzędnych sąsiednich komórek, komórki o współrzędnych (i,j): N = (i,j-1), E = (i+1,j), S = (i,j+1), W = (i-1,j). To samo dotyczy indek-sów ir, jr

Streszczenie

W pierwszej części pracy przedstawiono przegląd metod modelowania krystalizacji żeliwa. Przegląd ten ma charakter krótkiego rysu historycznego począwszy od modeli ana-litycznych a skończywszy na modelach w pełni numerycznych. Opisano w skrócie cechy charakterystyczne poszczególnych modeli. W trakcie opisu zaprezentowano najważniejsze prace, jakie ukazywały się od lat 60-tych ubiegłego wieku.

W rozdziale „Opis modelu” przedstawiono szczegółowy opis opracowanego modelu matematycznego krystalizacji żeliwa z grafitem sferoidalnym. Model ten oparty jest na technice automatu komórkowego (AK), który wykorzystywany jest do symulacji tworzenia się struktury żeliwa podczas krystalizacji. Jego cechą charakterystyczną jest to, że kształt rosnących ziaren austenitu i grafitu, nie jest z góry założony, lecz jest wynikiem obliczeń modelowych. W modelu brane są pod uwagę takie zjawiska jak: przepływ ciepła, dyfuzja węgla w cieczy i austenicie, wydzielanie się ciepła przemiany fazowej i jego wpływ na warunki panujące na froncie krystalizacji, zarodkowanie ziaren zarówno austenitu, jak i grafitu. W modelu uwzględniono nierównowagowy charakter przemian fazowych oraz wpływ krzywizny frontu na temperaturę równowagi termodynamicznej. Uwzględniane zjawiska ujmowane są w skali mezo*

W części doświadczalnej wykonano modelowanie krystalizacji żeliwa o składzie podeutektycznym, eutektycznym oraz nadeutektycznym. Wyniki obliczeń skonfrontowano z wynikami doświadczalnymi uzyskanymi zarówno dla żeliwa z grafitem kulkowym, jak i dla tzw. przeźroczystych substancji modelowych. Porównano ze sobą rzeczywistą krzywą chłodzenia z krzywą uzyskaną podczas obliczeń. W celu jakościowego porównania struk-tury, a w szczególności wielkości i dystrybucji kulek grafitu, zestawiono ze sobą mikro-strukturę odlewów doświadczalnych i komputerowe mapy modelowanej struktury. Analiza wyników potwierdza poprawność założonej tezy pracy.

.

* _{Mezo w tym przypadku oznacza skalę wymiarową pomiędzy wielkością ziarna a odległością}

Wstęp

Żeliwo sferoidalne jest materiałem o bardzo dużym znaczeniu technologicznym. Jest używane do odlewania detali o wadze od kilkuset gramów do nawet ponad stu ton. Wyko-rzystywane jest ono m.in. w przemyśle motoryzacyjnym, metalurgicznym, energetycznym, w tym dla energetyki odnawialnej, chemicznym i petrochemicznym. Jego udział w świa-towej produkcji odlewów w stosunku do innych stopów za rok 2009 przedstawiony jest na poniższym wykresie kołowym (rys. 1) [1].

Rys. 1. Światowa produkcja odlewów w 2009 roku wg [1]

Jak widać z tego rysunku, żeliwo sferoidalne w 2009 roku miało drugi co do wielko-ści udział w ogólnej ilowielko-ści produkcji odlewów. Jak podaje to samo źródło, wzrost jego pro-dukcji w dekadzie poprzedzającej 2009 rok wynosił 39.6 %. Dla porównania, produkcja odlewów staliwnych w tym samym okresie wzrosła o 50%, odlewów ze stopów metali nieżelaznych o 41.1% a z żeliwa szarego tylko o 9.3%.

Produkcja odlewów z żeliwa sferoidalnego (Ductile Iron) jako jedyna wskazuje na stały wzrost od 2004 r. (rys. 2), podczas gdy dla odlewów staliwnych (Steel) od 2006 stwierdza się zahamowanie wzrostu, a dla aluminium (Aluminium) i żeliwa szarego (Gray Iron) od roku 2007 notuje się spadek produkcji.

Rys. 2. Produkcja odlewów w świecie z kilku rodzajów stopów w okresie od 2004 do 2008 roku

wg [2

Tak duża produkcja odlewów z tego stopu związana jest z kilkoma jego cechami. Jest to materiał stosunkowo tani w otrzymywaniu, odznacza się bardzo dobrymi właściwo-ściami odlewniczymi np. lejnością, posiada bardzo dobre właściwości wytrzymałościowe i plastyczne. Technologia jego otrzymywania jest znana od lat i nie nastręcza większych problemów; zwłaszcza sam zabieg sferoidyzacji jest dobrze opanowany.

]

Ciekawe zastosowanie żeliwa sferoidalnego można zaobserwować w ostatnim czasie w wykorzystaniu go do prób wykonywania odlewów tzw. supercienkościennych [3, 4

Przedstawienie znaczenia, jakie odgrywa żeliwo sferoidalne we współczesnym świa-towym odlewnictwie, ma na celu uzasadnienie wyboru stworzenia modelu matema-tycznego krystalizacji eutektycznej na przykładzie właśnie tego materiału.

]. Chociaż gęstość żeliwa jest prawie trzykrotnie większa, niż na przykład stopów na bazie aluminium, to ma ono większą wytrzymałość w przeliczeniu na jednostkę masy, niżeli stopy nieżelazne. Dlatego odlewy ze stopu aluminium można w niektórych przypadkach zastąpić w pełni funkcjonalnymi bardziej wytrzymałymi, a równocześnie lżejszymi odpo-wiednikami z żeliwa sferoidalnego, zmniejszając w odpowiedni sposób grubość ścianek.

1.

Przegląd metod modelowania krystalizacji

żeliwa

1.1. Modele analityczne

Jeden z pierwszych modeli analitycznych, który w sposób ilościowy opisuje krystali-zację eutektyki płytkowej i włóknistej można znaleźć w pracy [5

Ograniczenia związane z stosowaniem modelu tylko do krystalizacji kierunkowej i warunków izotermicznych, panujących na płaskim froncie, zostały pokonane w pracy [

]. Praca ta wprawdzie nie dotyczy bezpośrednio żeliwa i została zweryfikowana na przykładzie krystalizacji orga-nicznych przeźroczystych substancji modelowych oraz dwufazowego stopu Pb-Sn, jednak w literaturze specjalistycznej uważana jest za jedną z przełomowych. Zaproponowany mo-del krystalizacji eutektyki płytkowej i włóknistej jest oparty na rozwiązaniu równania dy-fuzji składnika w ruchomym układzie współrzędnym dla przypadku krystalizacji kierun-kowej. W modelu zakłada się, że front krystalizacji znajduje się w warunkach izotermicz-nych i ma płaski kształt. Podejście takie sprawdza się jedynie dla przypadku krystalizacji kierunkowej.

6

Pierwszy analityczny model opisujący wzrost eutektyki w żeliwie sferoidalnym zo-stał zaproponowany w 1972 roku przez Wetterfalla i współautorów [

]. Zaproponowano model sprzężonego wzrostu eutektyki płytkowej w sferycznym ziarnie eutektycznym dla przypadku krystalizacji żeliwa szarego. Model ten potrafi przewidzieć zależność pomiędzy prędkością wzrostu a odległością międzypłytkową.

7

1.1

]. Pozwala on na obliczenie wzrostu otoczki austenitu (γ) wokół kulki grafitu (gr) z fazy ciekłej (L), co sche-matycznie zostało pokazane na rys. . Proces ten jest kontrolowany dyfuzją węgla w fazie stałej. Równanie opisujące zależność prędkości wzrostu otoczki austenitu ma postać

(

)

γ γ γ γ γ γ γ γ − − − = τ _{/ /} / / L L L gr gr gr C C C C r r r r D d dr (1.1) gdzie:

Dγ – współczynnik dyfuzji węgla w austenicie, C – stężenie węgla na odpowiednich granicach.

Znaczenie indeksów dolnych z powyższego równania wyjaśniono schematycznie na rys. 1.1.

Rys. 1.1. Schemat rozkładu stężenia węgla w otoczce austenitu i cieczy

Model ten przetrwał próbę czasu i w przedstawionej lub zmodyfikowanej formie używany jest do dzisiaj w większości rozwiązań, które pozwalają na obliczenia mikro-struktury np. w żeliwie z grafitem sferoidalnym. Stosowany jest on w modelach typu mi-kro, których jedną z cech jest m.in. to, że kształt rosnących ziaren jest z góry założony [8

1.2. Komputerowe modele obliczeniowe

-11].

Do tej grupy będą zaliczane modele, w których przynajmniej jedno z zadań wymaga zastosowania rozwiązania numerycznego za pomocą komputera. Reszta zadań może być rozwiązywana w sposób analityczny. Niektóre problemy są zbyt skomplikowane i czasochłonne do pokonania bez pomocy metod numerycznych. Przypadek taki ma miej-sce, gdy zachowanie się układu opisane jest równaniami różniczkowymi, które to dla danej sytuacji i warunków brzegowo-początkowych nie mają rozwiązania analitycznego. W sytuacji takiej niezbędne jest zastosowanie metod numerycznych rozwiązywania rów-nań różniczkowych, które ogólnie sprowadzają się do wykonywania dużej ilości powtarza-jących się obliczeń.

1.2.1. Modele numeryczne z

analitycznym opisem przepływu ciepła

Uważa się, że pierwszą pracą zaliczaną do tej grupy i zarazem przełomową był arty-kuł [12] opublikowany w 1966 roku przez W. Oldfielda. Opracował on komputerowy

mo-del pozwalający na wyliczenie krzywych chłodzenia (rys. 1.2) żeliwa szarego przy zadanej szybkości chłodzenia. Jest to pierwszy model pozwalający na prognozowanie mikrostruk-tury poprzez komputerowe obliczenia. Przedstawiono również wyniki weryfikacji do-świadczalnej tych obliczeń. Zgodnie z założeniami Oldfielda liczba ziaren eutektycznych zależy jedynie od stopnia przechłodzenia poniżej temperatury równowagi eutektycznej i nie zależy od czasu. Zaproponowano też równanie na prędkość rozrastania się kulistego ziarna eutektycznego. Prędkość ta jest proporcjonalna do kwadratu przechłodzenia.

Rys. 1.2. Eksperymentalne i obliczeniowe krzywe chłodzenia oraz równania opisujące liczbę

ziaren i prędkość wzrostu zaproponowane przez Oldfielda [12]

Kolejnymi pracami, o których warto wspomnieć, są artykuły Fredrikssona i Svenssona, opublikowane odpowiednio w 1985 [13] i w 1988 roku [14

1.2.2. Modele numeryczne

]. Autorzy tych prac stworzyli model, który rozwiązywał przepływ ciepła w sposób analityczny z parabolicznym prawem wzrostu eutektyki w żeliwie szarym i białym, uwzględniał dyfu-zję węgla w otoczce austenitu w żeliwie sferoidalnym. Fredriksson i Svensson przeanali-zowali wpływ na kształt obliczeniowej krzywej chłodzenia takich czynników jak: rodzaju krystalizującej eutektyki, szybkości chłodzenia oraz ilości podkładek zarodkowania.

Pierwszym modelem, który wykorzystywał metodą różnic skończonych do rozwią-zywania wszystkich równań podczas symulacji krystalizacji żeliwa sferoidalnego, był za-proponowany przez Su i in. w pracy [15]. Użyli oni modelu zarodkowania Oldfielda [ ] 12 oraz prawa wzrostu grafitu kontrolowanego dyfuzją węgla poprzez otoczkę austenitu. W pracy tej przeprowadzono walidację modelu poprzez porównanie obliczeniowej krzy-wej chłodzenia z krzywą eksperymentalną. Ocena wyliczonej wielkości i dystrybucji kulek

grafitu okazała się mało zgodna z eksperymentem. Jak sami autorzy podkreślali, ta część modelu wymagała poprawy.

Model Su po raz pierwszy został poprawiony przez Rappaza i współautorów [16

(1.1)

]. Uwzględniano tam zmianę objętości podczas krystalizacji eutektyki globularnej modyfiku-jąc równanie do postaci

(

)

(

)

(

)

1 / / 2 3 3 0 / / / 3 − γ γ γ γ γ γ γ γ γ γ γ γ ₋         − τ ρ ρ + − − = τ L L L L gr L gr gr C C r r r d dC r r r C C r D d dr (1.2) gdzie:

r0 – tzw. równoważny promień ziarna eutektycznego.

Kolejne poprawy tego równania można znaleźć w pracach [17, 18

Pierwszą próbę wizualizacji tworzenia się struktury żeliwa sferoidalnego podjęto w pracy [

].

19]. Model ten oparty jest na wcześniejszej pracy tych samych autorów [20

Kolejną generacją modeli matematycznych służących do modelowania krystalizacji żeliwa sferoidalnego są modele wykorzystujące automat komórkowy nie tylko do wizuali-zacji wyników, jak to ma miejsce we wspomnianej pracy [

]. Założono w nim, że każdy sferoid grafitu otoczony jest jedną otoczką austenitu tworząc jedno ziarno eutektyczne. Ziarno takie rozrastało się sferycznie do momentu zetknięcia się z innym ziarnem. Cały obszar obliczeniowy znajdował się w jednorodnym polu temperatu-ry. Ciepło z obszaru odprowadzane było na podstawie odpowiednich warunków założo-nych z góry. W celu wizualizacji wyników użyto automatu komórkowego z trójwymiaro-wą siecią sześciennych komórek.

19], ale do tworzenia się struktu-ry. Bazują one na wcześniejszych rozwiązaniach służących modelowaniu krystalizacji sto-pów za pomocą automatów komórkowych w układach dwufazowych typu: ciecz-faza sta-ła. W modelach tych, cechą charakterystyczną jest to, że kształt ziaren nie jest z góry na-rzucony. Jest on wynikiem obliczeń numerycznych.

Należy zaznaczyć, że modele tego typu są stosunkowo nowe i poza kilkoma publika-cjami w literaturze, nie spotyka się innych. Przykładem mogą być wyniki opublikowane w [21]. W pracy tej, przedstawiony jest skrótowo model krystalizacji żeliwa sferoidalnego z wykorzystaniem automatu komórkowego. Cechą charakterystyczną tego modelu jest zakładanie z góry, że w każdym momencie w komórkach dwufazowych panują warunki równowagowe. Oznacza to, że w przypadku np. komórki zawierającej austenit i ciecz, przyrost austenitu w danym kroku czasowym dobierany jest w taki sposób, aby stężenie

węgla w cieczy było równe wartości równowagowej. Temperatura jest jednakowa w całym obszarze obliczeniowym, a szybkość jej zmiany narzucona jest z góry. Przedstawiono tam wyniki obliczeń dla składu nadeutektycznego oraz podeutektycznego. Obliczeniowa krzy-wa chłodzenia została porównana z krzywą uzyskaną w eksperymencie.

Innymi pracami wykorzystującymi podobny aparat matematyczny, są prace opubli-kowane z udziałem autora niniejszej rozprawy [22-24].

2. Cel i zakres pracy

Głównym celem pracy jest opracowanie modelu matematycznego krystalizacji żeli-wa z grafitem sferoidalnym, algorytmu rozwiązywania równań modelowych oraz progra-mu komputerowego służącego do celów syprogra-mulacji dwuwymiarowej procesu krystalizacji tego stopu.

Kształt i ilość ziaren rosnących faz, szybkość zarodkowania i wzrostu ziaren oraz ich wzajemne rozmieszczenie przestrzenne nie będą zakładane z góry lecz stanowić będą wy-niki modelowania.

W modelowaniu wykorzystano metodę różnic skończonych (wyznaczenie pola tem-peratury i pola stężenia składnika rozpuszczonego) oraz technikę automatu komórkowego (modelowanie wzrostu i kształtu ziaren). Zastosowane w pracy metody i algorytmy podda-no testowaniu w celu zwiększenia dokładności rozwiązań.

Zasadniczą tezę pracy można sformułować w następujący sposób:

Prezentowany model umożliwia symulację krystalizacji dwufazowej w stopach o składzie eutektycznym, podeutektycznym i nadeutektycznych w przypadku, gdy jed-na z rosnących faz zostaje odizolowajed-na od cieczy prze drugą fazę stałą. Model uwzględnia stochastyczny charakter zarodkowania ziaren obu faz oraz nierównowa-gowy charakter procesów termodynamicznych wzrostu. Jest on poprawny jakościo-wo, z wynikami uzyskanymi w weryfikacji doświadczalnej, zarówno pod względem zgodności mikrostruktury jak i krzywych stygnięcia.

3. Opis modelu

W tej części pracy podano szczegółowy opis modelu matematycznego krystalizacji eutektycznej stopu dwuskładnikowego Fe-C na przykładzie żeliwa z grafitem sferoidal-nym. Przedstawiono znane metody obliczania kierunku normalnego frontu, krzywizny frontu i szybkości przemiany. Spośród tych metod, na podstawie zaproponowanych przez autora odpowiednich testów, wybrano najlepsze rozwiązania, wprowadzając w niektórych pewne modyfikacje lub proponując nowe.

3.1. Definicja automatu komórkowego

Zgodnie z [25, 26

− sieć komórek {i} przestrzeni D-wymiarowej,

] automat komórkowy (AK) jest pojęciem matematycznym, które można zdefiniować jako obiekt matematyczny, który składa się z trzech następujących składowych:

− zbiór {si} stanów pojedynczej komórki, zwykle ten sam dla wszystkich komórek, zawierający k elementów,

− reguła F określającą stan komórki w chwili t + 1 w zależności od stanu w chwili t tej komórki i komórek ją otaczających; si(t + 1) = F({sj(t)}), j ∈ O(i) gdzie O(i) jest oto-czeniem i-tej komórki.

Jeżeli funkcja F nie zależy od zmiennej losowej, to taki automat nazywa się determi-nistycznym, w przeciwnym przypadku nazywany jest probabilistycznym.

W niniejszej pracy użyto sześć stanów komórek k = 6. Zbiór {si} jest jednakowy dla każdej komórki AK i zawiera następujące elementy: {si}={0, 1, 2, 3, 4, 5}. Komórce przy-pisuje się element zbioru {si} w zależności od fazy, jaką zawiera. W danym kroku czaso-wym możliwa jest tylko jedna wartość. Poszczególnym komórkom w zależności od zawar-tości przypisuje się następujące warzawar-tości si:

0 – dla komórki w stanie ciekłym,

1 – komórka jest w stanie stałym (zawiera austenit – faza γ ), 2 – komórka jest w stanie stałym (zawiera grafit – faza gr),

4 – stan przejściowy pomiędzy cieczą a fazą stałą ze stanu 2,

5 – stan przejściowy pomiędzy fazą stałą ze stanu 1 a fazą stałą ze stanu 2. Stany od 3 do 5 w dalszej części pracy będą nazywane interfejsem.

Proponowany model nie uwzględnia tzw. punktu potrójnego, czyli komórki zawiera-jącej trzy fazy.

Każda komórka ma swoje otoczenie, przy czym ona sama nie jest zaliczana do swo-jego otoczenia. Znane są różne rodzaje otoczenia. Na przykład, jeśli do otoczenia komórki i zaliczane są komórki sąsiadujące z nią tylko wzdłuż jej boków (rys. 3.1 a), to takie oto-czenie nazywane jest otooto-czeniem von Neumana. Gdy do tego otoczenia dołączane są ko-mórki stykające się narożnikami, takie otoczenie nazywane jest otoczeniem Moore'a (rys. 3.1 b). Mówi się, że promień otoczenia wynosi odpowiednio a lub a 2, gdzie a jest dłu-gością boku kwadratowej komórki AK.

a) b)

Rys. 3.1. Schemat oznaczania sąsiednich komórek automatu komórkowego:

a) otoczenie von Neumana, b) otoczenie Moore'a

Oprócz przedstawionych powyżej otoczeń, w prezentowanej pracy w przypadku wy-liczania kierunku normalnego do frontu używane jest otoczenie o promieniu a 5.

3.2. Typy granic międzyfazowych w stosowanym AK

W celu wyznaczenia lokalnej konfiguracji geometrycznej granicy ziarna fazy stałej, do którego należy komórka interfejsu, wykorzystuje się następujący algorytm. Wszystkim 8 komórkom w otoczeniu Moore'a (rys. 3.1 b)danej komórki interfejsuprzypisuje się war-tości, stanowiące potęgę liczby 2, zgodnie ze schematem, pokazanym na rys. 3.2. Następ-nie kod konfiguracji interfejsu jest wyznaczany sumowaNastęp-niem

E

S

i

W

N

S-W

S-E

N-E

N-W

i

E

S

W

N

∑

= δ ⋅ = 8 1 i i i Int m K (3.1) gdzie

mi – wartość ze zbioru pokazanego na rys. 3.2,

δi – zmienna wynosząca 1, jeżeli sąsiednia komórka i znajduje się w stanie sta-łym i należy do tego samego ziarna, co komórka analizowana, lub 0 w każdym innym przypadku.

Ilość wszystkich kombinacji wynosi 256. Na podstawie kodu konfiguracji interfejsu można określić, po której stronie komórki interfejsu znajduje się faza stała, co ułatwia im-plementację algorytmów, np. algorytmu do wyznaczania krzywizny granicy międzyfazowej.

Rys. 3.2. Maska służąca do wyznaczania lokalnej konfiguracji interfejsu

Granice międzyfazowe, na podstawie kodu interfejsu, sklasyfikowano na następują-ce typy – płaskie, wypukłe, wklęsłe, cienkie kanaliki, dna kanalików, wnęki zamknięte oraz złożone. Przykład takiej klasyfikacji dla wybranych kodów przedstawiono na rys.3.3. Pełna klasyfikacja znajduje się w tabeli 3.1.

W powyższym schemacie przedstawiono typy interfejsu tylko pomiędzy fazą stałą a ciekłą. W przypadku granicy pomiędzy dwoma fazami stałymi (np. austenit – grafit) schemat klasyfikacji jest ten sam. Na schematach interfejsu przedstawionych w niniejszej pracy zastosowano następujące zasady: komórki interfejsu ciekło-stałego oddzielane są od komórek ciekłych podwójną cienką linią, natomiast od fazy stałej, pojedynczą grubą linią. W przypadku interfejsu zawierającego dwie fazy stałe, jest on oddzielony od nich po obu stronach grubą linią. Przerywaną linią zaznaczono możliwą granicę rzeczywistego frontu,

32

4

64

2

Int

8

16

1

128

gdyż jego faktyczne położenie w modelowaniu nie jest określane. Znany jest jedynie udział fasy w danej komórce.

a) b) c) d) e) f) g)

Rys. 3.3. Schemat klasyfikacji typów interfejsu na podstawie ich kodów: a) granica płaska,

b) narożnik, c) kąt, d) kanalik, e) dno kanalika, f) wnęka zamknięta, g) złożone (linia

pogrubiona – granica pomiędzy komórkami fazy stałej a komórkami interfejsu; linia

podwójna – granica pomiędzy komórkami fazy ciekłej a komórkami interfejsu; linia

przerywana – prawdopodobna rzeczywista pozycja frontu rosnącego ziarna; S – komórka

w stanie stałym, L – komórka w stanie ciekłym, I – komórka interfejsu)

Istnieją prace np. [27-30], w których stosuje się tzw. metodę śledzenia frontu (ang. front tracking), która wyznacza położenie frontu za pomocą punktów kontrolnych bezpo-średnio na granicach komórek, łącząc te punkty odcinkami linii prostych.

Tab. 3.1. Klasyfikacja typów interfejsu na podstawie ich kodów

Kod interfejsu na podstawie maski z rysunku 3.2

T ypy i nt er fe js u: zarodek 0 granice pła-skie 1, 2, 4, 8, 17, 18, 34, 36, 50, 68, 72, 100, 129, 136, 145, 200 narożniki 16, 32, 64, 128 kąty 3, 6, 9, 12, 19, 22, 25, 35, 38, 44, 51, 54, 70, 73, 76, 86, 89, 102, 108, 118, 131, 137, 140, 147, 153, 163, 172, 179, 201, 204, 217, 236 kanaliki 5, 10, 21, 26, 37, 42, 53, 58, 68, 74, 85, 90, 101, 106, 117, 122, 133, 138, 149, 154, 165, 170, 181, 186, 197, 202, 213, 218, 229, 234, 245, 250 dna kanali-ków 7, 11, 13, 14, 23, 27, 29, 30, 39, 43, 45, 46, 55, 59, 61, 62, 71, 75, 77, 78, 87, 91, 93, 94, 103, 107, 109, 110, 119, 123, 125, 126, 135, 139, 141, 142, 151, 155, 157, 158, 167, 171, 173, 174, 183, 187, 189, 190, 199, 203, 205, 206, 215, 219, 221, 222, 231, 235, 237, 238, 247, 251, 253, 254 wnęki za-mknięte 15, 31, 47, 63, 79, 95, 111, 127, 143, 159, 175, 191, 207, 223, 239, 255 złożone 20, 24, 28, 33, 40, 41, 48, 49, 52, 56, 57, 60, 65, 66, 67, 80, 81, 82, 83, 84, 88, 92, 96, 97, 98, 99, 104, 105, 112, 113, 114, 115, 116, 120, 121, 124, 130, 132, 134, 144, 146, 148, 150, 152, 156, 160, 161, 162, 164, 166, 168, 169, 176, 177, 178, 180, 182, 184, 185, 188, 192, 193, 194, 195, 196, 198, 208, 209, 210, 211, 212, 214, 216, 220, 224, 225, 226, 227, 228, 230, 232, 233, 240, 241, 242, 243, 244, 246, 248, 249, 252

3.3. Transport masy

Równanie różnicowe opisujące dyfuzję składnika wyprowadzono korzystając z bilansu masy. Dla przykładowej komórki będącej w stanie 3 tj. zawierającej austenit γ i ciecz L, tak jak to pokazano na rys. 3.4, równanie bilansu masy, w przypadku gdy przy-rost udziału fazy stałej w jednym kroku czasowym ∆τ wyniesie ∆fγ, ma następującą postać

(

)

(

)

(

)

ij N E S W t j i L L j i t j i L j i j i t j i L L j i j i t j i L f kC f C f f kC f f C Φ + Φ + Φ + Φ = ρ − ρ − − − ρ ∆ + + ρ ∆ − − γ γ γ γ γ γ + γ γ + , , , , , , , , , , , , 1 , , , , , , 1 , , 1 1 (3.2) gdzie:

CL – stężenie węgla w cieczy wyrażone w ułamkach wagowych, fγ – udział austenitu w komórce,

k – równowagowy współczynnik rozdziału węgla pomiędzy austenitem a cieczą, ρ , ρ gęstość odpowiednio cieczy i austenitu,

t – górny indeks oznaczający numer kroku czasowego,

Φ – zmiana zawartości na skutek strumienia dyfuzyjnego od strony sąsiedniej komórki.

Rys. 3.4. Przykładowy interfejs pomiędzy cieczą a austenitem: kolor szary – austenit; kolor biały -

ciecz

Lewa strona równania (3.2) jest równaniem bilansu masy węgla w komórce, nato-miast prawa strona (tzw. człony dyfuzyjne Φ) odpowiada za zmianę stężenia węgla w ko-mórce na skutek strumieni dyfuzyjnych w otoczeniu von Neumana, czyli w kierunkach N, E, S, i W. Równania różnicowe dla obliczenia wartości Φ mogą mieć różne postacie w zależności od stanów sąsiednich komórek i konfiguracji interfejsu (kodu interfejsu). Przykładowo, węgiel z komórki centralnej i,j na rys. 3.4 dyfunduje do sąsiedniej komórki w kierunku E poprzez fazę ciekłą, a w kierunku W przez fazę stałą (austenit), natomiast w kierunkach N i S zarówno przez ciekłą jak i stałą. W kierunku N powierzchnia, przez którą dyfunduje węgiel w austenicie, wyliczana jest jako średnia arytmetyczna w z udziału faz fγ,i,j i fγ,N. Podobne wartości średnie pola powierzchni wyliczane są dla danego kierun-ku, jeśli sąsiednia komórka jest interfejsem. Uwzględniając to, prawa strona równania (3.2), dla przypadku pokazanego na rys. 3.4, przyjmie następującą postać

− w kierunku N

(

)

(

)

(

)

(

, ,, 1 , ,,

)

2 a C C D w C C k D w L t j i L t N L L N t j i L t N L N N τ ∆ ρ − − + ρ − = Φ _{γ γ} (3.3) − w kierunku E

(

)

L t j i L t E L L E C C a D ∆τ − ρ = Φ 2 , ,, (3.4) − w kierunku S

(

)

(

)

(

)

(

, ,, 1 , ,,

)

2 a C C D w C C k D w t _L j i L t S L L S t j i L t S L S S τ ∆ ρ − − + ρ − = Φ _{γ γ} (3.5)

− oraz dla kierunku W

(

γ

)

γ γ − ρ τ ∆ = Φ t j i L t W W C kC a D _{2 , ,,} (3.6) gdzie:

a – długość boku kwadratowej komórki AK,

DL, Dγ – współczynnik dyfuzji węgla odpowiednio w cieczy i austenicie, Cγ – stężenie węgla w austenicie wyrażone w ułamkach wagowych, ∆τ – krok czasowy, 2 , , ,N ij N f f w = γ + γ 2 , , ,S i j S f f w = γ + γ

Z równania (3.2) można określić nowe stężenie węgla w następnym kroku czasowym t+1

(

)

(

)

(

γ γ

)

(

γ γ

)

γ γ γ γ + ρ ∆ + + ρ ∆ − − ρ + ρ − + Φ + Φ + Φ + Φ = j i j i L j i j i j i L j i t j i L W S E N t j i L f f k f f f k f C C , , , , , , , , , , , , , , 1 , , 1 1 (3.7)

Biorąc pod uwagę że w graficie nie występuje dyfuzja węgla, a jego stężenie wynosi C = 1, można analogicznie wyprowadzić równanie dla interfejsu grafit-ciecz. Dla konfigu-racji z rys. 3.4 wystarczy zastąpić austenit grafitem, aby otrzymać równanie określające nowe stężenie węgla w komórce interfejsu grafit-ciecz (stan 4) dla następnego kroku cza-sowego

(

)

(

)

(

grij

)

gri j gr L t j i L j i gr j i gr N S t j i L t E L t S L S t N L N L t j i L f f f f C f f w w C C C w C w a D C , , , , , , , , , , , , , , , 2 1 , , 1 1 1 1 ∆ − − ρ ρ ∆ − − + + ∆ − − + + − + + τ ∆ = + (3.8)

Dla komórek dwufazowych zawierających austenit i ciecz równanie (3.7) pozwala obliczyć stężenie węgla w cieczy. Stężenie węgla w austenicie w kontakcie z cieczą wy-znacza się na podstawie równowagowego współczynnika rozdziału k. Dla komórek zawie-rających grafit i ciecz, wynik obliczenia z równania (3.8) daje stężenie węgla w cieczy.

Równanie opisujące dyfuzję węgla dla komórek jednofazowych (oprócz stanu 2, tj. grafitu) przyjmuje prostszą postać. Można je otrzymać bezpośrednio z drugiego prawa Ficka, zamieniając pochodne cząstkowe na różnice skończone. Dla przykładu zilustrowa-nego na rys. 3.5, w przypadku schematu jawzilustrowa-nego otrzymano równanie różnicowe (3.9)

(

)

t j i fa fa t j i fa t W fa t S fa t E fa t N fa t j i fa C a D C C C C C C +1_,, = _, + _, + _, + _, −4 _,, ∆₂τ+ _,, (3.9)

gdzie indeks dolny fa oznacza jedną z faz tj. austenit albo ciecz.

Przedstawione równania (3.7) i (3.8) otrzymano dla szczególnych przypadków zilu-strowanego interfejsu. Dla innych konfiguracji w równaniach tych będą zmieniać się jedy-nie człony gradientowe, natomiast lewa strona równań dla danej komórki interfejsu pozo-staje niezmienna niezależnie od tego, po której stronie znajdują się sąsiednie komórki in-terfejsu tego samego typu.

Rys. 3.5. Przykład komórki jednofazowej (współrzędne i,j) wraz z komórkami sąsiednimi, służący

do wyjaśnienia równania (3.9) (fa jedna z dwóch faz – ciecz albo austenit; Int jest komórką interfejsu o stanie 3, 4 albo 5)

3.4. Przepływ ciepła

3.4.1. Zmiana pola temperatury

Przepływ ciepła w układzie opisany jest znanym równaniem Fouriera-Kirchoffa z uwzględnieniem ciepła przemiany fazowej. Jest to paraboliczne równania różniczkowe z pochodnymi cząstkowymi, które dla przestrzeni dwuwymiarowej ma postać

T v q y T y x T x T c _+      ∂ ∂ λ ∂ ∂ +       ∂ ∂ λ ∂ ∂ = τ ∂ ∂ (3.10) gdzie: T – temperatura, τ – czas,

λ – współczynnik przewodzenia ciepła,

cv – ciepło właściwe przypadające na jednostkę objętości, x, y – współrzędne,

qT – wydajność źródła ciepła.

Człon źródłowy qT powoduje nieliniowość powyższego równania. Jest on różny od zera jedynie w komórkach interfejsu. W innych komórkach (tj. jednofazowych) przyjmuje się wartość zerową.

3.4.2.

Rozwiązanie numeryczne

Rozwiązanie numeryczne równania (3.10), wraz z uwzględnieniem warunków po-czątkowych i brzegowych, nie nastręcza problemów. Sposoby otrzymywania układów równań modelowych są znane od lat i dobrze opisane w literaturze np. [31-33]. W niniejszej pracy zastosowano schemat różnicowy jawny. Stosowanie takiego schematu dla dwuwymiarowej siatki z kwadratowymi komórkami narzuca z góry maksymalny do-puszczalny krok czasowy ∆τmax,T [34-36] dla pola temperatury, przy którym rozwiązanie jest stabilne pod względem numerycznym. Wielkość dopuszczalnego kroku czasowego zależy od wartości współczynnika przewodzenia ciepła, ciepła właściwego oraz od wielko-ści kroku siatki, na której wykonywane są obliczenia. Dla pola bezźródłowego krok ten jest jednoznacznie określony i dany jest wzorem [37]

λ ∆ = τ ∆ d c x_{T v} T 2 2 max, (3.11)

gdzie:

d – liczba wymiarów przestrzennych.

W przypadku pól źródłowych (czyli z uwzględnieniem źródła ciepła przemiany fa-zowej zależnej od temperatury, a taka sytuacja ma miejsce w komórkach interfejsu) po-wyższe ograniczenie powinno być bardziej ostrzejszym. Przykładowo w [34] dla modelo-wania w przestrzeni dwuwymiarowej ograniczenie na maksymalny krok czasowy ma po-stać λ ∆ = τ ∆ 5 . 4 2 max, v T T c x (3.12)

Niestety nie jest możliwe jednoznaczne określenie wartości Δτmax,T [38

35

]. Zazwyczaj dopuszczalny krok czasowy dobiera się w sposób doświadczalny, przyjmując jego stałą wartość w obliczeniach, tak jak to ma miejsce w prezentowanej pracy. Może on być rów-nież zmienny w trakcie obliczeń, lecz wymaga to sprawdzania dodatkowych kryteriów stabilności [ , 38, 39]. Kryteria te mogą być bezpośrednio związane z maksymalnym do-puszczalnym przyrostem udziału fazy stałej w komórce w jednym kroku czasowym, tak jak to wyznaczono doświadczalnie w [40

Przekształcenie równania

].

(3.10) do jawnego schematu różnicowego przyjmuje po-stać

(

)

źrir jr t jr ir T t jr ir jr ir T t W W T t S S T t E E T t N N T t jr ir T T x T a T a T a T a T a T _,1 _{, , , ,} 4 _{, , , 2} + _, +∆ _{, ,} ∆ τ ∆ − + + + = + (3.13)

gdzie a_T jest średnim współczynnikiem wyrównywania temperatury określonym dla

każ-dej pary sąsiednich komórek N, E, S, W z komórką ir,jr. Człon ∆Tźr,ir,jr odpowiada za zmianę temperatury spowodowaną wydzielaniem się ciepła przemiany fazowej. Wartość

T

a wyliczana jest na podstawie średniej harmonicznej średnich współczynników przewo-dzenia ciepła dla komórki centralnej i sąsiedniej, oraz na podstawie średniego ciepła wła-ściwego tej komórki

(

)

gdzie

{

N E S W

}

c a jr ir v jr ir jr ir T , , , , 2 , , , , , _{λ +λ} = λ λ =     (3.14) gdzie:

v

c – średnie ciepło właściwe.

Człon 4aT,ir,jr z równania (3.13) równy jest

W T S T E T N T jr ir T a a a a a , , , , , , 4 = + + + (3.15)

3.4.3.

Dopuszczalny krok czasowy rozwiązania numerycznego

Podobne do (3.12) ograniczenie dopuszczalnego kroku w przypadku pola dyfuzji ma postać [34] 1 2 max, 5 . 4 D a C = τ ∆ (3.16) gdzie:

D1 – współczynnik dyfuzji pierwiastka w odpowiedniej fazie.

W przypadku zastosowania siatek z jednakowym krokiem przestrzennym dla pola temperatury i dyfuzji oraz dla typowych parametrów materiałowych jak dla stopu Fe-C, dopuszczalny krok czasowy schematu jawnego rozwiązania numerycznego, wyliczony z zależności (3.12), jest 103÷104 razy mniejszy niż dla pola stężenia zgodnie z ograniczeniem (3.16).

W niektórych pracach stosowano inne metody rozwiązywania równań różnicowych, np. w [40] zastosowano metodę iteracji prostej, a w [41

39 ] metodę kierunków naprzemien-nych. Podejście takie pozwala na dobieranie większych kroków czasowych, niżeli w przypadku stosowania metod jawnych. Innym sposobem zwiększenia kroku czasowego dla schematu jawnego jest zastosowanie podwójnej dyskretyzacji obszaru [ ]. W niniejszym modelu zdecydowano się na zastosowanie tego podejścia.

Idea jest następująca. Na sieć pola stężenia i AK o wspólnym kroku przestrzennym a (nazywana dalej siecią gęstą) nakłada się sieć pola temperatury (zwana dalej siecią rzad-ką) o kroku będącym wielokrotnością nk kroku sieci gęstej. Wartości temperatury dla sieci rzadkiej, wyliczane za pomocą równania (3.13), przypisane są do jej węzłów (czerwone punkty na rys. 3.6). Stężenie wyliczane za pomocą odpowiednich równań przedstawionych w rozdziale 3.3 przypisywane jest do środka komórek sieci gęstej.

Oprócz zwiększenia dopuszczalnego kroku czasowego, zastosowanie podwójnej dyskretyzacji obszaru można uzasadnić tym, iż wartość temperatury w mikroobszarze

zmienia się nieznacznie w porównaniu ze zmianami składu chemicznego. Często w publikacjach rozważa się wzrost ziaren w warunkach jednorodnego pola temperatury, np. w [42 - 44].

Rys. 3.6. Przykład węzła sieci rzadkiej wraz ze stosowanym otoczeniem, w którym oblicza się

średnie parametry materiałowe (sieć pola temperatury – czerwone linie; siatka AK i pola stężenia – czarne linie)

Wyliczenie średnich parametrów materiałowych, występujących w zależności (3.14), przebiega następująco. Dla każdego węzła sieci rzadkiej wyznacza się jego objętość kon-trolną będącą kwadratem o boku nk⋅a (patrz zielony przerywany prostokąt na rys. 3.6). Na-stępnie wewnątrz tego kwadratu wylicza się średnią ważoną danego parametru materiało-wego. Wagami w tym przypadku są udziały poszczególnych faz, przy czym za jednostko-wy udział przyjmuje się jedną komórkę AK. Suma udziału wszystkich faz wewnątrz oto-czenia równa jest nk2. Poszczególne średnie wartości parametrów materiałowych dla węzła sieci rzadkiej wylicza się z następującej formuły

2 k L L gr gr n f f f

∑

∑

∑

+ϑ +ϑ ϑ = ϑ γ γ (3.17)

gdzie:

ϑ – odpowiedni parametr materiałowy tj. ciepło właściwe lub współczynnik przewodzenia ciepła.

Znaki sum w powyższym równaniu oznaczają sumowanie udziałów poszczególnych faz we wszystkich komórkach AK znajdujących się wewnątrz objętości kontrolnej.

W podobny sposób wylicza się człon ∆Tźr,ir,jr z równania (3.13)

2 , , k v gr gr gr L gr L L L jr ir źr n c f H f H f H T = ∆ →γ

∑

∆ →γ +∆ →

∑

∆ → +∆ γ→

∑

∆ γ→ ∆ (3.18) gdzie:

∆H – entalpia przemiany fazowej,

∆f – przyrost udziału fazy wyliczany przy pomocy równania (3.24).

Indeksy dolne w powyższym równaniu oznaczają następujące przemiany fazowe: L→γ – wydzielanie się austenitu z cieczy,

L→gr – wydzielanie się grafitu z cieczy, γ→gr – wydzielanie się grafitu z austenitu.

W proponowanym modelu nie uwzględnia się przemian odwrotnych do wyżej wy-mienionych tzn. nadtapiania się austenitu lub rozpuszczania się grafitu w cieczy lub auste-nicie. Przykład modelowania, które uwzględnia nadtapianie się gałęzi dendrytu austenitu na skutek zmiany stężenia węgla i wydzielania się ciepła krystalizacji, można zobaczyć w pracy [45

3.4.4. Interpolacja temperatury

].

Temperatura w obszarze obliczeniowym wyliczana z równania (3.13) przypisywana jest tylko do węzłów sieci rzadkiej (rys. 3.6). W celu określenia temperatury dla komórki sieci gęstej stosuje się interpolację wartości temperatury z węzłów sieci rzadkiej. Schemat interpolacji jest następujący. Dla każdej komórki sieci gęstej wyznacza się cztery najbliż-sze węzły sieci rzadkiej tak, aby dana komórka leżała wewnątrz kwadratu zbudowanego z tych węzłów. Następnie określa się współrzędne (x, y) komórki w lokalnym układzie współrzędnych (rys. 3.7). Jako długość jednostkową w tym układzie przyjęto długość kro-ku ∆xT sieci rzadkiej. Węzły sieci rzadkiej o temperaturze T0,0, T1,0, T0,1, T1,1, mają współ-rzędne odpowiednio (0,0), (1,0), (0,1), (1,1), natomiast x i y ∈(0;1).

Rys. 3.7. Schemat interpolacji temperatury dla sieci gęstej

Jako bazowe funkcje interpolacji ξ wybrano [39]

       = ξ − = ξ − = ξ + − − = ξ xy xy y xy x y x xy 1 , 1 1 , 0 0 , 1 0 , 0 1 (3.19)

Wartość temperatury T(x,y) dla komórki sieci rzadkiej oblicza się z następującej formuły interpolacyjnej

( )

x,y 0,0T0,0 1,0T1,0 0,1T0,1 1,1T1,1

T =ξ +ξ +ξ +ξ (3.20)

3.5. Szybkość przemiany

3.5.1.

Znane rozwiązania

Znane są różne sposoby wyznaczania zmiany udziału fazy stałej w kwadratowej ko-mórce interfejsu w jednym kroku czasowym. Na przykład w pracy [46] zastosowano za-leżność

(

)

(

p

)

a u u u u a f _{x y x y} ⋅ 1+η1−2      _{+ −} ∆τ τ ∆ = ∆ (3.21)

gdzie:

ux, uy – rzuty wektora prędkości na osie układu współrzędnych, Δτ – krok czasu modelowania,

η – amplituda czynnika losowego (rzędu 0.1), p – liczba losowa z zakresu (0; 1).

Czynnik (1 + η(1–2p)) w powyższym równaniu uwzględnia fluktuacje prędkości wzrostu.

Inny sposób wyznaczania zmiany udziału fazy stałej zaproponowano w [34]. Wpro-wadzono tam tzw. długość wzrostu L. Jest to długość odcinka, który łączy środek bieżącej komórki ze środkiem sąsiedniej komórki w kierunku aktualnego wzrostu. Przyrost fazy liczony tą metodą realizowany jest za pomocą poniższego równania

1 1 sin cosθ + θ ⋅ τ ∆ = ∆ u L f (3.22) gdzie:

θ1 – kąt pomiędzy kierunkiem krystalograficznym najszybszego wzrostu ziarna a wybraną linią L (łączącą środek komórki bieżącej ze środkiem sąsiedniej). Parametr L może przyjmować wartość L = a dla wzrostu w kierunku najbliższych są-siednich komórek lub L=a 2 dla wzrostu w kierunku komórek po przekątnej.

Nieco inną metodę oceny przyrostu Δf w jednym kroku czasowym zaprezentowano w pracy [40]. Metoda ta może być stosowana również dla komórek prostokątnych. Zapro-ponowane równanie ma postać

y y x xu a u a u f + τ ∆ = ∆ 2 (3.23) gdzie:

ax, ay – długości odpowiednich boków prostokątnej komórki AK.

W przypadku, gdy siatka automatu komórkowego ma kwadratowe komórki, równa-nie (3.23) przekształca się do postaci

(

θ+ θ

)

τ ∆ ⋅ = ∆ sin cos a u f (3.24)

gdzie:

θ – kąt pomiędzy wektorem normalnym na froncie krystalizacji a osią X układu współrzędnych.

Niezależnie od wybranego sposobu obliczania Δf, dla każdej komórki interfejsu przyjmuje się początkową wartość f0 = 0 i wykonuje się całkowanie numeryczne po czasie do momentu uzyskania wartości końcowej fk >= 1.

3.5.2.

Testowanie równań

W pracy [47

Kierunek wektora prędkości w komórkach interfejsu wyznacza się wyłącznie na podstawie współrzędnych bieżącej komórki frontu, dla której aktualnie liczony jest przy-rost fazy stałej

] zaproponowano tzw. test "sztywnego koła" służący do testowania po-prawności równań do wyznaczania szybkości przemiany w modelach opartych na wyko-rzystaniu automatów komórkowych. Opracowane zadanie testowe powinno zapewnić wzrost ziarna o kształcie okrągłym. W tym celu modeluje się wzrost okręgu z zarodka, dla którego początkowa wartość udziału fazy stałej wynosi f0 = 1. Zarodek umieszczono w punkcie środkowym (o współrzędnych i0,j0) siatki z kwadratowymi komórkami.

0 0 arctan j j i i − − = θ (3.25) gdzie:

i, j – numery komórki frontu na siatce wzdłuż osi współrzędnych, i0, j0 – numery komórki zarodka.

Składowe wektora prędkości równe są

θ = θ

= cos , sin

u_x u u_y u (3.26)

Schemat ideowy testu sztywnego koła przedstawiono na rys. 3.8. Prędkość migracji frontu założono na podstawie równania

τ ∆ =h a

u (3.27)

Wyniki symulacji zależą od bezwymiarowego kryterium h, które równe jest stosun-kowi przemieszczenia się frontu w ciągu jednego kroku czasowego do długości boku

ko-mórki interfejsu. Wybrano je na poziomie 0.01, 0.02, 0.05 i 0.1. Oznacza to, że w komórkach interfejsu o współrzędnych (i,j0) i (i0,j) przemiana trwa dokładnie 1/h kro-ków czasowych. Dla pozostałych komórek (gdy i ≠ i0˄ j ≠ j0) w ostatnim kroku czasowym przemiany udział fazy stałej liczony w taki sposób może przekroczyć 1. W takich przypad-kach w zadaniu testowym wartość f jest obcinana do 1.

Rys. 3.8. Rysunek ilustrujący test sztywnego koła

Ilość kroków czasowych podczas przeprowadzania testów wynosiła 100/h, co ozna-czało nominalne zwiększenie promienia cząstki o 100⋅a. Wartość testową, czyli bezwymia-rowy promień R ziarna w każdej komórce interfejsu po zakończeniu testu określano na podstawie równania

(

₀

) (

2 ₀

)

2 _, 0.5

,j = u − + u − + i j −

i i i j j f

R (3.28)

Wyniki testów zaprezentowano w tabeli 3.2. Pod ocenę wzięto dwie wartości: średni bezwymiarowy promień oraz średnie względne pole powierzchni cząstki. Pierwsza z nich jest to stosunek średniego promienia do promienia nominalnego Rn. Średni promień ziarna wyliczono jako średnią arytmetyczną z wszystkich Ri,j komórek interfejsu w momencie zakończenia testu. Wartość końcowa promienia nominalnego dla przyjętych warunków zakończenia we wszystkich testach wynosiła Rn = 100.5. Druga wartość oceniania jest równa stosunkowi sumy wszystkich wartości fi,jdo pola powierzchni koła o promieniu Rn.

Tab. 3.2. Wyniki testu sztywnego koła testowanych równań na przyrost udziału fazy stałej

Równanie h

Względne Odchylenie promienia pole po-wierzchni ziarna średni promień (R Rminn) – (Rmax – Rn) standardowe (3.21) 0.01 2.44 1.581 31 99 19.9 0.02 2.39 1.561 29 97 19.5 0.05 2.21 1.506 25 89 18.6 0.10 1.96 1.252 19 33 3.97 (3.22) 0.01 0.923 0.959 –7.4 –0.02 2.43 0.02 0.919 0.957 –7.7 –0.02 2.50 0.05 0.907 0.951 –8.1 –0.05 2.55 0.10 0.823 0.907 –10.8 –0.68 1.19 (3.24) 0.01 0.996 0.998 –0.51 0 0.13 0.02 0.993 0.996 –0.88 0 0.21 0.05 0.981 0.990 –2.53 0 0.63 0.10 0.873 0.936 –9.00 –3.51 1.58

Na rys. 3.9 pokazano kształty ziaren uzyskane podczas testowania równań. Ze względu na czterokrotną symetrię ograniczono się do pokazania tylko 1/4 powierzchni. Widać, że przebieg granic uzyskany za pomocą równań (3.21) i (3.22) znacznie odbiega od założonego okrągłego kształtu. Jako miarę anizotropii odchylenia komórek interfejsu obli-czono różnicę pomiędzy minimalnym a nominalnym promieniem, maksymalnym a nominalnym promieniem oraz odchylenie standardowe promieni na obwodzie. Wyniki w tabeli 3.2 pokazują, że stosowanie równania (3.24) zapewnia również minimalną anizo-tropię wyników modelowania.

a) b) c)

Rys. 3.9. Kształty ziaren uzyskane w teście sztywnego koła: 1/4 powierzchni:

3.6. Kierunek normalny do frontu

3.6.1. Znane metody

W celu określenia kierunku interfejsu w pracy [46] użyto tzw. F-wektor. Wartość te-go wektora równa jest całkowitej ilości udziału fazy stałej f wewnątrz okręgu o promieniu R. Wektor ten zorientowany jest wzdłuż linii łączącej środek masy obszaru fazy stałej, znajdującej się w kole o promieniu R, ze środkiem analizowanej komórki. Jednostkowy wektor n może być wyliczony z poniższej zależności

(

) (

)

[

_{2 2}

]

12 m m m m m m m m m f g y f g x f g m m m Σ + Σ Σ − = r n (3.29) gdzie:

rm – wektor łączący środek komórki, dla której liczony jest kierunek normalny, ze środkiem komórki m,

fm – udział fazy stałej w komórce m, xm, ym – współrzędna wektora rm,

gm – współczynnik określający, jaka część komórki znajduje się wewnątrz koła o promieniu R.

Sumowanie wykonuje się po wszystkich komórkach m w otoczeniu R.

Kąt, jaki tworzy wektor n o składowych (nx; ny) z osią X, wyliczany jest z prostej za-leżności trygonometrycznej x y n n arctan = θ (3.30)

Inną metodę wyznaczania kierunku frontu, zaproponowano w pracy [36]. Metoda ta bazuje na koncepcji rozmytej granicy międzyfazowej. Udział fazy stałej f w komórce inter-fejsu oraz w komórkach z jej otoczenia von Neumana (rys. 3.1) może przyjmować warto-ści od 0 do 1. Traktując rozkład udziału fazy stałej na siatce AK, jako funkcję współrzęd-nych x i y można wyliczyć gradient pola zmiennej f. Jednostkowy wektor prostopadły do frontu przybliża się jako

f f ∇ ∇ − = n (3.31)

Analogicznie jak we wzorze (3.30), kąt jaki tworzy wektor normalny do frontu moż-na wyliczyć jako x f y f ∂ ∂ ∂ ∂ = θ arctan (3.32)

Schemat wyjaśniający sposób wykorzystania równań (3.30) i (3.32) przedstawiony jest na rys. 3.10. W środku komórki (zaznaczonej czerwoną przerywaną linią), dla której wyznaczany jest kąt θ, zaczepiono dwa lokalne pokrywające się układy współrzędnych. Na pierwszym z nich, o osiach X i Y, zaznacza się składowe wektora n, które mają charakter ciągły. Drugi układ o osiach i i j służy do określania współrzędnych sąsiednich komórek AK. Współrzędne te są liczbami całkowitymi. Udział fazy stałej f w danej komórce ma charakter ciągły i może przyjmować wartości od 0 do 1.

Rys. 3.10. Schemat ilustrujący sposób wyznaczania kierunku normalnego do frontu

z wykorzystaniem równań (3.29), (3.32) i (3.33) (oznaczenie granic zgodnie z rysunkiem 3.3; żółta

przerywana linia – przykładowe otoczenie stosowane w teście miękkiego koła)

Określenie kierunku normalnego przy użyciu F-wektora wymaga określenia współ-czynników gm ze wzoru (3.29). Dla komórek, które leżą wewnątrz koła i nie są przecięte okręgiem o promieniu R, wartość ta wynosi 1. Dla wszystkich komórek przeciętych

okrę-giem, współczynniki przyjmują natomiast wartości równe ułamkowi części pola komórki AK, jaka znajduje się wewnątrz tego okręgu. Dokładne określenie w takim przypadku war-tości poszczególnych współczynników wymaga na przykład całkowania numerycznego, co jest operacją dość czasochłonną. W niniejszej pracy zmodyfikowano metodę wykorzysta-nia F-wektora, uwzględniając w obliczeniach tylko te komórki, których środek znajdował się wewnątrz lub na brzegu analizowanego otoczenia, przyjmując dla nich wartość gm = 1. Po uwzględnieniu tej modyfikacji i wyliczeniu składowych wektora n (z równania (3.29)) oraz podstawieniu ich do wzoru (3.30) otrzymano prostą zależność dla obliczenia kąta θ

∑

∑

= θ m m m m m m f i f j arctan (3.33)

3.6.2. Testowanie metod

Zadanie testowe polegało na umieszczeniu zarodka na środku kwadratowego obszaru AK z kwadratowymi komórkami. Początkowy udział fazy stałej w komórce centralnej wy-nosił f0 = 1, natomiast w komórkach jej otoczenia Moore'a udziały fazy stałej wybrano tak, aby odpowiadały cząstce okrągłej o promieniu a. Do wyznaczania przyrostu fazy stałej w jednym kroku czasowym wykorzystano równanie (3.24), zapewniające minimalny po-ziom anizotropii. Bezwymiarowe kryterium h kroku czasowego przyjęto na poziomie 0.01. Zamiast równania (3.25) służącego do obliczenia kąta θ, jaki tworzy wektor kierunkowy frontu, wykorzystano testowane równania (3.32) i (3.33). Idealna metoda służąca do wy-znaczania kierunku normalnego do frontu, pomijając wpływ błędów z metody (3.24), po-winna zapewnić wzrost ziarna o okrągłym kształcie. Test ten opisany w [48

Testowane równania były sprawdzane dla kilku różnych promieni otoczenia. Ko-mórka brana była pod uwagę tylko wtedy, jeśli jej środek znajdował się wewnątrz lub na brzegu okręgu. W zależności od promienia ilość m uwzględnianych komórek wynosiła: 5, 9, 13, 21 i 25. Implementacja sprawdzanych równań przedstawiona jest w tabeli

] został nazwa-ny testem miękkiego koła.

3.3. Gdy promień otoczenia równy jest a lub a 2_{, implementacja obu testowanych równań}

wyglą-da tak samo.

Niedokładność wyznaczenia kierunku normalnego do interfejsu za pomocą powyż-szych równań prowadzi do rozbieżności pomiędzy uzyskanymi kształtami ziaren a kształtem uzyskanym w teście sztywnego koła, przedstawionym na rysunku 3.9 c.

Wyni-ki testu miękkiego koła, gdy promień nominalny wynosił Rn = 101a, zaprezentowano na rysunku 3.11. Wizualnie minimalne zniekształcenie daje równanie f pokazane w tabeli 3.3.

Tab. 3.3. Implementacja testowanych równań na kierunek interfejsu w zależności od promienia

otoczenia Promień (Ilość komórek m)

Rów-nanie Argument funkcji arctan Indeks

a (5) (3.32) (3.33) _{1,0 1,0} 1 , 0 1 , 0 f f f f − − − − a 2 a (9) (3.32) (3.33) _{1,1 1,0 1, 1 1,1 1,0 1, 1} 1 , 1 1 , 0 1 , 1 1 , 1 1 , 0 1 , 1 − − − − − − − − − − − − − + + − − − + + f f f f f f f f f f f f b a 2 (13) (3.32) 0 , 2 1 , 1 0 , 1 1 , 1 1 , 1 0 , 1 1 , 1 0 , 2 2 , 0 1 , 1 1 , 0 1 , 1 1 , 1 1 , 0 1 , 1 2 , 0 − − − − − − − − − − − − − − − − + + + − − − − + + + f f f f f f f f f f f f f f f f c (3.33)

(

(

)

)

1 , 1 0 , 1 1 , 1 1 , 1 0 , 1 1 , 1 0 , 2 0 , 2 1 , 1 1 , 0 1 , 1 1 , 1 1 , 0 1 , 1 2 , 0 2 , 0 2 2 − − − − − − − − − − − − − − − + + + − − − − + + + − f f f f f f f f f f f f f f f f d

5

a

(21) (3.32)         − − − − − + + + + + + − − − + +         − − − − − + + + + + + − − − + + − − − − − − − − − − − − − − − − − − − − − − − − − − − − 2 , 1 1 , 1 0 , 1 1 , 1 2 , 1 2 , 1 1 , 1 0 , 1 1 , 1 2 , 1 1 , 2 0 , 2 1 , 2 1 , 2 0 , 2 1 , 2 1 , 2 1 , 1 1 , 0 1 , 1 1 , 2 1 , 2 1 , 1 1 , 0 1 , 1 1 , 2 2 , 1 2 , 0 2 , 1 2 , 1 2 , 0 2 , 1 f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f e (3.33)

(

)

(

)

      − − − − − + + + + + + − − − + +       − − − − − + + + + + + − − − + + − − − − − − − − − − − − − − − − − − − − − − − − − − − − 2 , 1 1 , 1 0 , 1 1 , 1 2 , 1 2 , 1 1 , 1 0 , 1 1 , 1 2 , 1 1 , 2 0 , 2 1 , 2 1 , 2 0 , 2 1 , 2 1 , 2 1 , 1 1 , 0 1 , 1 1 , 2 1 , 2 1 , 1 1 , 0 1 , 1 1 , 2 2 , 1 2 , 0 2 , 1 2 , 1 2 , 0 2 , 1 2 2 f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f 2 2a (25) (3.32)         − − − − − + + + + + + − − − − − + + + +         − − − − − + + + + + + − − − − − + + + + − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − 2 , 1 1 , 1 0 , 1 1 , 1 2 , 1 2 , 1 1 , 1 0 , 1 1 , 1 2 , 1 2 , 2 1 , 2 0 , 2 1 , 2 2 , 2 2 , 2 1 , 2 0 , 2 1 , 2 2 , 2 1 , 2 1 , 1 1 , 0 1 , 1 1 , 2 1 , 2 1 , 1 1 , 0 1 , 1 1 , 2 2 , 2 2 , 1 2 , 0 2 , 1 2 , 2 2 , 2 2 , 1 2 , 0 2 , 1 2 , 2 f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f g (3.33)

(

)

(

)

      − − − − − + + + + + + − − − − − + + + +       − − − − − + + + + + + − − − − − + + + + − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − − 2 , 1 1 , 1 0 , 1 1 , 1 2 , 1 2 , 1 1 , 1 0 , 1 1 , 1 2 , 1 2 , 2 1 , 2 0 , 2 1 , 2 2 , 2 2 , 2 1 , 2 0 , 2 1 , 2 2 , 2 1 , 2 1 , 1 1 , 0 1 , 1 1 , 2 1 , 2 1 , 1 1 , 0 1 , 1 1 , 2 2 , 2 2 , 1 2 , 0 2 , 1 2 , 2 2 , 2 2 , 1 2 , 0 2 , 1 2 , 2 2 2 f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f f h

W celu ilościowej oceny testowanych równań, po zakończeniu testu dla każdej ko-mórki interfejsu na podstawie równania (3.28) obliczono promień, a następnie określono jego względną maksymalną, minimalną oraz średnią wartość. Wyliczono też jego względ-ne odchylenie standardowe odniesione do długości boku komórki. Jako ocenę ilościową posłużono się również bezwzględną wartością różnicy pomiędzy kątem θ, wyliczanym za pomocą testowanych równań, a kątem θ z testu sztywnego koła uzyskiwanym za pomocą wzoru (3.25). Dla różnicy kątów wyliczono również odchylenie standardowe. Wszystkie te parametry statystyczne zebrane są w tabeli 3.4. Implementacja równania f (tabela 3.3) daje najlepszą precyzję przybliżenia promienia (najmniejsze względne odchylenie standardowe promienia o wartości 0.24 w porównaniu z wartością bazową równą 0.13 uzyskaną w teście sztywnego koła). Ta sama implementacja równania daje również najlepsze wyniki pod względem różnicy kątów i odchylenia standardowego tej różnicy. To równanie prowa-dzi do uzyskiwania najbardziej izotropowego kształtu, co jest przedstawione na rysunku 3.11 f, spośród wszystkich testowanych.

a) b) c) d)

e) f) g) h)

Rys. 3.11. Kształty ziaren uzyskane w teście miękkiego koła za pomocą równań (3.32) i (3.33)

zgodnie z implementacją przedstawioną w tabeli 3.3

Analizując dane zawarte w tabeli 3.4 oraz rysunek 3.11 można stwierdzić, że naj-lepsze rezultaty uzyskuje się z wykorzystaniem równania (3.33) dla promienia otoczenia równego a 5. Dalsze jego zwiększanie zarówno pogarsza wyniki, jak i zwiększa ilość