NUMERYCZNY MODEL PRZEWIDYWANIA NAPRĘŻEŃ W HARTOWANYCH ELEMENTACH ZE STALI NARZĘDZIOWEJ

(1)

NUMERYCZNY MODEL PRZEWIDYWANIA NAPRĘŻEŃ W HARTOWANYCH

ELEMENTACH ZE STALI NARZĘDZIOWEJ

Adam Bokota

^1a

, Tomasz Domański

^2b

, Alieksiey Mikhaylovich Guriev

^3c

, Andrey Mikhaylovich Markov

^4d

1Instytut Informatyki Teoretycznej i Stosowanej, Politechnika Częstochowska

2Instytut Mechaniki i Podstaw Konstrukcji Maszyn, Politechnika Częstochowska

3Head of Department, I.I. Polzunov Altai State Technical University, Barnaul, Russia

4Department of Management of Technology, I.I. Polzunov Altai State Technical University, Bar- naul, Russia

abokota@icis.pcz.pl, ^bdomanski@imipkm.pcz.pl, ^cgurievam@mail.ru, ^dandmarkov@inbox.ru

Streszczenie

W pracy przedstawiono modele zjawisk mechanicznych procesu hartowania stali. Jako priorytetowe przyjęto zjawiska mechaniczne. Równania równowagi i związki konstytutywne wykorzystano w formie prędkościowej. Od- kształcenia plastyczne determinuje stowarzyszone prawo nieizotermicznego plastycznego płynięcia z warunkiem plastyczności Hubera-Misesa. Założono, że umocnienie materiału może być izotropowe i/lub kinematyczne. Oprócz odkształceń cieplnych, strukturalnych i odkształceń plastycznych uwzględniono również odkształcenia transformacyjne. Wielkości termofizyczne, takie jak moduł Younga, moduł styczny i granicę plastyczności uzależniono od temperatury i składu fazowego materiału. Zagadnienie termosprężystoplastyczności rozwiązano metodą elementów skończonych. Wykorzystując zaimplementowane algorytmy, wykonano symulacje hartowania elementu ze stali na- rzędziowej niestopowej. Dokonano analizy numerycznej wpływu efektów termicznych, przemian fazowych i rodzaju umocnienia materiału na naprężenia i odkształcenia w hartowanym elemencie ze stali narzędziowej.

Słowa kluczowe: hartowanie, przemiany fazowe, naprężenia, odkształcenia, modelowanie

NUMERICAL MODEL FOR PREDICTION OF STRESS IN HARDENED ELEMENTS MADE OF TOOL STEEL

Summary

The paper presents models of mechanical phenomena in steel hardening process. In the modeling as a priority was given to mechanical phenomena. Equilibrium equations and constitutive relationships were adopted in the rate form. Plastic strain determines the associated law of non-isothermal plastic flow with Huber-Misses yield condition. The model assumes that the strengthening of the material can be isotropic and/or kinematic. Besides the thermal strains, structural and plastic strains transformation plasticity was also taken into account.

Thermophysical values, such as Young's modulus and tangent modulus of the yield point were made dependent from the temperature and phase composition of the material. The issue of thermo-elastic-plasticity was solved by finite element method. Using implemented algorithms the simulations of hardening an element made of tool steel was performed. The effect of thermal phenomena, phase transformation and type of hardening of the material, i.e.

stresses and strains in the hardened element made of tool steel, was analysed.

Keywords: hardening, phase transformations, stresses, strains, modeling

(2)

Adam Bokota, Tomasz Domański, Alieksiey Mikhaylovich Guriev, Andrey Mikhaylovich Markov

1. WSTĘP

Różnorodność warunków pracy narzędzi stwarza ko- nieczność zróżnicowania wymagań stawianych stalom narzędziowym. Dążąc do osiągnięcia największej trwało- ści narzędzia, stosuje się odpowiednią obróbkę cieplną z zabiegiem hartowania.

Obecnie to metody numeryczne są wiodące w modelowaniu procesów technologicznych. Dają one możliwość analizy rozważanego zagadnienia w obszarach o skom- plikowanych kształtach, z dowolnymi warunkami po- czątkowymi i brzegowymi. Ponadto modelowanie numeryczne cechuje duża elastyczność zarówno z punktu widzenia rozbudowy algorytmu numerycznego, jak również możliwości uwzględniania zmiennych parame- trów termofizycznych analizowanych zjawisk.

Aby można było przeprowadzić pełną analizę obrób- ki cieplnej, konieczne jest dysponowanie odpowiednimi modelami matematycznymi i numerycznymi, mogącymi dawać informacje o chwilowych polach temperatury, zmianie w czasie udziałów poszczególnych składników fazowych materiału, rozkładach naprężeń chwilowych, a w konsekwencji naprężeń własnych.

Analiza numeryczna procesów obróbki cieplnej jest ważnym problemem stojącym przed współczesnymi pracowniami zajmującymi się projektowaniem dla prze- mysłu, niekoniecznie stalowego. Zjawiska towarzyszące obróbce cieplnej są złożone i do tej pory niekompletnie opisane. Szczególny nacisk na rozwój gałęzi metod numerycznych jest inspirowany przez przemysł, który ze względu na nowoczesne technologie i dążenie do ograni- czenia kosztów wymaga narzędzi usprawniających procesy obróbki cieplnej.

Określenie ilościowych udziałów poszczególnych faz, ich kinetyki oraz rodzaju uzyskiwanej struktury w procesie chłodzenia stopów żelaza jest nieodzowne do wyznaczania naprężeń towarzyszących takiemu zabiego- wi. Obecnie prace nad modelami numerycznymi obróbki cieplnej idą w kierunku uwzględnienia coraz większej liczby parametrów wejściowych takiego procesu [3,8].

Elementem mającym decydujący wpływ na wyniki symulacji numerycznej hartowania jest również właściwy dobór warunków chłodzenia, które modeluje się warunkami brzegowymi.

W procesie obróbki cieplnej generują się znaczące naprężenia. Aby zapewnić wiarygodność wyników symulacji numerycznych zjawisk mechanicznych, oprócz odkształceń termicznych, strukturalnych i plastycznych, uwzględniono również odkształcenia transformacyjne [5,9].

Do modelowania numerycznego zjawisk obróbki cieplnej wykorzystuje się metodę elementów skończo- nych, metodę różnic skończonych oraz metodę elemen- tów brzegowych [6,10]. Od wyboru metody obliczeniowej zależy często wiarygodność uzyskiwanych wyników symulacji numerycznych, a dotyczy to przede wszystkim

możliwości uwzględniania w modelu, bądź nie, mnogości zjawisk procesu hartowania oraz ich sprzężeń. Z tego też powodu w pracy zastosowano metodę elementów skoń- czonych.

Dotychczas brak jest kompleksowych modeli numerycznych hartowania grupy węglowych stali narzędzio- wych. Istniejące modele są w jakimś stopniu wycinkowe.

Dotyczą najczęściej przemian fazowych w procesie chłodzenia, nie kojarząc ich z generującymi się w takich procesach naprężeniami. Często też w modelach już opracowanych nie uwzględnia się odkształceń transformacyjnych, które mają miejsce w obróbce cieplnej (szczególnie w hartowaniu).

2. ZJAWISKA MECHANICZNE

W modelowaniu zjawisk mechanicznych równania równowagi przyjęto w formie prędkościowej. Daje to możliwość uwzględniania zmian wielkości termofizycznych od temperatury i składu fazowego w kolejnych przyrostach obciążenia [1,2]

( )

^x^t ⁰ ^σ ^σ^T

σ& & &

o = =

∇ , , ⁽¹⁾

gdzie: σ=σ(σαβ) jest tensorem naprężenia, znak (°) oznacza niepełny iloczyn wewnętrzny, x=x(xα) jest wektorem położenia rozważanej cząstki (punktu).

Równania (1) uzupełniają związki konstytutywne, które przy założeniu addytywności odkształceń mają postać:

(

^ε ^ε ^ε ^ε

)

^E ^ε^e

E

σ&= o &−&^Tph−&^p−&^tp +&o ⁽²⁾ gdzie: E=E(T, Ση^k) jest tensorem sprężystości zależnym od temperatury (T) i składu fazowego (Ση^k), ε^e jest tensorem odkształceń sprężystych, ε^Tph izotropowym tensorem odkształceń cieplnych i odkształceń strukturalnych, ε^p jest tensorem odkształceń plastycznych, natomiast ε^tp tensorem odkształceń transformacyjnych.

Wielkości termofizyczne, takie jak moduł Younga (E) i moduł styczny (Et), uzależnione są od temperatury i składu fazowego liniowo:

( ) ( )

(

^,

) ( )

^, ¹^,...,⁵

, ,

=

∑

k T E T

E

T E T

E

k tk k

t

k k k

η η

η

η (3)

gdzie: ηk jest udziałem k-tej fazy w strukturze materiału.

Równania równowagi (1) uzupełnia się warunkami początkowymi

(

^x

)

^σ

( )

^x ⁰ ^ε

(

^x

)

^ε

( )

^x ⁰

σ ,t0 = 0 = , ^e ,t0 = 0 = (4)

oraz warunkami brzegowymi na częściach brzegów (Γ) gdzie odebrane są stopnie swobody, tzn.:

(

^x

)

^U

(

^U

(

^x

)

⁰

)

U = =

Γ Γ

t

t , ,

, & &

& (5)

(3)

gdzie: U jest wektorem przemieszczenia.

W wyznaczaniu odkształceń plastycznych stosuje się model nieizotermicznego plastycznego płynięcia z warunkiem plastyczności Hubera-Misesa i wzmocnieniem izotropowym lub wzmocnieniem kinematycznym. Zatem funkcje płynięcia plastycznego (f=f( ,Y)) przyjmowane są w postaci:

( )

(

^,

)

⁰

, 0 , ,

0

p

=

−

=

−

=

∑

k ef

ef k ef

T Y f

η σ

ε η

σ ₍₆₎

gdzie

σ

_ef jest naprężeniem efektywnym,

ε

_ef^p _efektyw-

nym odkształceniem plastycznym, Y=Y(T,Ση^k,

ε

_ef^p _{) jest}

naprężeniem uplastyczniającym dla materiału o składzie fazowym Ση^k w temperaturze T i uplastycznieniu

ε

_ef^p _,

Y0=Y0(T,Ση^k) jest granicą plastyczności.

Odkształcenia plastyczne determinuje stowarzyszone prawo płynięcia plastycznego [8]:

0 , 0

p ,

=

= Λ

= & f f& f

&

σ ε

∂

∂ (7)

gdzie:

Λ

jest skalarnym mnożnikiem plastyczności.

Odkształcenia plastyczne wynikające z (7) dla wzmocnienia izotropowego i wzmocnienia kinematycznego wyznaczane są następująco:

Y 2

p 3S

ε& = &Λ ⁱ

0 p

2 3

Y α

ε S−

Λ

= &

& ⁽⁸⁾

gdzie: S jest dewiatorem tensora naprężenia (S= – Ikk/3), α jest tensorem przemieszczania się środka powierzchni plastycznego płynięcia [10].

Do wyznaczania odkształceń transformacyjnych zastosowano zmodyfikowaną formułę Leblonda wprowa- dzając funkcję obniżającą szybkość zmian odkształceń transformacyjnych z upływem trwania przemiany [9], tzn.:

( ) ( )







≥

−

≤

=

∑

⁼₌ ¹ ^ln ^, ⁰^.⁰³

03 . 0 , 0

5

2 1

1

k k k k

k k k

k tp

Y η K η η η

η

&

& S

ε ⁽⁹⁾

gdzie: ^ph

k

K1k=3ε1 są objętościowymi odkształceniami strukturalnymi przy przejściu z fazy wyjściowej „1” na k-tą fazę, Y1 jest granicą plastyczności fazy wyjściowej (miękkiej).

Naprężenie uplastyczniające, w przypadku stosowania modelu materiału ze wzmocnieniem izotropowym, wyznaczane jest zależnością:

( )

0

( ) (

^p

)

p

, , ,

,

, _k _ef Y T _k Y_H T _k _ef

T

Y

∑

η ε =

∑

η +

∑

η ε ⁽¹⁰⁾

Po wprowadzeniu modułów: wzmocnienia ( = (T,Ση^k,

p

ε

ef )), osłabienia termicznego (^T=^T(T,Ση^k,

ε

_ef^p _{)) oraz}

osłabienia bądź wzmocnienia strukturalnego (

κ

^η^r=^η(T,Ση^k,

ε

_ef^p )) prędkość zmian granicy plastycz- ności wyznacza zależność:

( ∑ )

= + +

∑

k k T

p p k

k ef

ef T

T

Y& , η ,ε κε& κ & κ^ηη& (11)

gdzie

( )

( ) ( )

k p ef k H

k k

k

p ef k H

T k

p ef

p ef k H

p ef

T Y T

Y Y

T T Y T

T Y T Y

T Y Y

k

η ε η η

η η

κ

ε η η

κ

ε ε η ε

κ

η

∂

∂ +

∂

=

∂

= ∂

∂

∂ +

∂

=

∂

=∂

∂

=

∂

= ∂

∑

, , ,

, ,

0

0 (12)

W przypadku stosowania modelu materiału ze wzmocnieniem anizotropowym (kinematycznym) pręd- kość zmian granicy plastyczności określana jest wzorem:

( ) ( )

∑

∑ ∑

∑

+

=

∂

∂ +

∂

=

k k T

k

k k

T k

T T Y

T T Y Y

η κ κ

η η

η &

&

& ₀ , ₀ ,

0 (13)

W iteracyjnym procesie wyznaczania odkształceń plastycznych wykorzystuje się algorytm Newtona- Raphsona [5,10].

3. WYNIKI SYMULACJI

Wykorzystując przedstawione modele zjawisk hartowania, zbudowano program komputerowy służący do symulacji zjawisk hartowania. Wykonano symulacje numeryczne hartowania elementu wykonanego z węglo- wej stali narzędziowej C80U. Strukturą wyjściową był perlit (sferoidyt). Symulacje dotyczyły hartowania po niepełnym austenityzowaniu (hartowanie przypo- wierzchniowe).

Jak już wspomniano, w symulacji zjawisk mechanicznych moduł Younga i moduł styczny uzależniono od temperatury, natomiast granicę plastyczności - od temperatury i składu fazowego hartowanego obiektu.

Moduł Younga i moduł styczny (E i Et) były równe odpowiednio E=2×10⁵ i Et=1,1×10⁴ [MPa] (Et=0,055E), natomiast granice plastyczności (Y0k): 150, 450, 1100 i 300 [MPa], odpowiednio dla austenitu, bainitu, mar- tenzytu i perlitu (oraz ferrytu), w temperaturze 300 K.

W temperaturze solidusu (1700 K) i wyższej moduł Younga i moduł styczny przyjęto o wartościach 100 i 5.5 [MPa], natomiast granice plastyczności były równe 5 [MPa]. Wielkości te przyjęto na podstawie danych zamieszczonych w pracach [1,4,7]. Wykorzystując dys- kretne wartości wielkości termofizycznych zależne od temperatury, dokonano ich aproksymacji sklejonymi funkcjami kwadratowymi.

Założone wartości modułu Younga, modułu styczne- go oraz granicy plastyczności dla perlitu (sferoidytu),

(4)

w temperaturze 300 K potwierdzają wyniki uzyskane ze statycznej próby rozciągania próbki z rozważanej stali.

W zależności od przykładu założono, że materiał charak- teryzuje się wzmocnieniem izotropowym, bądź wzmocnieniem kinematycznym.

Symulacji hartowania poddano obiekt osiowosyme- tryczny o wymiarach φ25×50 mm (rys. 1).

Rys. 1. Schemat rozważanego układu oraz przyjęte warunki brzegowe

Nagrzewanie obiektu realizowano warunkiem brzegowym typu Newtona symulującym nagrzewanie w złożu fluidalnym o temperaturze 1550 K. Temperatu- ra początkowa nagrzewanego obiektu była równa 300 K, a strukturą wyjściową był perlit (sferoidyt). Współczyn- nik przejmowania ciepła złoża fluidalnego przyjęto stały (niezależny od temperatury) o wartości 3800 W/(m²K).

Na powierzchniach czołowych obiektu przyjęto mniejszy współczynnik przejmowania ciepła o wartości 2000 W/(m²K). Uwzględniono w ten sposób utrudniony (gorszy w praktyce) opływ warstwy fluidalnej na powierzchniach czołowych nagrzewanego w ten sposób obiektu.

Rozkład temperatury i uzyskana strefa austenitu po zakończeniu nagrzewania przedstawiono na rys.2.

Rys. 2. Rozkłady po zakończeniu nagrzewania; a) temperatury [K], b) udziału fazowego austenitu

Symulację nagrzewania kontynuowano do osiągnięcia maksymalnej temperatury 1350 K w otoczeniu punktu 3 (rys. 1). W symulacji przemiany fazowej nagrzewania (przemiany austenitycznej) wykorzystano wykres CTPA

odczytując z niego Ac1(t), Accm(t).

Chłodzenie nagrzanego obiektu realizowano również warunkiem brzegowym typu Newtona, przyjmując wartość współczynnika przejmowania ciepła równą 3800 W/(m²K) (chłodzenie w złożu fluidalnym). Tempe- ratura medium chłodzącego (T∞) była równa 300 K.

Entalpie przemian fazowych przyjęto o wartościach:

H2=HA-B=314×10⁶, H4=HA-M=630×10⁶, H3=H5=HA-P=800×10⁶, [J/m³]. Symulację chłodzenia kontynuowano do uzyskania w całym obiekcie temperatury równej 300 K.

Rys. 3. Kinetyka przemian w punktach 1 i 2 przekroju A-A Rozkłady naprężeń własnych, po procesie hartowania, przedstawiono na rysunkach 4÷8. Na rysunkach porównano naprężenia uzyskane z symulacji przy zało- żeniu wzmocnienia izotropowego (rysunki po lewej stronie), bądź kinematycznego (rysunki po prawej stronie).

Rys. 4. Rozkłady naprężeń własnych promieniowych

Rys. 5. Rozkłady naprężeń własnych stycznych

0 3 6 9 12

r, mm

0 3 6 9 12

r, mm

0 5 10 15 20 25

z, mm

T

950 1000 1050 1100 1150 1200 1250 1300 1350 1400

a) b)

0 3 6 9 12

r, mm

0 3 6 9 12

r, mm

0 5 10 15 20 25

A

T A

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

0.1 1.0 10.0 100.0

Czas, s

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Ułamek fazy

perlit

austenit martenzyt

bainit pkt 2, r=12.5 mm

Ac₁(t), Ac_cm(t)

prz. A-A, z=12.5 mm perlit

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25 b)

r

σ a)

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25

r

σ c)

-700-600

-500

-400

-300

-200

-100

0

100

200300

400

r

σ

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25 b)

a)

rz

σ

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25

rz

σ c)

-250

-225-200

-175

-150-125

-100-75

-50

-250

25

σrz

(5)

Rys. 6. Rozkłady naprężeń własnych obwodowych

Rys. 7. Rozkłady naprężeń własnych osiowych

Rys. 8. Rozkłady naprężeń własnych efektywnych Na rys .9 przedstawiono rozkłady naprężeń własnych wzdłuż promienia. Ac1(t)i Accm(t) oznaczają temperatury początku i zakończenia austenityzowania z uwzględ- nieniem szybkości nagrzewania (wykres CTPa).

Rys. 9. Rozkłady naprężeń własnych wzdłuż promienia, prze- krój A-A.

W celu ułatwienia analizy generujących się naprężeń w procesie hartowania rozważanego obiektu, na rys. 10 (w wyróżnionych punktach) przedstawiono historię generowania się izotropowych odkształceń termicznych i odkształceń strukturalnych (ε^Tph). Końcowe wartości tych odkształceń (czas 200 s), to odkształcenia struktu- ralne (ε^ph).

Rys. 10. Historia generowania się izotropowych odkształceń termicznych i strukturalnych (^Tph) w punktach 1 i 2 przekroju

A-A (rys.1)

Na kolejnych rysunkach przedstawiono historię generowania się naprężeń hartowniczych w wyróżnionych punktach przekroju A-A (rys.1).

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25 b)

a)

σϕ

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25

ϕ

σ c)

-1200-1100-1000-900-800-700-600-500-400-300-200-1000100200300400

σϕ

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25 b)

a)

z

σ

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25

z

σ c)

-1250

-1050

-850

-650

-450

-250

-50

150

350

550

z

σ

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25 b) a)

σef

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25

σef

c)

0

100200

300400

500600

700800

900

10001100

1200

σef

(6)

Rys. 11. Historia generowania się naprężeń w otoczeniu punk- tów 1 i 2 przekroju A-A. Wpływ rodzaju wzmocnienia Rozkłady efektywnych odkształceń plastycznych przedstawiono na rys.12, natomiast rozkłady tych od- kształceń wzdłuż promienia w przekrojach A-A i B-B - na rys. 14. Przedstawiono również rozkłady efektywnych odkształceń transformacyjnych (rys.13).

Rys. 12. Rozkłady efektywnych odkształceń plastycznych (×10^-3)

Rys. 13. Rozkłady efektywnych odkształceń transformacyjnych (×10^-3)

Rys. 14. Rozkłady efektywnych odkształceń plastycznych wzdłuż promienia, przekroje A-A i B-B (rys.1)

4. PODSUMOWANIE

Rozkłady naprężeń uzyskane z zastosowaniem w modelu zjawisk mechanicznych wzmocnienia izotropowego lub kinematycznego są podobne, a ich wartości można przyjąć za porównywalne (rys. 4÷8). W obydwu przypadkach rozkłady znaczących normalnych naprężeń hartowniczych są korzystne (ściskanie w warstwach przypowierzchniowych). Stosując wzmocnienie kinematyczne, otrzymano nieznacznie niższy poziom naprężeń normalnych i efektywnych (rys. 11). Odzwierciedla to również historia generowania się naprężeń hartowniczych. Uzyskany rozkład odkształceń plastycznych jest również porównywalny (rys. 12), natomiast rozkład odkształceń transformacyjnych jest nieco inny w przypadku wzmocnienia izotropowego, a inny w przypadku wzmocnienia kinematycznego (rys. 13). Różnice pomię- dzy uzyskiwanymi wynikami symulacji zjawisk mechanicznych, w zależności od wyboru rodzaju umacniania

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25

z, mm

b)

a)

p

εef

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25 c)

p

εef

0

2

4

6

8

10

12

14

16

18

20

p

εef

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25

z, mm

b) a)

0 3 6 9 12

r, mm

0 3 6 9 12

0 5 10 15 20 25 tp c)

εef ^tp

εef

0123456789101112131415

tp

εef

(7)

się materiału hartowanego obiektu, to skutek sposobu hartowania. Symulacja zarówno nagrzewania, jak i chłodzenia sprawia, że w procesie hartowania mają miejsce zarówno dociążenia, jak i odciążenia w historii generowania się naprężeń hartowniczych. Również

rozkłady efektywnych odkształceń są korzystne [rys. 12].

Umocnienie materiału ma miejsce w warstwach przypowierzchniowych, a to oznacza, że warstwa wierzchnia po takim hartowaniu ma dobre własności mechaniczne (podwyższona granica plastyczności). .

Literatura

1. Bokota A., Domański T.: Numerical analysis of thermo-mechanical phenomena of hardening process of elements made of carbon steel C80U. “Archives of Metallurgy and Materials” 2007, nr 2, 52, .. 277-288.

2. Bokota A., Domański T.: Model numeryczny posuwowego hartowania elementów maszyn wykonywanych ze stali narzędziowych. W: Informatyka w technologii metali. Materiały XII Konferencji KomPlasTech2005, nr 9. Ustroń 16-19 stycznia 2005, s. 65-72.

3. Bokota A., Iskierka S.: Effect of phase transformation on stress states in surface layer of laser hardened carbon steel. ISIJ International 1996, t. 36, nr 11, p. 1383-1391.

4. Coret M., Combescure A.: A mesomodel for the numerical simulation of the multiphasic behavior of materials under anisothermal loading (application to two low-carbon steels). “International Journal of Mechanical Sci- ences” 2002, p. 1947-1963.

5. Gür H. C., Tekkaya E. A.: Numerical investigation of non-homogeneous plastic deformation in quenching process. “Materials Science and Engineering” 2001, p. 164–169.

6. Liu C., Xu X., Liu Z., A FEM modeling of quenching and tempering and its application in industrial engineering. “Finite Elements in Analysis and Design” 2002, p. 1053-1070.

7. Pietrzyk M.: Through-process modelling of microstructure evolution in hot forming of steels. “Materials Process- ing Technology” 2002, p. 53-62.

8. Raniecki B., Tanaka K.: On the thermodynamic driving force for martensitic phase transformation. ICRS 3, Residual Stresses III, Science and Technolgy, t. 1, 1992, p. 196-201.

9. Taleb L., Sidoroff F.: A micromechanical modelling of the Greenwood-Johnson mechanism in transformation induced plasticity. “International Journal of Plasticity” 2003, p. 1821-1842.

10. Zienkiewicz O.C. and Taylor R.L.: The finite element method. T1,2. Fifth ed. Oxford: Butterworth-Heinemann, 2000.