MODEL HARTOWANIA ELEMENTU OSIOWO-SYMETRYCZNEGO Z UWZGLĘDNIENIEM RUCHÓW CHŁODZIWA A

(1)

41, s. 213-220, Gliwice 2011

MODEL HARTOWANIA ELEMENTU OSIOWO-SYMETRYCZNEGO Z UWZGLĘDNIENIEM RUCHÓW CHŁODZIWA

ADAM KULAWIK

Instytut Informatyki Teoretycznej i Stosowanej, Politechnika Częstochowska e-mail:adam.kulawik@icis.pcz.pl

Streszczenie. W pracy zaproponowano model matematyczny oraz numeryczny procesu chłodzenia zbudowany z wykorzystaniem uogólnionej metody różnic skończonych. Ruchy chłodziwa, zarówno konwekcyjne jak i wymuszone, zamodelowano, rozwiązując równanie Naviera-Stokesa metodą rzutowania.

Rozwiązanie równania przewodzenia ciepła z członem konwekcyjnym uzyskano na podstawie stabilizowanej bezsiatkowej metody różnic skończonych. Do wyznaczenia udziałów przemian fazowych zastosowano makroskopowy model bazujący na analizie wykresów CTPc.

1. WSTĘP

Zjawiska towarzyszące hartowaniu są trudne do modelowania numerycznego z powodu ich złożoności oraz wzajemnych ich powiązań. Do najważniejszych z nich zalicza się zjawiska cieplne, mechaniczne i przemiany fazowe. W pracy prezentowane są modele matematyczne oraz numeryczne zjawisk towarzyszących stygnięciu elementów stalowych, a w szczególności rozbudowany model zjawisk cieplnych oraz model przemian fazowych w stanie stałym.

W modelowaniu zjawisk cieplnych najczęściej przyjmuje się, że warunki chłodzenia przybliżane są za pomocą odpowiednich warunków brzegowych, które w większym lub mniejszym stopniu przybliżają warunki rzeczywiste procesu. Jednak takie podejście, zwłaszcza w modelowaniu procesów obróbki cieplnej, może być obarczone dużym błędem przybliżenia. Ma to miejsce wówczas, gdy chłodziwem są ciecze, a także dla skomplikowanych geometrycznie części maszyn. Nie można w takich przypadkach przyjmować warunków chłodzenia identycznych dla całego brzegu, pomijać w modelowaniu wpływ grawitacji oraz wymuszony ruch chłodziwa. Dlatego też model zjawisk cieplnych dla procesów obróbki termicznej powinien uwzględniać ruchy cieczy chłodzących.

W modelowaniu zjawisk cieplnych chłodzenia oprócz ruchów cieczy należy także uwzględnić przemiany fazowe w stanie stałym. Przemiany te wpływają na pola temperatury poprzez ciepła transformacji oraz mają znaczący wpływ na pola naprężeń chwilowych i własnych.

(2)

2. MODELOWANIE PROCESU HARTOWANIA

Obróbka termiczna elementu stalowego jest procesem łączącym zjawiska cieplne, przemiany fazowe i zjawiska mechaniczne (rys. 1). Poprawnie zbudowane modele hartowania powinny uwzględniać te zjawiska oraz związki pomiędzy nimi.

Rys.1. Schemat ideowy modelu hartowania

W pracy przedstawiono modele matematyczne oraz numeryczne procesu hartowania cieczą osiowo-symetrycznego elementu stalowego.

2.1. Modelowanie zjawisk cieplnych

W prezentowanym modelu do wyznaczania pól temperatury wykorzystano równanie różniczkowe przewodzenia ciepła z członem konwekcyjnym we współrzędnych osiowo-symetrycznych w postaci

V z

r Q

t C T z V T r V T z C

T z r T r r T

r =

∂

− ∂

⎟⎠

⎜ ⎞

⎝

⎛

∂ + ∂

∂

− ∂

⎟⎠

⎜ ⎞

⎝

⎛

∂

∂ + ∂

⎟⎠

⎜ ⎞

⎝

⎛

∂

∂ + ∂

∂

∂ λ λ ρ ρ

λ (1)

gdzie: T [K] - temperatura, t [s] - czas, λ = λ(T) [W/mK] – współczynnik przewodzenia ciepła, ρ [kg/m³] - gęstość, C [J/kgK] - ciepło właściwe, V [m/s] – wektor prędkości, QV [W/m³] - objętościowe źródło ciepła.

Równanie (1) uzupełniono warunkami brzegowymi Dirichleta, Neumanna oraz Newtona- Robina, a także odpowiednimi warunkami początkowymi.

Do rozwiązania równania (1) wykorzystano ugólnioną metodę różnic skończonych (GFDM) w schemacie niejawnym całkowania po czasie [4,5]. Macierz współczynników wyznacza się z zależności.

( ) ( ) ( ) ( )

( )

j

(

j

( )

r i j

( )

z i

)

i j j j i

n

j

i z j n

j

i r j n

j j i n

j

j j i

i

V s V s C r s

s s Am

V s V

s C r s

s s Am

) 2 ( )

1 ( ) 1 ( ) 4 ( ) 3 ( ,

1 ) 2 ( 1

) 1 ( 1

) 4 ( ) 3 ( ,

1

+

− + +

=

⎟⎟⎠⎞

⎜⎜⎝⎛

+

⎟⎟⎠+

⎜⎜⎝ ⎞

⎛− + −

=

∑ ∑ ∑ ∑

=

λ ρ λ

ρ λ

(2)

Wektor wyrazów wolnych wyznaczany jest następująco

Vi s i

i T Q

t

D C −

−Δ

= ρ ⁻¹ (3)

gdzie

(3)

⎟⎠

⎜ ⎞

⎝

⎛ + + + +

= m i j i j i j i j i j j

j i

j g h g k g h g k g hk

s⁽⁾ l₂ ₁ ₂ ₃ ² ₄ ² ₅

2 1 2

1

1 (4)

Współczynniki gij,to elementyodwróconej macierzy [4,5]

1

1 2

2 2

1 2

3

1 2 3

1 2

2

1 2 2

1 2

3

1 2 4

1 2

2 2

1 2 3

1 2

2

1 2 3

1 2

2 2

1 2 4

1 2 2

1 2 3

1 2

2

1 2 3

1 2 2

1 2

2

1 2

1 2 2

1 2

2

1 2 3

1 2 1

2 2

2 4

4 2

2

2 4

4 2

2

2 2

−

=

⎥⎥

⎦

⎤

⎢⎢

⎣

⎡

=

∑

n

j m j

j n j

j m j

j n j

j m j

j n j

j m j

j n j

j m j

j j

n

j m

j j n j

j m

j n j

j m

j j n j

j m

j n j

j m

j j j

n

j m

j j n j

j m

j j n j

j m

j n j

j m

j j n j

j m

j j

n

j m j

j n j

j m j n j

j m j

j n j

j m j n j

j m j

j j

n

j m

j j n j

j m

j j n j

j m

j n j

j m

j j n j

j m

j j

ij

l k h l

k h l

k h

l k h l

k l

k h l

k l

k h

l k h l

k h l

h l

k h l

h

l k h l

k l

k h l

k l

k h

l k h l

k h l

h l

k h l

h

g (5)

gdzie: h_j =r_j−r_i, k_j=z_j−z_i, lj=

(

h²j +k²j

)

, „i” oznacza węzeł wewnętrzny siatki węzłowej, „j” oznacza węzły otaczające węzeł „i”, parametr m przyjmuje się w zależności od liczby węzłów siatki węzłowej.

Do modelowania przepływu chłodziwa wykorzystano równanie Naviera-Stokesa w postaci

( )

^⎟

⎠

⎜ ⎞

⎝

⎛

∂ + ∂

=

−

∂ +

−∂

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂ +∂

∂ + ∂

∂

⎟⎠

⎜ ⎞

⎝

⎛

∂ + ∂

∂ =

−∂

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ −

∂ +∂

∂ + ∂

∂

z V V r V V dt T dV

T z a

p z

V r V r r

V

z V V r V V dt dV r

p r V z

V r V r r

V

z z z r z ref

z z

r z r r r r

r r r

ρ β

ρ μ

2 2 2

2

2 2 2 2

2

1 1

(6)

gdzie: V [m/s] - wektor prędkości, μ [kg/ms] - lepkość dynamiczna, az [m/s²] - przyspieszenie, β [K^-1] - objętościowy współczynnik rozszerzalności cieplnej, Tref [K] - temperatura odniesienia.

Równanie (6) uzupełniono równaniem ciągłości

=0

∂ +∂

∂ +

∂

z V r V r

V_r _r _z

(7)

Do rozwiązania równań (6) oraz (7) wykorzystano metodę rzutowania (CBS) [3,7,11], w której wyznacza się prędkości pomocnicze rozwiązując równanie pędu w postaci

( )

¹ ^*

2 2 2

2

* 1 2

2 2 2

2

1 1

z s z ref z

z z z

r s r r r r r

V V T T z a

V r V r r

t V

V r V

V z

V r V r r

V t

=

⎟+

⎟

⎠

⎞

⎜⎜

⎝

⎛ ⎟⎟⎠+ −

⎜⎜ ⎞

⎝

⎛

∂ +∂

∂ + ∂

∂

∂ Δ

=

⎟+

⎟

⎠

⎞

⎜⎜

⎝

⎛ ⎟⎟⎠

⎜⎜ ⎞

⎝

⎛ −

∂ +∂

∂ + ∂

∂ Δ ∂

−

β ρ ρ μ

ρ μ

(8)

Korektę wektora prędkości przeprowadza się na podstawie rozwiązania równania w postaci

(4)

α α

α ρ ^,

* t p

V

V − =Δ (9)

gdzie wartość ciśnienia otrzymuje się z rozwiązania równania

⎟⎟⎠

⎜⎜ ⎞

⎝

⎛

∂ +∂

∂ + ∂

∂

=Δ

∂ +∂

∂ +

∂

2 2 2

* 2

*

* 1

z p r p r r

p t z V r V r

V_r _r _z

ρ (10)

Ostatecznie wartość prędkości otrzymuje się z sumowania prędkości pomocniczych oraz poprawki wyznaczonej na podstawie ciśnień.

( )

α α

α ^, ρ ^,

*

* t p

V V V

V = +Δ Δ =−Δ (11)

Macierz współczynników dla rozwiązania równania pędu za pomocą GFDM ma postać

( )

_⎟⎟⎠^⎞

⎜⎜⎝⎛

− + Δ −

=

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

⎟⎟⎠⎞

⎜⎜⎝⎛

+ +

− Δ

−

=

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛

⎟⎟⎠⎞

⎜⎜⎝⎛

+ +

−

− Δ

−

=

∑

=

) 1 ( ) 4 ( ) 3 ( ,

1 ) 1 ( 1

) 4 ( ) 3 ( ,

1 ) 1 ( 1

) 4 ( ) 3 ( , 2

1 1 1

j i j j j

i

n

j j i n

j

j j y

i i

n

j j i n

j

j j i

x i i

rs s t s

Am

r s s s t

Am

r s s r s

t Am

ρμ ρ μ ρ μ

α

(12)

Wektor wyrazów wolnych ma postać

( ) ( )

_⎟^⎟

⎠

⎞

⎜⎜

⎝

⎛ ⎟⎟

⎠

⎜⎜ ⎞

⎝

⎛ −

+

− Δ

−

= ^s⁻ ⁱ ⁱ ^∞

∑

_jⁿ₌ ^j ^s^j⁻

∑

_jⁿ₌ ^j ^sⁱ⁻

i V t a T T z V z V

Dw

1 1

1 ) ( 1 ) ( 1

γ α γ α

α

α α β (13)

Rozwiązanie równania (10) za pomocą GFDM w schemacie niejawnym ma postać

( ) ( )

⎟⎟

⎠

⎞

⎜⎜

⎝

⎛ + +

=

⎟⎟⎠⎞

⎜⎜⎝⎛

− +

−

=

∑

₌

∑

₌

i j j j j i

n

j j i n

j

j j i

i

r s s s Am

) 1 ( ) 4 ( ) 3 ( ,

1 ) 1 ( 1

) 4 ( ) 3 ( ,

1

1 1

ρ

ρ (14)

( ) ( ) ( ) ( ) ( )

_⎟⎟⎠^⎞

⎜⎜⎝⎛

+

− +

Δ −

=

∑

₌

∑

₌

∑

₌

∑

₌

i r i n

j z i j n

j z j j n

j r i j n

j r j j

i r

V V s V

s V

t s D

1

* ) 2 ( 1

* ) 1 ( 1

* ) 1

1 (

(15)

Wartość korekty wektora prędkości wyznaczona jest zależnością

⎟⎟⎠⎞

⎜⎜⎝⎛ Δ −

−

=

Δ^V^α^*ⁱ _ρ^t

∑

_jⁿ₌₁^s⁽^j^α⁾^p^α^j

∑

_jⁿ₌₁^s⁽^j^α⁾^p^αⁱ ⁽¹⁶⁾

(5)

Rozwiązania równań pędu oraz przewodzenia ciepła stabilizowano za pomocą siatek kierunkowych [9]. Prezentowane układy równań rozwiązano z wykorzystaniem stabilizowanej metody sprzężonych gradientów [2].

3. PRZEMIANY FAZOWE

W prezentowanym modelu przemian fazowych wykorzystany został wykres ciągłego chłodzenia CTPc stali C45 [10]. Z diagramu tego odczytano chwilowe wartości czasu, temperatury oraz maksymalne udziały faz dla danego cyklu chłodzenia. Kinetykę faz określano na podstawie empirycznego równania Avramiego dla procesu chłodzenia [1,8]:

( ) ( ( ( )

ⁿ^{( )}^T

) )

i j

j i

i T t ⋅ − −bTt

⎭⎬

⎫

⎩⎨

⎧ −

=^min ^,^~

∑

_≠ ¹ ^exp

, ₍_%)

)

( η η η

η _γ (17)

gdzie η^~_γ jest objętościowym udziałem austenitu, η jest objętościowym udziałem faz _j powstających podczas procesu chłodzenia, η_{( )}_i_% jest oszacowanym na podstawie wykresu CTPc końcowym udziałem fazy, funkcje n(.) i b(.) określają równania [8]

( ) ( )

n_{( )}T s

f bT

T t

T T t

n

ts

0,01005 ,

ln

6,12733/ =

⎟⎟⎠⎞

⎜⎜⎝⎛

= (18)

Udział powstałego martenzytu wyznacza się na podstawie rozszerzenia empirycznego równania Koistinena i Marburgera [6]

( )

^, ^~ ⎟

(

¹−^exp

(

−

(

−

) ) )

^, =^0.01537

⎠

⎜ ⎞

⎝

⎛ −

=

∑

_≠ ^k ^M ^T ^k

t

T _S

M i

i

M η η

η γ (19)

Podczas przemian fazowych zachodzą widoczne zmiany objętości obrabianego cieplnie materiału. Zjawisko to jest następstwem odkształceń termicznych i strukturalnych. Wyznacza się je wzorami

( ) ( )

i

i ph i ph i

i i

T T dT dε T d

dε ⁼

∑

α η ^, ⁼

∑

ε η ⁽²⁰⁾

gdzie: αi

( )

^T są współczynnikami rozszerzalności termicznej poszczególnych faz, εi^ph

( )

^T są współczynnikami zmian objętości od przemian fazowych odpowiednio: austenitu w bainit, austenitu w ferryt, austenitu w martenzyt oraz austenitu w perlit.

4. PRZYKŁAD OBLICZEŃ

Przeprowadzono symulację numeryczną hartowania elementu w kształcie walca wykonanego ze stali C45 o wymiarach r=0.04 m, h=0.05 m. Element stalowy był wstępnie nagrzany do temperatury 1200 K, a następnie umieszczony w pionowym kanale chłodzącym wypełnionym chłodziwem o parametrach ciekłego sodu (rys. 2b). Założono, że prędkość chłodziwa na wlocie jest równa Vz=0,05 m/s. Wstępne prędkości wynikające z zanurzania zostały pominięte. Symulację procesu obróbki termicznej zakończono po wyrównaniu się temperatur elementu hartowanego oraz chłodziwa. Badano wpływ ruchów wymuszonych oraz konwekcyjnych chłodziwa na uzyskiwaną strukturę hartowanej części.

(6)

W prezentowanym przykładzie znaczącą rolę w procesie chłodzenia odgrywają ruchy konwekcyjne cieczy, które z powodu przyjętych założeń dla przepływu są sumowane razem z ruchami wymuszonymi, co prezentuje rys. 2b. Rola ruchów konwekcyjnych jest znacząca zwłaszcza w początkowym okresie procesu. Po kilku sekundach proces się normuje i główną rolę w chłodzeniu odgrywa ruch wymuszony cieczy.

0 0.02 0.04 0.06 0.08 0.1 r [m]

0 0.05 0.1 0.15 0.2

z [m]

250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200

0 0.025 0.05 0.075

r [m]

0 0.03 0.06 0.09 0.12 0.15

z [m]

V [m/s]

0.24

Rys.2. Czas procesu 8 s: a) pole temperatury [K], b) pole prędkości chłodziwa [m/s]

Wyniki obliczeń udziałów fazowych po zakończeniu procesu chłodzenia przedstawiono na rys. 3. Zauważalny jest wpływ ruchów chłodziwa na zaleganie poszczególnych faz. Można też zauważyć znaczący udział przemiany perlitycznej oraz ferrytycznej w miejscu o najmniejszej prędkości chłodzenia (rys. 2 i 3).

Na podstawie uzyskanych wyników można stwierdzić, że proces hartowania elementu stalowego jest bardzo złożony i trudny do modelowania za pomocą zazwyczaj przyjmowanych uproszczeń. Modelowanie obróbki termicznej elementów stalowych chłodzonych cieczami nie powinno być realizowane bez uwzględniania ruchów cieczy z powodu złożoności zjawiska. Nawet dla tak prostych geometrii (jak np. wałek) ruch cieczy powoduje bardzo złożony stan warunków termicznych na brzegach obrabianych części.

a) b)

(7)

0 0.01 0.02 0.03 0.04 r [m]

0.05 0.06 0.07 0.08 0.09 0.1

z [m]

0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 0.016 0.018 0.02 0.022 0.024 0.026 0.028 0.03 0.032 0.034 0.036

0 0.01 0.02 0.03 0.04

r [m]

0.05 0.06 0.07 0.08 0.09 0.1

z [m]

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65

0 0.01 0.02 0.03 0.04

r [m]

0.05 0.06 0.07 0.08 0.09 0.1

z [m]

0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 0.22 0.24 0.26 0.28 0.3 0.32 0.34 0.36 0.38 0.4

0 0.01 0.02 0.03 0.04

r [m]

0.05 0.06 0.07 0.08 0.09 0.1

z [m]

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1

Rys.3. Pola udziałów faz po zakończeniu procesu chłodzenia: a) ferryt, b) perlit, c) bainit, d) martenzyt

5. WNIOSKI

Zaproponowane modele umożliwiają przeprowadzenie symulacji hartowania elementów osiowo-symetrycznych. Uwzględnienie prędkości chłodziwa wpływa w dużym stopniu na dokładność oszacowania prędkości chłodzenia. Uwzględnienie prędkości chłodziwa pozwala na określenie warunków chłodzenia elementów o skomplikowanych kształtach, a zastosowanie złożonych modeli hartowania - w znacznym stopniu - na dokładniejsze modelowanie warunków rzeczywistych procesów. Zastosowana metoda numeryczna (GFDM) pozwala w łatwy sposób na zastosowanie zmiennej w czasie siatki węzłowej. Zastosowana metoda stabilizacji umożliwia modelowanie zjawisk z dużymi wartościami prędkości.

a) b)

c) d)

(8)

LITERATURA

1. Avrami M.: Kinetics of phase change. „J. Chem. Phys.”1939, Vol. 7, p. 1103.

2. Axelsson O.: Iterative solution methods. Cambridge University Press 2000.

3. Chorin A. J.: Numerical solution of the Navier-Stokes equation. „Math. Comput.” 1968, 23, p. 745-762.

4. Benitoa J.J., Ureńaa F., Gaveteb L.: Solving parabolic and hyperbolic equations by the generalized finite difference method. „Journal of Computational and Applied Mathematics” 2007, 209, p. 208-233.

5. Liszka T.: An interpolation method for an irregular net of nodes. „Int. J. Numer. Methods Eng.”1984, 20, p. 1599–1612.

6. Koistinen D. P., Marburger R. E.: A general equation prescribing the extent of the autenite-martensite transformation in pure iron-carbon alloys and plain carbon steels.

„Acta Metallica” 1959, 7, p. 59-60.

7. Kulawik A., Skrzypczak T., Węgrzyn-Skrzypczak E.: Influence of coolant motion on structure of hardened steel element. “Archives of Foundry Engineering” 2008, 8, p. 39- 44.

8. Kulawik A.: Analiza numeryczna zjawisk cieplnych i mechanicznych w procesach hartowania stali 45. Praca doktorska. Częstochowa 2005.

9. Kulawik A.: Modelling of forced convection using a stabilized meshless metod.

„Scientific Research of the Institute of Mathematics” 2010, 2(9), p. 123-130.

10. Wever F., Rose A.: Atlas zur Wärmebehandlung von Stähle, I Zeit Temperatur Umwandlungs Schaubilder. Düsseldorf: Verlag Stahl Eisen MBH, 1961.

11. Zienkiewicz O. C., Codina R.: A general algorithm for compressible and incompressible flow. Part I: The split characteristic based scheme. „International Journal for Numerical Metods in Fluids”1995, 20, p. 869-885.

MODEL OF HARDENING PROCESS OF THE AXISYMMETRIC ELEMENT WITH THE MOVEMENT OF COOLANT

Summary. In the paper a mathematical and numerical model for the cooling processes built on the basis of the generalized finite difference method have been proposed. The movements of the coolant, both convective and forced, have been solved using the Navier-Stokes equation with the characteristic based split scheme (CBS). The solution of the heat transport equation with the convective term has been obtained by a stabilized meshless finite difference method. In the paper the macroscopic model of the phase transformations based on an analysis of CCT diagrams has been presented.