Index of /rozprawy2/10807

Akademia Górniczo-Hutnicza

im. Stanisława STASZICA w Krakowie

WYDZIAŁ ODLEWNICTWA

Katedra Inżynierii Stopów i Kompozytów Odlewanych

Rozprawa doktorska

Opracowanie modelu matematycznego

rzeczywistego pola mikrodyfuzji

dla krystalizacji ziaren równoosiowych

mgr inż. Jacek Początek

Promotor:

Praca została wykonana w ramach projektu badawczego nr DEC-2011/01/B/ST8/01689 finansowanego przez Narodowe Centrum Nauki

(z wyjątkiem rozdziałów 4.5 pt. „Symulacja krystalizacji stopu Pb-Bi w warunkach izotermicznych” i 5.2 pt. „Modelowanie krystalizacji eutektycznej w żeliwie z grafitem kulkowym”, które zostały wykonane w ramach pracy statutowej AGH nr 15.11.170.449).

WYKAZ WAŻNIEJSZYCH OZNACZEŃ ... 4

STRESZCZENIE ... 5

WSTĘP... 6

1. MODELOWANIE MATEMATYCZNE KRYSTALIZACJI RÓWNOOSIOWEJ LIMITOWANEJ PROCESAMI DYFUZYJNYMI. STAN ZAGADNIENIA ... 7

1.1. MECHANIZM RÓWNOOSIOWEJ KRYSTALIZACJI PERYTEKTYCZNEJ IEUTEKTYCZNEJ... 7

1.2. MODELOWANIE SEGREGACJA SKŁADNIKA... 10

1.3. MIKROMODELOWANIE KRYSTALIZACJI... 13

1.4. DIAGRAMY VORONOIA... 21

1.5. METODY ANALIZY POZYCJI FRONTU KRYSTALIZACJI... 23

2. CEL I TEZA PRACY... 29

3. OPIS MODELU... 31

3.1. APROKSYMACJA RÓŻNICOWA GEOMETRII EPMD– UŚREDNIONY WIELOBOK VORONOIA... 31

3.2. NUMERYCZNA METODA ROZWIĄZYWANIA... 34

3.2.1. Parametry EPMD... 35

3.2.2. Transport masy... 36

3.3. RUCH GRANICY MIĘDZYFAZOWEJ... 40

3.4. FUNKCJA ŹRÓDŁA (WYDZIELANIE SIĘ UTAJONEGO CIEPŁA KRYSTALIZACJI) ... 41

4. WYNIKI SYMULACJI KRYSTALIZACJI STOPÓW PERYTEKTYCZNYCH I WERYFIKACJA DOŚWIADCZALNA ... 44

4.1. OPIS UKŁADU TERMODYNAMICZNEGO PB-BI... 44

4.2. KRYSTALIZACJA STOPU DOŚWIADCZALNEGO W WARUNKACH CIĄGŁEGO CHŁODZENIA... 46

4.3. PARAMETRY MODELOWANIA... 48

4.4. WYNIKI SYMULACJI... 48

4.5. WYNIKI SYMULACJI KRYSTALIZACJI IZOTERMICZNEJ STOPU PB –BI... 52

4.6. WALIDACJA WYNIKÓW SYMULACJI... 53

5. WYNIKI SYMULACJI KRYSTALIZACJI STOPÓW FE-C ... 55

5.1. OPIS UKŁADU TERMODYNAMICZNEGO FE-C ... 55

5.2. MODELOWANIE KRYSTALIZACJI EUTEKTYKI W ŻELIWIE SFEROIDALNYM... 56

5.2.1. Parametry modelowania ... 56

5.2.2. Wyniki symulacji ... 56

5.3. MODELOWANIE KRYSTALIZACJI PERYTEKTYKI... 56

5.3.1. Parametry modelowania ... 56

5.3.2. Wyniki symulacji ... 56

5.4. DYFUZJA I CZĘŚCIOWE MIESZANIE SIĘ CIECZY... 56

6. WNIOSKI... 56

Wykaz wa

żniejszych oznaczeń

a, b – współczynniki C – stężenie

cv – objętościowe ciepło właściwe

D – współczynnik dyfuzji f – udział fazy stałej

i, j, k – wektory jednostkowe równoległe do osi współrzędnych

J – strumień dyfuzji

k – współczynnik rozdziału rozpuszczonego pierwiastka pomiędzy fazami stałą

a ciekłą

NA – liczba ziaren na jednostkę powierzchni

n – liczba ziaren przypadająca na jednostkę objętości

q – moc cieplna źródła (człon źródłowy)

r – promień odpowiedniej granicy

r0 – równoważny promień ziarna eutektycznego

T – temperatura

T’ – szybkość stygnięcia

TO – temperatura otoczenia

TL – temperatura likwidus

TP – temperatura przemiany perytektycznej

TE – temperatura przemiany eutektycznej

V – objętość

x1, x2, x3 – współrzędne kartezjańskie

∆H – entalpia przemiany fazowej

f

∆ – przyrost udziału fazy stałej

τ

∆ – krok czasu

λ – współczynnik przewodzenia ciepła

ρ – gęstość

τ – czas

Streszczenie

Pierwsza część pracy zawiera krótki przegląd metod modelowania krystalizacji, począwszy od prostych metod pozwalających na modelowanie rozkładu stężenia za pomocą metod analitycznych, a zakończywszy na rozbudowanych modelach wykorzystujących zaawansowane rozwiązania numeryczne. Przedstawione zostały możliwości zastosowania teselacji Dirichleta i diagramów Voronoia dla opisu kształtu ziaren równoosiowych w stopach metalicznych i geometrii elementarnego pola mikrodyfuzji (EPMD) podczas krystalizacji. Opisano również metody pozwalające na śledzenie przemieszczającego się frontu krystalizacji.

Kolejny rozdział przedstawia w sposób szczegółowy zaproponowany model geometrii rzeczywistego pola mikrodyfuzji podczas krystalizacji ziaren równoosiowych. Model ten oparty został na statystycznej teorii krystalizacji Kolmogorova. Dzięki temu pozwala on na statystyczny opis rzeczywistej geometrii ziaren równoosiowych. Zaproponowany model pozwala, w symulacji procesów krystalizacji limitowanych przez dyfuzję w stanie stałym, zrezygnować ze stosowanej powszechnie idealizowanej geometrii EPMD (kulistej lub walcowej). Model uwzględnia wpływ losowych kontaktów sąsiednich ziaren na kształt EPMD. Cechą charakterystyczną proponowanego rozwiązania jest również to, że nie wymaga przeprowadzania wstępnej teselacji EPMD przed rozpoczęciem obliczeń, co znacznie zmniejsza czas potrzebny na uzyskanie rozwiązania.

W modelu została uwzględniona m.in. dyfuzja pierwiastka stopowego zarówno w zanikającej fazie, jak również w powstających produktach przemiany oraz jej wpływ na warunki panujące na froncie krystalizacji. Uwzględniono efekt dyfuzji wstecznej, co pozwala na dokładniejszy opis niejednorodności składu chemicznego powstających produktów przemian fazowych.

W części eksperymentalnej przedstawione zostały wyniki modelowania krystalizacji nadperytektycznego stopu Pb-Bi (wraz z jego weryfikacją) oraz modelowania krystalizacji stopów Fe-C krystalizujących z przemianą perytektyczną (staliwo) i eutektyczną (żeliwo z grafitem kulkowym).

Wst

ęp

Modelowanie matematyczne, jako proces użycia języka matematyki do opisu zachowania układu (np. termodynamicznego) [1], ma coraz większe znaczenie w przygotowywaniu tech-nologii wykonywania odlewów. Przyczynę coraz częstszego wykorzystania metod matema-tycznych w modelowaniu procesów odlewniczych upatrywać można we wzroście wiedzy teoretycznej dotyczącej podstaw krystalizacji metali i ich stopów. Istotnie znaczenie ma rów-nież ciągłe zwiększanie się mocy obliczeniowej współczesnych komputerów osobistych. Dzięki temu możemy coraz dokładniej i szybciej prognozować cechy makro- i mikrostruktury krzepnącego metalu, a co za tym idzie – właściwości odlewu wykonanego za pomocą dopiero projektowanej technologii.

Poznawane na przestrzeni wieków prawa i zasady z takich nauk jak matematyka, fizyka, lub chemia pozwalają na w miarę dokładny opis podstawowych procesów fizyko-chemicz-nych towarzyszących produkcji odlewów. Niestety, nawet bardzo skomplikowany aparat matematyczny pozwala na uzyskanie rozwiązania analitycznego jedynie dla bardzo upro-szczonych przypadków. W rozwiązaniu zadań bardziej skomplikowanych geometrycznie (bliższych praktyce odlewniczej) niezbędne jest zastosowanie metod numerycznych.

Przedstawiony w pracy mikro-model krystalizacji kontrolowanej przez dyfuzję w przy-padku wzrostu ziaren równoosiowych zajmuje pośrednią pozycję pomiędzy znanymi i rozpowszechnionymi modelami krystalizacji ziaren idealizowanych o kształcie kulistym z jednej strony, a znacznie bardziej szczegółowymi i wymagającymi pod względem mocy obliczeniowej sprzętu komputerowego modelami „ziarnistymi”.

Proponowany model (w odróżnieniu od modeli kulistych) uwzględnia (podobnie jak w przypadku modeli „ziarnistych”) odchylenie kształtu rzeczywistego ziaren równoosiowych od idealizowanego kształtu kulistego. Zastosowanie natomiast mniej szczegółowego, niż w modelach ziarnistych, podejścia statystycznego do opisu rzeczywistego kształtu EPMD stawia zdecydowanie mniejsze wymagania co do sprzętu komputerowego.

1.

Modelowanie matematyczne krystalizacji

równoosiowej limitowanej procesami

dyfuzyjnymi. Stan zagadnienia

1.1. Mechanizm równoosiowej krystalizacji perytektycznej

i eutektycznej

Wzrost ziaren równoosiowych w stopach metali dwu- lub wieloskładnikowych jest możliwy podczas krystalizacji objętościowej zarówno w przypadku wydzielania się jednej fazy stałej z roztworu ciekłego, jak i w przypadku krystalizacji perytektycznej bądź eutektycznej.

Sekwencję krystalizacji w układzie perytektycznym zostanie przedstawiona wykorzystując schemat wykresu równowagi perytektycznej przedstawionej na rys. 1.1 na przykładzie stopu o stężeniu C (Cα < C < CL). W zakresie od temperatury likwidus TL do

temperatury przemiany perytektycznej Tp w roztworze ciekłym zarodkują i wzrastają tylko

ziarna fazy pierwotnej α. Po dalszym obniżeniu temperatury tego stopu poniżej temperatury perytektycznej Tp rozpoczyna się wzrost fazy perytektycznej β, który odbywa się kosztem

zużycia pierwotnej fazy α i cieczy. W trakcie tego procesu można wyróżnić 2 zasadnicze etapy [2-4]:

• Reakcja perytektyczna. Podczas reakcji wszystkie trzy fazy są w kontakcie ze sobą, gdyż tylko wówczas ciecz i przedperytektyczna faza α mogą reagować bezpośrednio ze sobą, dając perytektyczną fazę β. W rezultacie reakcji perytektycznej czoło fazy β przemieszcza się wzdłuż powierzchni ziaren fazy α i szybko otacza ją całkowicie, izolując tym samym od cieczy. Od tego momentu rozpoczyna się 2. etap krystalizacji, czyli tzw. przemiana perytektyczna.

od siebie. Faza β rozrasta się zarówno w kierunku cieczy, jak i fazy pierwotnej, a szybkość jej wzrostu zależy od szybkości dyfuzji przez powstającą warstwę.

T e m p e ra tu ra Stężenie składnika B

Rys. 1.1. Schemat fragmentu wykresu równowagi fazowej krystalizującego z przemianą perytektyczną

W miarę upływu czasu ciecz zanika, natomiast faza β powiększa swoje wymiary. Zależnie od składu chemicznego stopu otrzymać można stop jednofazowy β lub dwufazowy α+β. Mechanizm krystalizacji, który w danej chwili przeważa, zależy od parametrów termofizycznych faz, składu chemicznego, jak również od różnicy w objętości właściwej poszczególnych faz [5]. Z reguły reakcja perytektyczna występuję na początku, a po szybkim odizolowaniu ziaren fazy pierwotnej od cieczy krystalizacja rozwija się według mechanizmu przemiany perytektycznej.

Podobny mechanizm krystalizacji równoosiowej można zaobserwować również w stopach eutektycznych. Na rodzaj uzyskanej eutektyki wpływają:

• entropia rozpuszczania faz eutektycznych,

• prędkość wzrostu i gradient temperatury,

• modyfikacja stopu.

Dokładniejszy opis poszczególnych czynników na rodzaj uzyskanej eutektyki można znaleźć w [4]. T e m p e ra tu ra Stężenie składnika B

Rys. 1.2. Schemat fragmentu wykresu równowagi fazowej krystalizującego z przemianą eutektyczną

Występowanie eutektyki globularnej stwierdzono tylko w nielicznej grupie stopów, z których największe znaczenie mają stopy oparte na bazie żelaza z węglem Fe-C. Stopem tym jest żeliwo z grafitem kulkowym (żeliwo sferoidalne). Podczas krystalizacji eutektycznej z cieczy wydzielają się ziarna grafitu w postaci kulistej. Następnie na ich powierzchni szybko pojawia się ciągła warstwa austenitu, oddzielająca grafit od cieczy. Pomimo odcięcia grafitu od cieczy, jego ziarna nie przestają rosnąć. Kontynuuje on swój wzrost wskutek dyfuzji węgla z cieczy poprzez otoczkę austenitu. Budowę ziarna eutektyki globularnej charakteryzują proste parametry geometryczne, takie jak: promień kulki grafitowej i promień otoczki austenitu, będący równocześnie promieniem ziarna eutektycznego [4, 6].

Analiza początkowych etapów wzrostu eutektyki globularnej i przemiany perytektycznej wykazuje, że są to zjawiska podobne pod względem geometrii EPMD. W obu przypadkach najpierw wzrasta jedna faza, a następnie zostaje ona otoczona przez drugą fazę. Różnica

krystalizacji opisanych powyżej stopów można zastosować wspólny model matematyczny krystalizacji równoosiowej, przewidujący modelowanie dyfuzji składnika przez rosnącą warstwę fazy stałej pomiędzy zanikającą fazą ciekłą a ziarnem fazy stałej, zlokalizowanym w środku ziarna równoosiowego.

1.2. Modelowanie segregacja składnika

Opis segregacji składnika (domieszki) w stopach podczas przemian fazowych wymaga określenia zmian w czasie pola stężenia składnika w fazie stałej oraz w cieczy, czyli wyznaczenia funkcji, którą ogólnie dla układu jednowymiarowego można zapisać w postaci [4]:

( )

x t f

C= , (1.1)

gdzie:

C – stężenie składnika w punkcie x i w czasie t.

Jedną z pierwszych, a zarazem jedną z najprostszych, metod pozwalających na analizę rozkładu stężenia podczas procesu krystalizacji była reguła dźwigni. Zakładała ona idealną (nieskończenie szybką) dyfuzję zarówno w fazie stałej, jak i w ciekłej [4]. W przypadku krystalizacji jednofazowej stężenie składnika w fazie stałej CS można wyznaczyć za pomocą

następującego równania wynikającego z bilansu masy tego składnika:

(

1

)

1 0 − + = k fs kC C_S (1.2) gdzie:

k – równowagowy współczynnik rozdziału (domieszki) pierwiastków w krystali-zującej cieczy (k = CS/CL),

fs – udział fazy stałej.

Niestety, z powodu małej szybkości dyfuzji większości składników w stopach metalowych, przybliżenie to nie pozwala na zastosowanie tego schematu w rzeczywistych warunkach. Metoda ta została rozwinięta przez Gullivera w 1913 r. [7], a następnie ponownie zastosowana przez Scheila w 1942 r. [8]. Model ten jest dziś nazywany jako model

Gullivera-Scheila (GS). Zakłada on nieskończenie szybką dyfuzję w cieczy oraz brak dyfuzji w stanie stałym, jak również brak ograniczeń co do warunków geometrycznych. Przyjmuje się również stężenia równowagowe na granicach faz. Zgodnie z modelem GS stężenie składnika w fazie stałej określa następująca zależność:

(

)

1 0 1 − − = k S S kC f C (1.3)

Największą wadą modelu było pominięcie dyfuzji w fazie stałej. Z tego powodu model ten nie mógł być zastosowany w stopach, w których występuje stosunkowo szybka dyfuzja w fazach stałych, ani w przypadkach, kiedy proces krystalizacji trwa bardzo długo.

Reguła dźwigni, tak samo jak model GS mogą być stosowane wyłącznie dla symulacji krystalizacji jednofazowej. Pomiędzy tymi dwoma przypadkami mogą występować pośrednie, które uwzględniają zarówno limitowaną szybkość dyfuzji w stanie stałym, jak i w stanie ciekłym. Przykładem takiej pracy może posłużyć publikacja Brodiego i Flemingsa (BF) [9], gdzie zaproponowano ogólniejszą formę dla modelu GS, w której uwzględniono całkowitą dyfuzję w cieczy oraz częściową dyfuzję wsteczną w fazie stałej za pomocą równania:

(

)

[

]

( ) (k k) S S kC k f C = ₀1− 1−β −1/1−β (1.4)

Parametr β został nazwany parametrem wstecznej dyfuzji i przez wielu kolejnych badaczy był wyznaczany na różne sposoby. Pierwotnie Brody i Flemings założyli, że:

α =

β 2 (1.5)

gdzie: α – bezwymiarowy współczynnik dyfuzji (bezwymiarowa liczba Fouriera 2. rodzaju dla dyfuzji masy), wyrażony jako:

2 L t DS = α (1.6) gdzie:

Ds – współczynnik dyfuzji segregującego pierwiastka w fazie stałej,

Równanie (1.4) upraszcza się do modelu Scheila przedstawionego równaniem (1.3), kiedy współczynnik dyfuzji w fazie stałej dąży do zera (brak dyfuzji w stanie stałym). Z drugiej strony, jeżeli założyć, że współczynnik dyfuzji w fazie stałej dąży do nieskończoności, wydawać się powinno, że model powinien przyjąć postać reguły dźwigni. Nie dzieje się tak, dlatego, że przy dużych DS (β > 1) model nie zachowuje bilansu masy. Aby zapewnić bilans

masy, Ohnaka [10] zaproponował, aby współczynnik β był wyznaczany jako:

α + α = β 2 1 2 (1.7)

Współczynnik ten został wyznaczony poprzez porównanie przybliżonego rozwiązania równania dyfuzji dla płytki dendrytu, zakładając kwadratowy profil stężenia fazy stałej. Ohnaka zaproponował również inną postać parametru β, uwzględniającą m.in. nieregularną mikrostrukturę, w postaci: α + α = β 4 1 4 (1.8)

Ohnaka wykazał, że zaproponowana przez niego nowa postać parametru β pozwala uzyskać lepsze przybliżenie wyników otrzymanych z doświadczeń eksperymentalnych. Clyne i Kurz w pracy [11] zaproponowali kolejną modyfikację modeli GS i BF w postaci:

      α − −             α − − α = β 2 1 exp 2 1 1 exp 1 (1.9)

która pozwala na asymptotyczne zbliżanie się zarówno do równania Scheila, jak i do reguły dźwigni dla odpowiednio nieskończenie małych oraz nieskończenie dużych współczynników dyfuzji. Model ten został jednak skrytykowany w [12] za brak dostatecznych podstaw fizycznych do modyfikacji parametru β. W pracy [13] został zaprezentowany numeryczny model krystalizacji potrójnego stopu Al-Li-Cu z użyciem modelu krystalizacji opracowanego przez Brodiego i Flemingsa [9]. Pomimo przyjęcia pewnych uproszczeń, takich jak płaski front krystalizacji i stała szybkość wzrostu fazy stałej, model pozwalał na dość dokładne wyznaczenie rozkładu stężenia w poszczególnych fazach. Model Brody-Flemingsa [9] był rozwijany również przez innych naukowców, takich jak Kobayashi [14], Nastac i Stefanescu [15], Solari i Biloni [16], Wołczyński i in. [17, 18].

Niestety, wszystkie wymienione wyżej modele są użyteczne tylko w początkowym etapie krystalizacji, kiedy ziarna fazy stałej jeszcze nie stykają się ze sobą. Nie uwzględniają one również zmian pola temperatury. Z tego powodu, niezbędne stało się stworzenie takiego modelu, który w sposób analityczny lub numeryczny rozwiązywałyby zagadnienia związane z przepływem ciepła, wzrostem faz (śledzenie frontu krystalizacji) oraz pozwalałby wyznaczyć rozkład stężenia danego pierwiastka w poszczególnych fazach, biorąc pod uwagę kontakty pomiędzy sąsiednimi ziarnami.

1.3. Mikromodelowanie krystalizacji

Jedną z pierwszych, przełomowych prac w zakresie modelowania mikrostruktury stopów powstającej podczas krystalizacji, jest praca [19], w której został opisany analityczny model krystalizacji eutektyki płytkowej i włóknistej za pomocą równania dyfuzji dla ruchomej granicy międzyfazowej (która znajduje się w warunkach izotermicznych), z uwzględnieniem płaskiego frontu krystalizacji. W celu walidacji wyników otrzymanych za pomocą modelowania posłużono się stopem Pb-Sn oraz organicznymi materiałami przeźroczystymi. Porównanie to wykazało dobrą zgodność modelu z eksperymentem. W kolejnych latach powstało szereg prac, których celem było wyeliminowanie niektórych ograniczeń tego modelu. Zgodnie z informacjami zawartymi w pracy [20], Tiller [21] rozwinął model Jacksona-Hunta tak, aby nie był on już ograniczony tylko do krystalizacji kierunkowej, natomiast Wetterfall i in. w [22] zaprezentowali model analityczny wzrostu kulki grafitu przez otoczkę austenitu w żeliwie z grafitem kulkowym. Szybkość wzrostu fazy γ została opisana jako:

(

)

γ γ γ γ γ γ γ γ − − − = L L L G G G C C C C r r r r D dt dr (1.10) gdzie:

D – współczynnik dyfuzji węgla w austenicie,

C – stężenie węgla na odpowiedniej granicy (indeksy dolne: G/γ – w graficie na granicy z austenitem, L/γ – w cieczy na granicy z austenitem, γ/L –

Model ten używany jest do tej pory w modelowaniu zmian mikrostruktury. Niestety, wymaga on założenia, że stosunek promienia austenitu do promienia grafitu pozostaje stały podczas krystalizacji. Magnin i Kurz w pracy [23] założyli nieizotermiczną powierzchnię interfejsu, czyli granicy pomiędzy fazą stałą i ciekłą. Również Cataliy i in. w [24]zrezygnowali z powierzchni izotermicznej na interfejsie oraz przyjęli dodatkowo różne gęstości dla fazy ciekłej i obu faz stałych.

Jedną z pionierskich prac związanych z numerycznym modelowaniem krystalizacji była praca Murraya i Landisa [25] opublikowana w 1959 r. Dotyczyła ona obliczeń przepływu ciepła podczas topienia oraz krystalizacji czystych metali z wykorzystaniem metody różnic skończonych w przestrzeni jednowymiarowej. Założono w niej również stałe, ale różne dla każdej z faz, właściwości termiczne. W celu uniknięcia numerycznych problemów związanych z migracją granicy międzyfazowej, autorzy zastosowali siatkę różnicową ze zmiennym krokiem przestrzennym. Metoda ta była stosowana w późniejszych pracach, m.in. w [26, 27]. Jej bardziej szczegółowy opis przedstawiono w rozdziale 1.5 niniejszej rozprawy.

W 1966 r. powstała druga, również bardzo ważna, praca z zakresu modelowania autorstwa Oldfielda [28]. Praca ta pozwalała na wyznaczenie krzywej chłodzenia żeliwa z grafitem płatkowym. Uważa się również, że jest to jeden z pierwszych modeli numerycznych, który umożliwił otrzymanie danych dotyczących mikrostruktury krzepnącego stopu. Praca Oldfielda wprowadzała szereg innowacyjnych rozwiązań. Zakładała ona m. in., że prędkość wzrostu ziaren eutektycznych jest wprost proporcjonalna do kwadratu przechłodzenia oraz, że liczba ziaren eutektycznych jest zależna tylko od kwadratu przechłodzenia poniżej temperatury równowagi i nie zależy od czasu (rys. 1.3). Na krzywych chłodzenia uzyskanych za pomocą tego modelu obserwuje się lokalny wzrost temperatury w początkowym etapie krystalizacji stopu (tzw. rekalescencja, patrz rys. 1.3). Jest to zjawisko typowe dla krystalizacji objętościowej wielu stopów odlewniczych.

Kolejni naukowcy rozwijali w kolejnych latach model Oldfielda. W 1974 r. Stefanescu i Trufinescu w [29] zbadali wpływ modyfikacji stopu na krzywą chłodzenia i stałą zarodkowania.

W pracy [30] Tanzilli i Heckel zaproponowali numeryczną metodę analizy rozkładu stężenia pierwiastka kontrolowanego przez dyfuzję w jednowymiarowym obszarze dwufazowym. W celu śledzenia płaskiego, sferycznego i cylindrycznego frontu krystalizacji została użyta metoda różnic skończonych. Określono wpływ poszczególnych parametrów wejściowych (współczynników dyfuzji, geometrii, stężeń równowagowych) zarówno na kinetykę krystalizacji, jak i na homogenizację stopu.

Kolejnymi ważnymi pracami są prace Fredrikssona i Svenssona [31,32] opublikowane w 1985 r. i w 1988 r., w których połączony został analityczny model przepływu ciepła z parabolicznym prawem wzrostu eutektyki grafitowej i węglikowej w żeliwie oraz została wzięta pod uwagę dyfuzja węgla z cieczy, przez otoczkę austenitu, do kulki grafitu w żeliwie z grafitem kulkowym. Uwzględniono także zderzania się sąsiednich ziaren ze sobą. W dalszych latach model Oldfielda był rozwijany przez kolejnych autorów, takich jak Stefanescu [33], Fraś [34]. T e m p e ra tu ra , o C Czas, s

Rys. 1.3. Eksperymentalne i obliczeniowe krzywe chłodzenia oraz równania opisujące gęstość ziaren i ich prędkość wzrostu wg. Oldfielda [28]; linia ciągła wyniki eksperymentu, linia

przerywana wyniki symulacji

Kolejnym etapem w rozwoju modeli krystalizacji były modele, w których zarówno przepływ ciepła jak i masy był wyznaczany w sposób numeryczny. Jednym z pierwszych takich modeli był model z dyfuzją ustaloną, opracowany w 1985 r. przez Su i in. [35]. Dotyczył on krystalizacji żeliwa z grafitem kulkowym i pozwalał na rozwiązanie wszystkich równań za pomocą metody różnic skończonych. W tym modelu zastosowano prawo zarodkowania Oldfielda oraz uwzględniono dyfuzję węgla przez warstwę austenitu. Powyższy model został rozwinięty przez Rappaza i in. w pracy [36], gdzie została zaproponowana następująca modyfikacja równania (1.10):

(

)

(

)

γ γ γ γ γ γ γ γ γ γ γ γ − − ρ ρ + − − = L L L L G L G G C C r r r dt dC r r r C C r D dt dr 2 3 3 0 3 (1.11) gdzie:

r₀ – równoważny promień ziarna eutektycznego.

Model ten uwzględnia niestacjonarną dyfuzję w fazie stałej, lecz pomija dyfuzję w cieczy. Chen i in. [37] zmodyfikowali równanie (1.11) aby uniknąć błędów, jakie wynikały z przyjętego założenia całkowitego mieszania się węgla w cieczy.

Kolejnymi pracami, które włączały równanie Fouriera do modelowania mikrostruktury, były prace [38, 39]. W pracy [40] zaprezentowany został szczegółowy model krystalizacji dendrytów równoosiowych, uwzględniający zarodkowanie, kinetykę wzrostu i przepływ ciepła. Model wprowadzał sferyczną otoczkę wokół dendrytu, która oddziela ciecz międzydendrytyczną (ang. interdendritic liquid) od cieczy znajdującej się na zewnątrz tej otoczki (ang. extradendritic liquid). Zastosowane zostały dwa różne kroki czasowe, jeden dla opisu zmian pola temperatury w skali makro, drugi, znacznie mniejszy, do opisu procesu zarodkowania i wzrostu fazy stałej w skali mikro. Niestety, model dotyczył tylko przestrzeni jednowymiarowej z uwzględnieniem symetrii osiowej.

W pracy [41] Stefanescu i in. zaproponowali model krystalizacji łączący makroskopowy przepływ ciepła z mikroskopową kinetyką krystalizacji. Wprowadzili również dokładniejszy opis ciepła krystalizacji oraz zaproponowali natychmiastowe zarodkowanie ziaren, gdy tylko temperatura cieczy obniża się poniżej wcześniej założonej temperatury zarodkowania.

Jednym z pierwszych modeli, pozwalającym opisać krystalizację równoosiową w skali mikro-makro, był już wcześniej wspomniany model opracowany przez Rappaza i Thevoza [39]. W kolejnych pracach [42, 43] zaproponowano wykorzystanie prawa normalnego rozkładu statystycznego (Gaussa) do modelowania procesu zarodkowania ziaren fazy stałej. Sposób ten polegał na podzieleniu obszaru na regularną siatkę sześcianów, podobnie jak ma to miejsce w automatach komórkowych. Wzrost każdego ziarna był opisywany przez zestaw parametrów takich jak: pozycja środka ziarna, wektor wiodący oraz status komórki (czy komórka graniczy z cieczą, czy z innym ziarnem). Technika ta pozwalała na otrzymanie dobrych wyników, lecz charakteryzowała się bardzo długim czasem obliczeń.

W 1995 r. Trivedi opublikował pracę [44], w której zaproponował teoretyczny model wzrostu struktur pasmowych, powstających w wyniku przemiany perytektycznej. Według

autora tej pracy szerokość warstw, powstałych podczas krystalizacji, jest odwrotnie proporcjonalna do prędkości przesuwania się frontu krystalizacji. Współczynnik ten związany jest z przechłodzeniem, wymaganym do zarodkowania zarówno fazy pierwotnej, jak i perytektycznej oraz z pierwotnym składem stopu. Model ten został rozszerzony w pracy [45] o zjawiska konwekcji. Uniknięto założenia stałej prędkości przemieszczania się interfejsu.

Analityczny model kinetyki przemiany perytektycznej, oparty na liniowym rozkładzie stężenia składnika rozpuszczonego w fazie perytektycznej, został zaprezentowany w pracy [46]. Das i współautorzy zaprezentowali prosty model matematyczny, który cechował się dobrą zgodnością z wynikami eksperymentu. Niestety, pod koniec procesu krystalizacji, różnica pomiędzy szybkością reakcji, otrzymaną z modelowania, a danymi eksperymentalnymi powiększała się. Według autorów, odchyłka modelu od rzeczywistej geometrii wynikała z odchylenia kształtu elementarnego pola dyfuzji od geometrii idealnej kuli. Ci sami autorzy, w pracy [47], zaproponowali bardziej rygorystyczny model numeryczny, opisujący kinetykę przemiany perytektycznej, bazujący na dyfuzji nieustalonej. Założyli oni również nieliniowy rozkład stężenia w fazie stałej, co pozwoliło na uzyskanie lepszego dopasowania do wyników eksperymentalnych.

W 1996 r. opublikowano artykuł [48], w którym został zaprezentowany model makro-segregacji podczas przemiany perytektycznej w stali niskowęglowej. Celem pracy był opis zależności, panujących pomiędzy typem reakcji perytektycznej a wielkością segregacji i rodzajem naprężeń, występujących podczas krystalizacji. Zostało udowodnione m. in., że skurcz termiczny generowany przez reakcję perytektyczną i naprężenia rozciągające generowane przez przemianę w stanie stałym są proporcjonalne do szybkości chłodzenia, oraz że makrosegregacja, szczególnie przy powierzchni wlewka, jest bardzo wrażliwa na szybkość odprowadzania ciepła.

Model opracowany przez zespól Frasia [49], oparty został na modelowaniu mikro-makro i dotyczył krystalizacji żeliwa z grafitem kulkowym. Zakładał on, że na sferycznej kulce grafitu zaczyna rosnąć otoczka austenitu. W pierwszym etapie krystalizacji wzrost ziaren jest kontrolowany przez przechłodzenie na powierzchni faza stała-ciecz i stacjonarną dyfuzję w cieczy. Podczas trwania tego etapu krystalizacji przyjęty został stały stosunek promienia austenitu do grafitu wynoszący 2,4 [22]. Zgodnie z [50] proces ten kończy się w momencie osiągnięcia maksymalnego przechłodzenia. Dalszy wzrost fazy stałej kontrolowany jest przez

ziaren i dyfuzji składników. Model pozwala na wyznaczenie zarówno krzywej stygnięcia, promieni poszczególnych faz, jak również rozkładu stężenia węgla w austenicie i cieczy.

Matsumiya i współautorzy [51] opracowali model krystalizacji, w którym badany obszar krystalizacji został podzielony na komórki w kształcie sześciokątów lub kwadratów. Analizę swoją oparli na uwzględnieniu pełnego mieszania cieczy i wykorzystaniu prawa Ficka w fazie stałej. Przyjęte zostało, że interfejs będzie poruszał się o jeden węzeł siatki, co wymagało odpowiedniego dostosowania kroku czasowego. Oznacza to, że dla danego kroku czasowego problem był iterowany w ten sposób, aby po upływie określonego czasu, temperatura w sąsiedniej komórce (składającej się z cieczy) odpowiadała temperaturze z krzywej chłodzenia. Dla rozwiązania układu równań wykorzystano metodę różnic skończonych, i w rezultacie uzyskano wyniki zbliżone do wyników z eksperymentu rzeczywistego. Model ten był rozwijany dalej, m. in. w pracy [52], gdzie posłużył do opisu zjawisk mikrosegregacji podczas przemiany perytektycznej z wykorzystaniem stopu wieloskładnikowego (Fe, C, Si, Mn, P, S) oraz w pracy [53], gdzie został połączony z komercyjnym programem ChemAPP w celu wyznaczenia równowagowego stężenia na interfejsie. W tym przypadku obszar obliczeniowy był przybliżany za pomocą kształtu gwiazdy (rys. 1.4), a do symulacji posłużył stop, w skład którego wchodziły takie pierwiastki jak Fe, C, Mn, Si, P, Mo. Natomiast w pracy [54] został przedstawiony jednowymiarowy model łączący makrosegregację (uwzględniającą segregację w skali całego odlewu) z mikrosegregacją (segregacja w skali jednego ziarna) w wieloskładnikowym stopie z przemianą perytektyczną. W modelu tym wykorzystano również wartość stężenia równowagowego danego pierwiastka na interfejsie, korzystając z termodynamicznych układów równowagi fazowej, obliczanych dzięki komercyjnemu programowi Thermo-Calc. Niestety, wadą tego modelu jest bardzo duża czasochłonność obliczeń, a dla jej uniknięcia niezbędne jest użycie jak najmniejszej liczby węzłów w siatce odpowiedzialnej za mikrosegregację.

W pracach [55, 56] dotyczących krystalizacji dendrytów kolumnowych jak i równo-osiowych Wang i Beckermann zaproponowali obszar tzw. objętości kontrolnej (ang. control

volume), która składała się zarówno z fazy stałej, jak i dwóch obszarów fazy ciekłej (cieczy

międzydendrytycznej i cieczy znajdującej się na zewnątrz dendrytu). Model zakładał, że obszar ten będzie jednowymiarowym sferycznym polem o promieniu R. Podejście to zostało ponownie wykorzystane w pracy [57] do modelowania krystalizacji wielofazowej. W pracy tej założono, że dla każdej nowo powstałej fazy tworzona jest nowa objętość, dla której są prowadzone obliczenia.

Rys. 1.4. Schemat ilustrujący podłużne i poprzeczne przekroje dendrytu wykorzystane w pracy [53]

Kiedy pojawia się drugi zarodek, automatycznie tworzony jest kolejny obszar, w którym centralną część zajmuje powstała faza stała. Model zakładał wymianę masy zarówno pomiędzy fazami (w jednej objętości), jak również pomiędzy sąsiednimi objętościami kontrolnymi. W wyniku uwzględnienia dyfuzji dla każdej z faz oraz kinetyki wzrostu każdego składnika mikrostruktury można było uzyskać zmianę udziałów objętościowych poszczególnych frakcji w funkcji czasu. Metoda uśrednionej objętości została ponownie użyta w pracy [58]. W tym przypadku dotyczyła ona modelowania krystalizacji stopu wieloskładnikowego z przemianą perytektyczną. W pracy zostało przyjęte, że zarodkowanie nowej fazy następuje w centralnej części uśrednionej objętości, po czym następował jego

sferyczny wzrost. Ze względu na dużą wartość liczby Lewisa1 temperatura była ujednorodniana w całym układzie.

Rys. 1.5. Schematyczne przedstawienie objętości kontrolnej składającej się fazy stałej, cieczy międzydendrytycznej i cieczy, znajdującej się na zewnątrz dendrytu [55]

Jak widać z przeprowadzonej powyżej analizy, w celu modelowania wzrostu ziaren równoosiowych limitowanego dyfuzją, krystalizujących z przemianą perytektyczną lub jako eutektyka globularna, mogą zostać wykorzystane modele, wykorzystujące tzw. sferyczne elementarne pole mikrodyfuzji. Podejście takie jest znane także w innych pracach, m. in. w: [59-63].

Niestety, rzeczywisty kształt krystalizujących ziaren różni się zdecydowanie od kształtu kulistego, nawet jeżeli weźmie się pod uwagę natychmiastowe zarodkowanie i symetryczny wzrost ziaren. Spowodowane jest to losowym rozkładem zarodków w przestrzeni. W tym celu, aby uwzględnić wpływ kontaktu pomiędzy sąsiednimi ziarnami, można skorzystać z podziału przestrzeni za pomocą teselacji Dirichleta na tzw. wielościany Voronoia.

1

Liczba Lewisa – liczba podobieństwa, określająca stosunek dyfuzji termicznej (Dth) do dyfuzji molekularnej (Dmol): Le=Dth Dmol

1.4. Diagramy Voronoia

W matematyce podział przestrzeni na diagramy Voronoia opisuje sposób pokrycia płaszczyzny wielokątami wzajemnie przylegającymi, ale nie zachodzącymi na siebie. W przypadku płaszczyzny, na której znajduje się zbiór n punktów, nazywanych dalej środkami komórek, dzieli się jej powierzchnię na n obszarów wielokątnych w taki sposób, że każdy punkt w dowolnym z tych wielokątów znajduje się bliżej określonego środka ze zbioru

n punktów, niż od pozostałych n-1 punktów. Utworzone w ten sposób wielokąty mają nazwę

wielokątów lub komórek Voronoia, a procedura podziału jest nazywana teselacją Dirichleta lub teselacją Voronoia. Aby wyznaczyć odcinki linii prostych, oddzielające poszczególne komórki Voronoia od siebie, należy prostą symetryczną podzielić odcinek pomiędzy położonymi najbliżej siebie punktami. Po wykonaniu tej operacji dla każdej najbliższej pary punktów powstaną odcinki, które w miejscach przecięć utworzą wierzchołki wielokątów.

Przykładowy rozkład komórek Voronoia dla zbioru losowych punktów został przedstawiony na rys. 1.6. Podobny podział (teselację) można zastosować dla przestrzeni trójwymiarowej. Wynikiem takiej operacji będą wielościany Voronoia.

Rozwiązanie to jest bardzo często stosowane w modelach matematycznych dla celów reprezentacji mikrostruktury równoosiowej metali i stopów [64-71] ze względu na jej topolo-giczne podobieństwo z mikrostrukturą metali. Zgodnie z [72] każda wygenerowana komórka Voronoia, tak samo jak pojedyncze ziarno, spełnia kryterium Eulera:

2

= + −E F

V (1.12)

jak również zależność:

4 2 − = F V (1.13) gdzie: V – liczba wierzchołków, E – liczba krawędzi,

Ponadto, zgodnie z [73], odpowiednie wartości kątów dwuściennych2 oraz kątów wiązania3 dla komórki Voronoia spełniają kryterium minimalnej energii powierzchniowej, natomiast w pracy [74] wykazano, że wielościany Voronoia są zawsze wypukłe, ale nie koniecznie muszą mieć ograniczoną powierzchnię.

Rys. 1.6. Komórki Voronoia dla losowego zbioru punktów na

płaszczyźnie [http://pl.wikipedia.org/wiki/Diagram_Woronoja. Pobrano: 2014-06-12]

Równoosiowa struktura ziaren, powstająca podczas krystalizacji stopu, może być przybliżona za pomocą wielościanów Voronoia. W miejscach, gdzie ziarna rosnące z jedna-kową prędkością (pierwotnie kuliste) stykają się ze sobą, powstają odcinki granicy płaskiej pomiędzy tymi ziarnami. Przykład takiego podejścia dla krystalizacji dwuwymiarowej został przedstawiony w pracach [69, 75]. Proces krystalizacji i mikrosegregacji w tych pracach był analizowany w obszarach trójkątnych, tworzonych przez połączenie wierzchołków komórki Voronoia z miejscem zarodkowania ziarna (jej środkiem). Powyższe podejście uzyskało nazwę „granular model” (model ziarnisty) i może być zastosowane również w przestrzeni trójwymiarowej, gdzie zamiast wielokątów będą użyte wielościany Voronoia [70]. Modele krystalizacji oparte na diagramach Voronoia można również znaleźć w pracach [76, 77].

2 _{Kąt pomiędzy dwiema półpłaszczyznami posiadającymi wspólną krawędź.} 3

W modelach ziarnistych dla poszczególnych obszarów trójkątnych (2D) lub czworościen-nych (3D) stosowana jest indywidualna analiza numeryczna jednowymiarowego pola mikrodyfuzji, co wymaga dużej mocy obliczeniowej i dłuższego czasu symulacji.

Wybrane cechy geometrii wielościanów Voronoia wykorzystano w niniejszej pracy.

1.5. Metody analizy pozycji frontu krystalizacji

Po ochłodzeniu cieczy poniżej temperatury likwidus rozpoczyna się krystalizacja fazy stałej. Ziarna fazy stałej są oddzielone od fazy ciekłej ruchomą granicą. Ruch granicy międzyfazowej, nazywanej również interfejsem „faza stała-ciecz”, powoduje zwiększanie się udziału objętościowego fazy stałej. Klasycznym zagadnieniem matematycznym, opisującym zmianę położenia granicy międzyfazowej, jest tzw. problem Stefana. Problem Stefana zwany również zadaniem Stefana jest opisem matematycznym zjawiska przemieszczania się ruchomej granicy międzyfazowej z uwzględnieniem procesów transportu ciepła w materiałach, w których zachodzą przemiany fazowe tj. topnienie i krystalizacja.

Analityczne rozwiązanie tego problemu jest skomplikowane, wymaga bowiem ustalenia nieznanego położenia granicy międzyfazowej z uwzględnieniem nieliniowych zmian warunków brzegowych na tej granicy. Pierwsze rozwiązanie tego problemu zaproponował Stefan [78] dla zadania przemarzania i roztapiania gruntu w warunkach stałej temperatury przemiany fazowej i braku strefy dwufazowej. Bardziej zaawansowane rozwiązanie analityczne zadania, w którym uwzględniono występowanie typowej dla stopów metali strefy dwufazowej, w której występują równocześnie dwie fazy: stała i ciekła, zostało podane w pracy [79].

Analityczne rozwiązanie problemu Stefana, jak już to zostało wspomniane wyżej, jest znane, ale tylko do prostej geometrii zadania. Niestety, wspomnianych rozwiązań nie da się zastosować w zadaniach praktycznych, występujących w przemyśle hutniczym i odlewnic-twie. Z tego powodu do celów symulacji krystalizacji metali i stopów z uwzględnieniem faktycznego kształtu geometrycznego produktów przemysłowych wykorzystuje się metody numeryczne.

rozwój metod symulacji procesu krystalizacji: uwzględnienie złożonej geometrii zadania oraz wprowadzenie modeli matematycznych opisujących przejście ze stanu ciekłego do stanu stałego na poziomie mikrostruktury.

Przedstawiony poniżej krótki przegląd metod numerycznych obejmuje tylko te, które najczęściej występują w literaturze, i nie obejmuje metod, które np. stanowią połączenie metod klasycznych.

Zgodnie z klasyfikacją zaprezentowaną w [80] rozwiązanie problemu ruchomej granicy może nastąpić za pomocą następujących metod:

• metoda unieruchomienia frontu (ang. front-fixing)

• metoda front-obszar (ang. front-domain)

• metoda śledzenia frontu (ang. front-tracking)

W metodzie unieruchomienia frontu zaprezentowanej m. in. w [81, 82] zastosowana jest transformacja współrzędnej krzywoliniowej, polegająca na takim przekształceniu fizycznych granic analizowanej przestrzeni, aby poruszający się front krystalizacji, znajdujący się pomiędzy tymi granicami, miał cały czas takie same współrzędne. Metoda ta dobrze sprawdza się w prostych konfiguracjach, np. dla powierzchni płaskich, cylindrycznych lub sferycznych. Niestety, dla bardziej skomplikowanych kształtów użycie tej metody jest nieefektywne, ze względu na skomplikowane równania transformacyjne, niezbędne do opisu ruchu węzłów siatki.

Jedną z najbardziej znanych metod należącej do grupy metod front-obszar jest metoda entalpowa [83, 84]. Metoda ta polega na takim przekształceniu równania energii, aby w opisie matematycznym procesów cieplnych zamiast temperatury pojawiła się entalpia, zdefiniowana jako [85-87]:

( )

=

_∫

ρ +ρ

(

−

( )

)

T T S S Od T f L dT c T H 1 (1.14)

gdzie Tod jest temperaturą odniesienia.

Różniczkując równanie (1.14) względem czasu, otrzymujemy:

t f L t T c t H _S S _∂ ∂ ρ − ∂ ∂ ρ = ∂ ∂ (1.15) lub t T c t H eqv _∂ ∂ = ∂ ∂ (1.16)

gdzie zmienna ceqv, nazywana zastępczym ciepłem właściwym, jest zdefiniowana następująco: T f L c c S S eqv _∂ ∂ ρ − ρ = (1.17)

Metody wyznaczania zastępczego ciepła właściwego dla symulacji krzepnięcia za pomocą metody entalpowej opisane zostały m. in. w [88].

Wykorzystanie powyższych zależności pozwala zastąpić rozwiązywanie nieliniowego równania różniczkowego przewodzenia ciepła w zadaniu krystalizacji rozwiązaniem liniowego równania różniczkowego:

(

T

)

t T ceqv =∇⋅ λ∇ ∂ ∂ (1.18)

Powyższe równanie, po jego przekształceniu do zbioru równań różnicowych, jest rozwiązywane numerycznie z wykorzystaniem odpowiednich metod (najczęściej – metody różnic skończonych lub elementów skończonych). Metoda ta nie wymaga ciągłego śledzenia pozycji interfejsu i pozwala na modelowanie nieizotermicznych procesów krzepnięcia, czyli krystalizacji w zakresie temperatury.

Zarówno pozycja, jak i szerokość strefy dwufazowej mogą zostać wyznaczone po zakończeniu obliczeń na podstawie analizy wyników symulacji pola temperatury. Wadą tej metody są również trudne do uniknięcia wahania pola temperatury, pojawiające się w pobliżu interfejsu. Głównym powodem wahań jest niska gęstość siatki obliczeniowej. Niestety, zwiększenie liczby węzłów powoduje wydłużenie czasu oczekiwania na wynik.

Wartą uwagi jest również metoda zbiorów poziomicowych (ang. level-set), która coraz częściej używana jest nie tylko w rozwiązywaniu problemu Stefana, ale również w opisie przepływu cieczy oraz do opisu kinetyki wzrostu kryształów [89, 90]. Jej idea polega na analizie i wyznaczeniu ruchu interfejsu Γ

( )

t obszaru Ω w polu prędkości υ. W 1988 roku Osher i Sethian [91] zdefiniowali gładką funkcję φ(x, t), nazwaną zbiorem funkcji poziomic (ang. level set function). Granica Γ

( )

t określona jest jako zbiór, gdzie funkcja Φ

( )

x,t =0,

(

x x x

)

R

x= ₁, ₂K, _n ∈ _{. Oznacza to, że zerowy poziom tej funkcji odpowiada interfejsowi} Γ

( )

_{t ,} natomiast dla punktów, znajdujących się w jednolicie dyskretyzowanym obszarze, przypisuje

( )

( )

( )

( )

     Γ = Ω ∂ ∈ = φ Ω ∈ < φ Ω ∈ > φ t x t x x t x x t x dla dla dla 0 , 0 , 0 , 2 1 (1.19)

gdzie: Ω₁∪Ω₁∪∂Ω=Ω. Zgodnie z (1.19) odległości w ciele stałym mają znak ujemny, natomiast odległości w cieczy są dodatnie.

Ostatnią z wymienionych metod wyznaczania pozycji granicy międzyfazowej jest metoda śledzenia frontu (ang. front-tracking). Jej ogólna zasada polega na podziale obszaru, na którym są prowadzone obliczenia, na dwa niezależne od siebie regiony, a następnie na określeniu położenia granicy pomiędzy tymi fazami w każdym kroku czasowym. Największymi zaletami tej metody są: możliwość dokładnego określenia położenia granicy międzyfazowej oraz dokładnego określenia wydzielonego się ciepła krystalizacji. Natomiast wśród jej największych wad należy wymienić wysoki nakład obliczeń, jakie trzeba wykonać, żeby uzyskać wynik, oraz złożoność programistyczna, która się zwiększa wraz z liczbą wymiarów przestrzeni [92].

Metoda śledzenia frontu jest przykładem metody, dla której potrzebny jest zmienny krok przestrzenny siatki różnicowej, przynajmniej w okolicach interfejsu. W celu rozwiązania tego problemu na przestrzeni lat opublikowano szereg prac, w których stosuje się różne podejścia pozwalające uzyskać numeryczne rozwiązanie problemu Stefana. Jedną z pierwszych, a zarazem jedną z ciekawszych była praca przedstawiona przez Murraya i Landisa [25], w której autorzy zaproponowali dwie różne metody śledzenia frontu.

Pierwsza z nich opierała się na wykorzystaniu dwóch fikcyjnych poziomów temperatury, będących wynikiem jej ekstrapolacji odpowiednio w fazie stałej i ciekłej. Pozwalało to, w połączeniu z temperaturą topnienia i pozycją interfejsu znaną z poprzedniego kroku czasowego, na wyznaczenie temperatury w komórkach sąsiadujących z interfejsem. W celu uzyskania formuły opisującej ruch frontu wykorzystane zostało rozwinięcie w szereg Taylora. Druga z metod zakładała zastosowanie „elastycznej” siatki. Zarówno w obszarze fazy stałej, jak i fazy ciekłej, zastosowano siatki różnicowe o równomiernym kroku przestrzennym. W trakcie krystalizacji krok siatki w obszarze fazy stałej zwiększał się, a fazy ciekłej malał, natomiast liczba kroków obu siatek pozostawała bez zmian. Niestety, takie podejście powodowało, że należało przyjąć na początku obliczeń minimalną niezerową grubość warstwy fazy zakrzepłej, a obliczenia należało zakończyć, zanim skończy się krystalizacja. Każda zmiana kroku siatki wymagała ekstrapolacji w celu obliczenia nowych wartości węzłowych pola stężenia.

W pracy [93] została zaprezentowana jawna metoda modelowania krystalizacji dwufazowej z wykorzystaniem metody różnic skończonych zarówno dla dwóch jak i dla trzech wymiarów. Metoda ta została oparta na wcześniejszych pracach [94, 95] i wykorzystywała zestaw pomocniczych, bardziej skomplikowanych równań różniczkowych. Dla przypadku dwuwymiarowego był to zestaw 4 dodatkowych równań, natomiast dla przypadku trójwymiarowego było tych równań aż 6. Dla punktów siatki oddalonych od interfejsu stosowane było standardowe przybliżenie metodą różnic skończonych. Natomiast w pobliżu granicy międzyfazowej (gdzie odległość między węzłem siatki, a pozycją interfejsu zmienia się w każdym kroku czasowym) zastosowano wcześniej przygotowane równania pomocnicze. Wybór, które z równań pomocniczych ma zostać użyte, zależał od odległości między pozycją interfejsu a węzłem siatki. Pomimo dość dobrych wyników, jakie można było uzyskać dzięki temu modelowi, nie był on często wykorzystywany ze względu na złożoność równań pomocniczych.

Nowe podejście do metody śledzenia frontu na podstawie krystalizacji czystych substancji, krystalizujących dendrytycznie, z użyciem metody różnic skończonych oraz prostego algorytmu iteracyjnego zostało pokazane w pracy [96]. Proponowane podejście wykorzystywało dwie metody: interpolację wielomianową (interpolację Lagrange’a), która służyła do opisu poruszającego się interfejsu, oraz metodę Eulera, która została użyta do opisu pola temperatury na nieruchomej siatce. Wyniki obu metod zostały połączone ze sobą za pomocą tzw. techniki zanurzonej granicy (ang. immersed boundary technique). Technika ta jest dość skomplikowana i składa się z kilku etapów:

• obliczenie wektora prędkości normalnej dla każdego punktu leżącego na interfejsie, a następnie na jego podstawie wyznaczanie nowego położenia każdego z punktów w nowym kroku czasu;

• wyznaczanie funkcji charakterystycznej zbioru (ang. indicator function), której zadaniem jest opis nieciągłości właściwości materiałów (na granicy faza stała – ciecz);

• wraz z odpowiednimi warunkami brzegowymi wyznacza się rozkład pola temperatury na nieruchomej siatce;

• temperatura w każdym punkcie interfejsu jest interpolowana na podstawie wcześniej wyznaczonego pola temperatury, jak i nowego położenia frontu;

Jeżeli warunek nie został spełniony, wyznaczana jest nowa wartość prędkości interfejsu, po czym następuje powtórzenie poszczególnych kroków.

Niestety, ze względu na dość problematyczny sposób wyznaczania współrzędnych interfejsu w przestrzeni trójwymiarowej, model ten ograniczony został tylko do dwóch wymiarów. Metoda ta, oprócz modelowania wzrostu dendrytów, mogła również być wykorzystana do modelowania napięcia powierzchniowego oraz zmian topograficznych w strukturze materiałów.

Pomimo wielu trudności, związanych ze śledzeniem pozycji frontu dla stałego kroku siatki w przestrzeni trójwymiarowej, w pracy [97] został przedstawiony sposób postępowania w takich przypadkach. Zaproponowana metoda, wykorzystująca jawny schemat śledzenia pozycji frontu w przestrzeni trójwymiarowej, podzielona została na trzy główne części: śledzenie pozycji interfejsu, wyznaczanie zarówno powierzchni normalnej, jak i prędkości normalnej interfejsu oraz rozwiązywanie równania zachowania energii w komórkach

międzyfazowych. Główną ideą proponowanej metody jest określanie położenia specjalnie

wybranych punktów – „znaczników” (tzw. marker points) będących punktami przecięcia się interfejsu z pionowymi liniami siatki. Następnie dla tych punktów wyznaczana jest prędkość w kierunku normalnym do granicy. W kolejnym kroku czasowym znaczniki przemieszczają się, stając się tzw. punktami dodatkowymi (ang. advected points), które najczęściej nie leżą już na liniach siatki. Aby wyznaczyć współrzędne nowych znaczników w nowym kroku czasowym oraz poprawić położenie granicy, została zaproponowana tzw. funkcja sklejana (ang. cubic-spline reconstruction). Model ten zakłada nierównomierny pole temperatury w całej objętości, ale nie uwzględnia przepływu ciepła. Zostało to poprawione w pracy [98], w której przedstawiono również walidację tej metody. W pracy tej wyniki obliczeń numerycznych zostały porównane z rozwiązaniem analitycznym dla cylindrycznego i kartezjańskiego układu współrzędnych oraz z wynikami doświadczenia przeprowadzonego z wlewkiem aluminiowym. Otrzymane wyniki wskazały, że zaproponowany algorytm umożliwia bardziej dokładne śledzenie frontu krystalizacji.

Metoda śledzenia frontu jest metodą bardzo uniwersalną i pomimo pewnych trudności, związanych z modelowaniem zjawisk w przestrzeni trójwymiarowej, jest ona bardzo często proponowana przez autorów kolejnych publikacji. Oprócz śledzenia frontu krystalizacji metali i stopów, nadaje się również do opisu zjawisk, jakie zachodzą podczas przepływów cieczy [99], śledzenia ruchu pęcherzyków gazu w cieczy [100] i wielu innych.

2.

Cel i teza pracy

Znane metody numeryczne symulacji wzrostu ziaren równoosiowych (w tym dla krystalizacji jednofazowej, eutektycznej lub perytektycznej) dostarczające uśrednione dane, odnośnie do parametrów powstającej mikrostruktury, wykorzystują z reguły idealizowaną geometrię elementarnego pola mikrodyfuzji – EPMD (kulistą lub walcową).

Metody numeryczne, polegające na bezpośrednim modelowaniu mikrostruktury krystalizującego stopu, reprezentują geometrię ziaren i EPMD za pomocą wstępnej teselacji Dirichleta (przeprowadzanej przed wykonaniem symulacji). Mikrostruktura stopu w takich modelach jest reprezentowana zbiorami wielokątów lub wielościanów Voronoia. Wadą tych metod jest konieczność rozwiązywania jednokierunkowego zadania dyfuzyjnego w poszczególnych segmentach wszystkich komórek Voronoia.

Celem niniejszej pracy jest opracowanie modelu matematycznego w skali mikro, dostarczającego dane uśrednione, odnośnie do mikrostruktury i pola dyfuzji, pozwalającego na symulację wzrostu ziaren równoosiowych w warunkach krystalizacji perytektycznej i eutektycznej, oraz opracowanie algorytmu rozwiązywania równań modelowych i programu komputerowego pozwalającego na przeprowadzenie symulacji numerycznej.

Model uwzględniać będzie zarówno niejednorodność chemiczną zanikającej fazy macierzystej, jak i powstających produktów przemiany. Uwzględniony zostanie także efekt dyfuzji wstecznej. Dzięki temu model będzie mógł precyzyjnie wyznaczać niejednorodność powstających produktów przemian fazowych.

Kształt i związane z nim parametry geometryczne EPMD będą wyznaczane na podstawie statystycznej teorii krystalizacji Kolmogorova.

W oparciu o przedstawioną w poprzednim rozdziale analizę danych literaturowych zasadniczą tezę pracy można sformułować w następujący sposób:

Istnieje możliwość symulacji numerycznej uśrednionych parametrów wzrostu ziaren równoosiowych w strefie ciekło-stałej, limitowanych dyfuzją składników stopowych w warunkach krystalizacji perytektycznej lub eutektycznej w strukturach typu żeliwa z grafitem kulkowym, w warunkach natychmiastowego zarodkowania i wzrostu, z uwzględnieniem odchylenia geometrii elementarnego pola mikrodyfuzji od geometrii idealnej.

3.

Opis modelu

3.1. Aproksymacja ró

ż

nicowa geometrii EPMD – u

ś

redniony

wielobok Voronoia

W znanych modelach matematycznych wzrostu ziaren równoosiowych najczęściej wykorzystywana jest idealna geometria sferycznego pola dyfuzji. Niestety, jak już wspomniano wyżej, rzeczywisty kształt ziarna w krystalizującym stopie różni się od kulistego. W celu uwzględnienia odchylenia kształtu rzeczywistego pola dyfuzji w przypadku ziaren równoosiowych wykorzystano zasady podziału przestrzeni (teselacji) Dirichleta i podstawowe właściwości wielościanów Voronoia, wspomnianych w rozdziale 1.4.

W układzie z losowym rozkładem zarodków, poszczególne wielościany Voronoia będą zawierały wszystkie punkty znajdujące się bliżej ich środków (zarodków ziaren), niż środka jakiegokolwiek innego ziarna. Zewnętrzne ściany takich wielościanów są fragmentami płaszczyzn prostopadłych do odcinków łączących sąsiednie zarodki i dzielące ten odcinek na dwie równe części. Dokładny kształt i objętość każdego z tak powstałych wielościanów, jak również liczba powierzchni i krawędzi będą zależały od pozycji sąsiednich zarodków. Struktura tego typu powstaje podczas jednoczesnego zarodkowania ziaren i ich sferycznego wzrostu, z jednakową (ale niekoniecznie stałą) prędkością. Schematyczne przedstawienie teselacji Voronoia na płaszczyźnie 2D na tle mikrostruktury żeliwa z grafitem kulkowym zostało zaprezentowane na rys.3.1.

W proponowanym modelu zakłada się, że wewnątrz obszaru EPMD, w dowolnym punkcie i w dowolnym momencie czasu stężenie danego pierwiastka C może być opisane jako funkcja odległości od środka EPMD (zarodka) r i czasu τ.

( ) ( )

r,τ = f r,τ

Rys. 3.1. Schemat wielościanów Voronoia dla płaszczyzny 2D na tle mikrostruktury żeliwa z grafitem kulkowym [SPCI 10, oddane do druku]

Obszar geometryczny EPMD, wykorzystany w niniejszym modelu, został nazwany

Uśrednionym Wielościanem Voronoia (UWV).

Właściwości topologiczne wielościanów Voronoia zostały przedstawione wcześniej w rozdziale Diagramy Voronoia 1.4. Zależności matematyczne i cechy tej bryły, niezbędne do rozwiązania problemu jednowymiarowego zjawiska mikrodyfuzji składników w EPMD, wyznaczono na podstawie statystycznej teorii krystalizacji Kolmogorova [101]. Zgodnie z tą teorią, w przypadku natychmiastowego zarodkowania n ziaren w jednostce objętości i ich sferycznego wzrostu z jednakową (ale niekoniecznie stałą) prędkością, objętość materiału, znajdującego się w odległości nie przekraczającej r od miejsca zarodkowania, można wyznaczyć za pomocą następującej zależności:

( )

            π − − = 3 3 4 exp 1 1 nr n r V (3.2) gdzie:

n – liczba ziaren przypadająca na jednostkę objętości.

Funkcja ta opisuje średnią objętość, która została nazwana jako objętość wewnętrzna UWV.

Pole powierzchni, oddzielającej wewnętrzną część UWV V

( )

r od jego części zewnętrznej, bardziej oddalonej od miejsca zarodkowania, może zostać opisane jako:

( )

( )

     _{− π} ⋅ π = = 2 3 3 4 exp 4 r nr dr r dV r F (3.3)

Wartość otrzymana z równania (3.3) została nazwana jako powierzchnia rozdzielająca

UWV. Dla uśrednionej komórki Voronoia, prawdopodobieństwo tego, że dowolnie wybrany punkt tego ziarna będzie się znajdował w odległości mniejszej niż r od środka ziarna opisuje

dystrybuanta:

( )

      π − − = 3 3 4 exp 1 nr r P (3.4)

Gęstość prawdopodobieństwa tego rozkładu statystycznego:

( )

     _{− π} ⋅ π = 2 3 3 4 exp 4 nr nr r p (3.5) ma maksimum w punkcie:

( )

13 2π − = n R_VP (3.6)

Parametr R może zostać nazwany promieniem charakterystycznym UWV. Można go VP

scharakteryzować w następujący sposób: spośród wszystkich losowo wybranych punktów przestrzeni, dla której przeprowadzono teselację Voronoia, najwięcej punktów będzie znajdować się w odległości RVP od najbliższego środka wielościanu. Powierzchnia

rozdzielająca o tym promieniu F(RVP)będzie miała również największą wartość.

Zgodnie z zależnością (3.2), dla małego promienia r objętość wewnętrzna UWV

zbliżona jest do objętości kuli o tym samym promieniu. Wraz ze wzrostem promienia ziarna

r, objętość kuli o tym promieniu będzie rosła do nieskończoności, natomiast objętość

wewnętrzna UWV będzie dążyć do 1/n (rys.3.2a).

Tak samo, jak w przypadku objętości wewnętrznej, dla małej wartości promienia r,

W zg lę d n a w e w n ę tr zn a o b ję to ść W zg lę d n a p o w ie rz ch n ia o d d zi e la jąc a

Względny promień EPMD R = r/RVP Względny promień EPMD R = r/RVP

a) b)

Rys. 3.2. Względna objętość wewnętrzna (a) i względna powierzchnia oddzielająca w UWV, dla kształtu sferycznego (linia przerywana) i UWV (linia ciągła)

w funkcji względnego promienia [102, 103]

3.2. Numeryczna metoda rozwi

ą

zywania

W celu rozwiązania równań opisujących dyfuzję składnika stopowego w krystalizującym materiale, została wykorzystana metoda bilansów elementarnych. Schemat wytyczania objętości kontrolnych, użyty w proponowanym rozwiązaniu z wykorzystaniem UWV, wraz z uwzględnieniem zastosowanej notacji zmiennych został pokazany na rys. 3.3a. Dla porównania, na rys. 3.3b zamieszczony został schemat dla pola dyfuzji o idealizowanym sferycznym kształcie pola dyfuzji.

Rys. 3.3. Schemat EPMD bazującego na Uśrednionym Wielościanie Voronoia (linie przerywane wskazują granice ziaren z zerowym strumieniem dyfuzyjnym) (a) i schemat sferycznego EPMD (b)

3.2.1. Parametry EPMD

Odpowiedni wybór wymiaru EPMD w przypadku założenia jego idealnego kształtu kulistego nie jest oczywisty. Przykłady sposobów wyznaczenia promienia stosowanych w znanych modelach przedstawiono w tabeli 3.1.

Wybór wielkości promienia obliczonego na podstawie bilansu objętości (poz. 1 w tab. 3.1) oznacza, że w przypadku losowego (nieregularnego) rozmieszczenia środków tych ziaren, w przestrzeni obszary zewnętrzne sąsiednich pól będą nakładały się na siebie. Równocześnie, zgodnie z równaniem (3.2), ponad 30% objętości próbki będzie znajdowało się poza obszarami EPMD. Zwiększenie promienia EPMD (na przykład według wzoru poz. 2 w tab. 3.1) skutkuje zaburzeniem bilansu objętości i masy pierwiastków w polu dyfuzji.

Tab. 3.1 Promienie sferycznego EPMD

L.p. Promień sferycznej

elementarnej komórki Wyjaśnienie oznaczeń Źródło

1. 3 1 4 3       π = N r_o N - gęstość ziaren, 1/m3 [30, 104, 105] 2. 3 1 0 4 3       = N_β r Nβ – gęstość ziaren, 1/m3 [47] 3. lRVE Ng d 3 1

= d −średnia średnica (wymiar) ziarna Ng – liczba ziaren, 1/m3 [106] 4. 0 ₃ 1 N r = N – gęstość ziaren, 1/m3 [49] 5. 3 1 3 4 3       π = N l r ₃ − l N

końcowa gęstość ziaren [42]

6. A L N N R π

= 2 NA – gęstość powierzchniowa ziaren,

NL – średnia liczba przecięć z prostą, [40]

7. 3 1 0 0 4 3       π = N V

R N0 – liczba komórek Voronoia przypadająca

na objętość V0, [65] 8. 3 4 3 n f R e π = n f_e

     _{− π} = ε 3 3 4 exp nr (3.7)

Wartość tę możemy nazwać stałą zaokrąglenia. Zakładając małą wartość stałej

zaokrąglenia, czyli poziom tego ułamka objętości, która nie będzie uwzględniona w modelowaniu (zwykle rzędu 10-3), promień EPMD możemy zdefiniować następująco:

3 4 ln 3 n R_EPMD π ε − = (3.8)

Zgodnie z (3.2), pomiędzy środkiem ziarna a promieniem charakterystycznymR _VP

znajdować się będzie 48,7% objętości danego ziarna, podczas gdy w obszarze znajdującym się w odległości większej niż podwójny promień charakterystycznym 2⋅R_VP będzie występować mniej niż 0.5% objętości danego ziarna.

Objętość stosowanego w rozwiązaniu numerycznym elementu kontrolnego, pokazanego na rys. 3.3a zgodnie z równaniem (3.2), wynosi:

           _{− π} −      _{− π} = ∆ 3₊ 1 3 3 4 exp 3 4 exp 1 i i i nr nr n V (3.9)

Dla stałego kroku siatki ∆r=r_i+1−r_i, wraz ze wzrostem wartości r od 0 do RVP,

wartość ∆V również wzrasta. Powyżej wartości RVP, na skutek kontaktów z sąsiednimi

ziarnami, wartość ∆V zmniejsza się i dąży do zera, pomimo dalszego wzrostu promienia.

3.2.2. Transport masy

Zjawisko dyfuzji opisane jest za pomocą dwóch praw Ficka. W proponowanym modelu, ze względu na zmianę strumienia dyfuzji w czasie, wykorzystuje się drugie prawo Ficka w postaci:

( )

₌

_[

_⋅

_{( )}

_τ

_]

τ ∂ τ ∂ , grad div , x C D x C f (3.10) gdzie: