Zagadnienia odwrotne w termodynamice procesów odlewniczych

P O L I T E C H N I K A Ś L Ą S K A

ZESZYTY NAUKOWE

RZYMKOWSKI

ZAGADNIENIA O DW RO TN E W TERMODYNAMICE

PROCESÓW ODLEW NICZYCH

G L I W I C E

1 9 9 1

S P I S TR EŚCI

1. W S T Ę P ... 9

.2. PRZEGLĄD LITERATURY ... 11

3. PRZYBLIŻONA METODA WYZNACZANIA PARAMETRÓW OPISUJĄCYCH PROCES WYMIANY CIEPŁA NA POWIERZCHNI WLEWKA . ... 14

3.1. Wstęp ... 14

3.2. Założenia ... 15

3.3. Metoda rozwiązania . .'... 17

3.4. Wybór przybliżenia początkowego ... 19

3.5. Model numeryczny ... 21

3.6. Przykłady obliczeń ... 25

3.7. Podsumowanie ... 32

4. WYZNACZANIE WARUNKÓW WYMIANY CIEPŁA NA POWIERZCHNI WLEWKA NA PODSTAWIE POŁOŻENIA FRONTU KRZEPNIĘCIA ... 33

4.1. Wstęp ... 33

4.2. Założenia ... 34

4.3. Metoda rozwiązania ... 36

4.4. Model numeryczny... • •• 38

4.5. Przykłady o b l i c z e ń ... ... ... 51

4.6. Podsumowanie ... *... • •* 56

5. DOBOR TEMPERATURY ZALEWANIA I PRĘDKOŚCI. WYCIĄGANIA ZAPEWNIAJĄCYCH WŁAŚCIWE STYGNIĘCIE WLEWKA CIĄGŁEGO ... 57

5.1. Wstęp ... 57

5.2. Model matematyczny ... -. 58

5.3. Metoda r o z w i ą z a n i a ... 60

5.4. Model numeryczny ... 63

5.5. Przykłady obliczeń • • • • 5.6. Podsumowanie ... 76

PRZEBIEGU PROCESU KRZEPNIĘCIA ... 77

6.1. W s t ą p ... 77

6.2. Założenia ... 77

6.3. Model matematyczny ... 78

6.4. Metoda rozwiązania ..;... 79

6.5. Przykłady obliczeń ... 82

6.6. Podsumowanie ... 89

7. NUMERYCZNA METODA WYZNACZANIA OPORU CIEPLNEGO SZCZELINY NA STYKU ODLEWU I FORMY ... 90

7.1. Wstęp ... 90

7.2. Założenia ... 90

7.3. Model m a t e m a t y c z n y ... 91

7.4. Metoda rozwiązania ... 93

7.5. Przykład obliczeniowy dla obiektu rzeczywistego ... 95

7.6. Podsumowanie .... ... 98

8. WYZNACZANIE ZASTĘPCZEJ POJEMNOŚCI CIEPLNEJ ... 99

8.1. Wstęp ... 99

8.2. Model matematyczny ... 100

8.3. Metoda rozwiązania ... 104

8.4. Przykłady o b l i c z e ń ... 107

6.5. Podsumowanie ... 116

9. LITERATURA ... 117

10. WNIOSKI ... 133

DODATEK: METODY MINIMALIZACJI FUNKCJI ... 136

D.l. Wstęp ... . ... 136

D.2. Klasyfikacja metod m i n i m a l i z a c j i ... 137

D.3. Metody poszukiwania minimum w kierunku ... 140

D.4. Metody kierunków poprawy ... 147

D.4.1. Metoda największego spadku ... 147

D.4.2. Metoda gradientu prostego ... 148

D.4.3. Metody gradientu sprzężonego ... 149

D.4.4. Metody zmiennej metryki ... 150

D.5. Literatura ... 152

S T R E S Z C Z E N I A 153

C ONTENTS

1. INTRODUCTION ... 9

2. R E F E R E N C E S ... .. il 3. APPROXIMATE METHOD OF DETERMINATION OF PARAMETERS DESCRIBING HEAT TRANSFER PROCESSES ON THE SURFACE OF CONTINUOUS CASTING .... 14

3.1. Introduction ... 14

3.2. Assumptions ... 15

3.3. The method of solution ... 17

3.4. The choice of initial approximation ... 19

3.5. Numerical model ... 21

3.6. Examples of computations ... 25

3.7. Summary ... 32

4. DETERMINATION OF HEAT EXCHANGE CONDITION ON THE SURFACE OF CONTINUOUS CASTING ON THE BASIS OF ASSUMED SOLIDIFICATION FRONT .. 33

4.1. Introduction ... 33

4.2. Assumptions ... 34

4.3. The method of solution ... 36

4.4. Numerical model ... 38

4.5. Examples of computations ... 51

4.6. Summary ... 56

5. THE CHOICE OF POURING TEMPERATURE AND PULLING RATE ASSURING THE PROPER CONDITIONS OF CONTINUOUS CASTING C O O L I N G ... 57

5.1. Introduction ... 57

5.2. Mathematical description ... 58

5.3. The method of solution ... 60

5.4. Numerical model .... 63

5.5. Examples of computations .... 73

5.6. Suramar y. ... 76

6. DETERMINATION OF THICKNESS OF MOULD ASSURING THE ASSUMED

SOLIDIFICATION PROCESS ... 77

6.1. Introduction ... 77

6.2. Assumptions ... 77

6.3. Mathematical model ... 78

6.4. The method of solution ... 79

6.5. Examples of computations ... 82

6.6. Summary ... 89

7. NUMERICAL METHOD OF AIR-GAP THERMAL RESISTANCE IDENTIFICATION ON THE CAST-MOULD CONTACT SURFACE ... 90

7.1. Introduction ... 90

7.2. Assumptions ... 90

7.3. Mathematical model ... 91

7.4. The method of solution ... 93

7.5. The example of callculations for real object ... 95

7.6. S u m m a r y ... 98

8. DETERMINATION OF SUBSTITUTE THERMAL CAPACITY ... 99

B.l. I n t r o d u c t i o n ... 99

8.2. Mathematical model ... 100

8.3. The method of solution ... 104

8.4. Examples of computations ... 107

8.5. Summary ... 118

9. REFERENCES ... 117

10. CONCLUSIONS ... 133

APPENDIX: THE METHODS OF FUNCTION MINIMIZATION ... 136

D.l. Introduction ... 136

D.2. Classification of minimization methods ... 137

D.3. The method of minimum searching in the direction ... 140

D.4. The methods directions correction ... 147

D.4.1. The method of the maximum decrease ... 147

D.4.2. The method of simple gradient ... 148

D.4.3. The method of coupled gradient ... 149

D.4.4. The method of variable metric ... 150

D.5. References ... 152

S U M M A R Y ... 153

C O Æ E P X A H E

1. B B E Æ E H M E ... ... 9

2 . OCMOTP JIHTEPATyPU ... 1 1

3 . flPHEUHSEHHblR METOfl OnPEflEAEHHH nAPAMETPOB, OnHCUBAJÛHHX I1POUECC

TEHJIOOEMEHA HA flOBEPXHOCTH C JU 1T K A ... 1 4

3 . 1 . B B e a e H H e ... 1 4

3 . 2 . OcHOBHbie n o a o i e H H n ... 1 5 3 . 3 . H e T o a p e a e H H « ... 1 7 3.4. Bhiöopiea H a s a / i b H o r o npHÖaHaeuHn ... 19 3.5. M H c a e H H a s H o a e a i> ... 21

3 . 6 . r i p u n e p u BkiH H caeH H tt ... 2 5

3 . 7 . P e 3 M n e ... 3 2

4 . OnPEAEJIEHHE yCAOBHR TEH-flOOEMEHA HA nOBEPXHOCTH CAHTKA HA

OCHOBAHHH PACnO^OXEHHH PPAHHUH 3ATBEPÆEBAHMS ... 3 3 4 . 1 . B B e n e H H e ... . . 3 3 4 . 2 . O cH O BH ue noTioseH M B ... 3 4 4 . 3 . MeTOA peneH M B ... 3 6 4 . 4 . H H C ueH H afl MOfle/ik ... 3 8 4 . 5 . npH M epbl BHMHCJieHH« ... 5 1

4 . 6 . P e s B M e ... 5 6

5 . I10ÆE0PKA TEM ÍlEPATyPH PA3JIHBKH H CKOPOCTH EHTJirHBAHH»,

OEECnEHHBAIOfflHX ilPABHJlbHOE OCTUBAHHE H EIlPEPltB H O rO C A H T K A ... 5 7 5 . 1 . ... 5 7 5 . 2 . H aT e«aT O M eC K aH M o a eJib ... 5 8

5 . 3 . MeTOfl p e S B H M B ... 6 0

5 . 4 . HKCJieH H aa w o a e /ib ... 6 3

5 . 5 . rip H M ep u B H M H C /ieH H ft... 7 3

5 . 6 . P e 3 » n e ... 7 6

nPEJUlOAATAEHOrO X O M riPOUECCA 3ATBEPJ1EBAHHÎÎ... 77

6.1. BBeaeHHe... 77

6.2. OcHOBHue nonoxeHHH ... 77

6.3. HaTeiiaTHHecKafl noaenb ... 78

6.4. HeToa pe*ennn ... 79

6.5. npHMepH BhlMHC/ieHHfl ... 82

6.6. P e 3 » n e ... 89

7. 4HCAEHHUR METOÆ OnPEHEJIEHHfl TEByiOBOTO COnPOTHBJIEHHS 3A30PA MEXay OTAHBKOR H ÍOPMOR ... 90

7.1. BseaeHHe ... 90

7.2. OcHOBHue noaoseHMB... ... ... 90

7.3. UaxeHaTHMecKafl Hoaeai> .... 91

7.4. MeToa p e s e H H H ... 93

7.5. SHCaeHHUfl npnuep ann peanbHoro o6beKTa ... 95

7.6. Pe3ioMe... 98

8. OflPEaEAEHME 34>*EKTHBH0r0 K0344HUHEHTA TEnAOEMKOCTH... 99

8.1. BseaeHHe... 99

8.2. MaTeMaTHHecKa« noaenb ... 100

8.3. MeTOa peaeHua ... 104

8.4. npHMepu BbiHHcaeHHR ... 107

8 5. Peoione... 116

9. AHTEPATYPA ... 117

1 0 . ilOABEÄEHHE H T O r O B ... 133

nPHAOIEHME: METOilU MHHMMA7IH3AÜHH *yHKUMH ... 136

D.I. BseaeHHe ... 136

D. 2. IOiaCCH<ł>MKaUMH HeTOaOB MHHHHaaHSaUHH... 137

D. 3. Heroau HaxoiaeHua HHHHMyHa no HanpaBaeHHn... 140

D.4. Meroau KoppeKTHpoBaHHoro HanpaBaeHHn ... 147

D. 4.1. MeToa HaHcKopeftiaero cnycKa ... 147

D. 4.2. MeToa rpaaeHTa ... 148

D.4.3. Me-roau conpaieHHoro rpaaeHTa ... 149

D.4.4. Heroau nepeaeHofl MerpHKH ... 150

D. 5. ÜHTepaTypa ... 152

PE3DME 153

i. WSTĘP

W ostatnich latach w ramach ogólnej teorii przepływu ciepła zostały sformułowane i są intensywnie rozwijane nowe kierunki badań wykorzystujące metody rozwiązywania tzw. zadań odwrotnych transportu ciepła. W teorii przewodnictwa cieplnego wyróżnia się kilka typów zagadnień odwrotnych, które umownie można pogrupować następująco:

- zadania początkowe, polegające na wyznaczeniu warunków początkowych opisujących stan obiektu w chwili przyjętej za początkową,

- zadania brzegowe, polegające na wyznaczeniu warunków brzegowych (np.

współczynników wymiany ciepła, absorpcyjności),

- zadania geometryczne, polegające na określeniu niektórych cech jooraetrycznych obiektu,

- zadania parametryczne (inwersyjne), polegające .na wyznaczeniu parametrów fizycznych materiałów (np. przewodności cieplnej, dyfuzyjności cieplnej, ciepła właściwego, ciepła przemian fizycznych),

- zadania mieszane, polegające na jednoczesnym wyznaczaniu cech różnych typów (np.współczynnika Wymiany ciepła i ciepła przemian fizycznych).

Każde z wymienionych zadań może być jeszcze sformułowane jako zadanie identyfikacyjne lub projektowe. Źródłem informacji dla zagadnień identyfikacyjnych są pomiary, a dla zagadnień projektowych wymogi technologiczne prowadzenia procesu. Możliwości wykorzystywania metod rozwiązywania zadań odwrotnych są jak dotąd mało wykorzystywane dla nieliniowych, wielowymiarowych zadań z ruchomymi brzegami, między innymi w termodynamice procesów odlewniczych.

Celem prezentowanej pracy wykonanej w latach 1981—90 w ramach Centralnego Programu Badań Podstawowych 02.09. (dawniej MR 20) była adaptacja i rozszerzenie pewnych metod rozwiązywania zadań, odwrotnych do obliczeń związanych z krzepnięciem odlewów w formie (w szerokim rozumieniu tego słowa). Obok wyników teoretycznych prezentowane będą również efektywne

algorytmy realizujące obliczenia numeryczne; wykorzystano je do rozwiązywania zadań spotykanych w teorii i praktyce procesów odlewniczych, przeprowadzono również weryfikację niektórych rezultatów. Rozważano prawie wszystkie wymienione wcześniej typy zagadnień odwrotnych, w tym między innymi następujące.problemy;

- identyfikacja warunków wymiany ciepła na powierzchni wlewka ciągłego (pionowe i łukowe wlewki prostokątne, pionowe i łukowe wlewki okrągłe), - identyfikacja oporu cieplnego szczeliny gazowej na styku odlewu i formy, - identyfikacja grubości warstw formy wielowarstwowej,

- wyznaczanie optymalnych prędkości wyciągania i temperatur zalewania dla wlewków produkowanych sposobem ciągłym,

- wyznaczenie zastępczej dyfuzyjności cieplnej dla krzepnącego w przedziale temperatur metalu.

Zagadnienia odwrotne wymagają znajomości termicznych odpowiedzi układu. W zależności od rozważanego problemu były to pola temperatury na określonych powierzchniach wewnątrz ciała, przebiegi temperatury w punktach kontrolnych, wartości temperatur w zbiorze izolowanych punktów, położenia powierzchni i linii izoterraicznych oraz profile stygnięcia.

Dużo miejsca w pracy poświęcono problemom odlewania ciągłego.

Technologia ta jest szczególnie wydajna i efektywna, zapewnia wysoką jakość produkowanych wyrobów i w chwili obecnej zdecydowanie wypiera inne metody wytwarzania wlewków. Okazuje się przy tyra, że właśnie metody rozwiązywania zadań odwrotnych transportu ciepła są wyjątkowo przydatne do analizy procesów cieplnych w procesie odlewania ciągłego, a do tej pory fakt ten pozostawał niezauważony.

2. PRZEGLĄD LITERATURY

Metody rozwiązywania odwrotnych zagadnień nieustalonego przewodzenia ciepła doczekały się bardzo obszernej bibliografii, chociaż problematyka ta stała się przedmiotem zainteresowania stosunkowo niedawno. Pierwsze prace z tej dziedziny ukazały się na przełomie lat pięćdziesiątych i sześćdziesiątych. Można tu wymienić prace N .V.Shuraakova (2.84], T.J.

Mirsepassiego (2.70,2.71], E.M. Sparrowa [2.85], G.J.Stolza (2.86], A.G.Temkina [2.89], czy znaną monografię H.S.Carslawa i J.C.Jaegara [2.17], w której wprawdzie autorzy nie wprowadzają pojęcia zadań odwrotnych, ale gdzie podane są pewne metody rozwiązywania zagadnień bezpośrednich, łatwe do zastosowania przy rozpatrywaniu niektórych zadań odwrotnych. Burzliwy rozwój metod rozwiązywania zagadnień odwrotnych nastąpił wraz z rozwojem metod numerycznych w latach siedemdziesiątych. Liczba publikacji w tyra zakresie jest bardzo wielka, przy czym najważniejszymi wydają się być monografie [2.2,2.3,2.4,2.10,2.55,2.57,2.60,2.64,2.79,2.88] oraz prace [2.1,2.5,2.6,2.7,2.8,2.11,2.12,2.14,2.15,2.18,2.19, 2.20, 2.25, 2.44, 2.45, 2.46,2.47,2.48,2.50,2.51,2.52,2.61,2.62,2.63,2.65,2.66,2.67,2.68,2.69,2.76, 2.80,2.81,2.98,2.101 ].

Obszerne przeglądy literatury można znaleźć w opracowaniach [2.4,2.10,2.24,2.55,2.64], a w szczególności w monografii [2.55], gdzie omówionych zostało kilkaset pozycji. W tej sytuacji omówię tylko te prace, które bezpośrednio wiążą się z tematem niniejszej rozprawy.

Zagadnienia odwrotne przepływu ciepła należą do zadań źle uwarunkowanych. Kryteria dobrego uwarunkowania zostały po raz pierwszy podane przez J.Hadamarda [2.43,2.44]. Według kryteriów Hadamarda zadanie jest dobrze uwarunkowane, jeżeli istnieje rozwiązanie zadania, rozwiązanie to jest jedyne i zależy w sposób ciągły od warunków granicznych. Ostatnia własność nazywana stabilnością rozwiązania względem warunków granicznych »a podstawowe znaczenie dla zastosowań, ponieważ chodzi o to, aby małe błędy w

danych wejściowych nie powodowały zbyt, dużych zmian w rozwiązaniach zagadnień.

Wyznaczenie poprawnego rozwiązania wiąże sio z wykorzystaniem specjalnych metod obliczeń tego typu problemów. Metody takie opracowano w ostatnich latach, a ich powstanie jest główną zasługą uczonych ze szkoły A . N .Tikhonova [ 2 . 7 5 , 2 . 9 0 , 2 . 9 1 , 2 . 9 3 , 2 . 9 4 ] . Do metod stabilnych należą głównie metody wariacyjne, takie jak np. metoda quasi-rozwiązania [2.53,2.54,2.56], odchyłek '[2.56,2.92] i regularyzacji [2.74,2.94]. Na podkreślenie zasługuje również metoda oparta na rachunku wyrównawczym [2.58), którą można stosować nie tylko w celu poprawy stabilności, ale przed© wszystkim w celu poprawy zbieżności wyników.

Chociaż wszystkie wymienione wyżej metody pozwalają budować stabilne algorytmy rozwiązania zadań żle uwarunkowanych, to jednakże ich realizacja jest zazwyczaj bardzo złożona. Dlatego też nie zawsze sie z nich korzysta.

Niejednokrotnie wystarczają metody tzw. naturalnej regularyzacji. Metody te w celu uzyskania stabilnego rozwiązania wykorzystują stabilizujące własności zastosowanych algorytmów obliczeniowych i możliwość sterowania niektórymi ich parametrami.

Dla ekstremalnych sformułowań zagadnień odwrotnych przepływu ciepła, a z takich głównie bgdg korzystał, właściwymi metodami są metody gradientowe z wyborem określonej klasy funkcji, w której poszukuje sie rozwiązania.

Przy doborze klasy funkcji uwzględnia sie wszystkie dodatkowe informacje o ich dynamice i strukturze. Stabilność rozwiązania polepsza sie wraz ze zmniejszeniem liczby parametrów, które charakteryzują wybraną klasę funkcji, a podlegają określeniu. Role parametru regularyzacji spełnia liczba iteracji. Metody te rozwijane przez O.M.Alifanova i jego współpracowników mają bogatą bibliografie» w tym monografie

[2.2,2.3,2.4,2.95].

Liczba prac bezpośrednio dotycząca zagadnień odwrotnych przepływu ciepła w procesach odlewniczych jest niewielka. Tematyką tą zajmuje sie też i niewielu badaczy. Autorowi znane są, nie licząc opracowań własnych [2.26,2.27,2.28,2.29,2.30,2.31,2.32,2.33,2.34,2.35,2.36] i opracowań w s p ó ł a u t o r s k i c h [ 2 . 3 7 , 2 . 3 8 , 2 . 3 9 , 2 . 4 0 , 2 . 4 1 , 2 . 4 2 , 2 . 7 2 , 2 . 7 3 , 2 . 8 7 ] , j e d y n i e

- 13 -

pojedyncze, nieliczne prace innych autorów [2.13,2.16,2.21,2.22,2.23,2.59, 2.78,2.82,2.83,2.99,2.100,2.101,2.102], Nie wyczerpują one tematu, tyra bardziej że dotyczą z reguły bardzo uproszczonych zagadnień.

Niewielki stan zaawansowania badart w zakresie wykorzystania zagadnień odwrotnych w termodynamice procesów odlewniczych, a przy tym liczne możliwości zastosowania tych metod w praktyce projektowej uzasadniają podjecie badart w tej dziedzinie. Z analizy dostępnych prac wynika, że najbardziej odpowiednimi metodami poszukiwania rozwiązań typowych zadań odwrotnych związanych z przepływem ciepła w układzie odlew—forma—otoczenie są metody gradientowe do sformułowań ekstremalnych. Tego typu podejście jest rozpatrywane w dalszej części pracy.

P R O C E S WYMIANY C I E P Ł A NA P O W I E R Z C H N I WLEWKA

3 . 1 . WSTĘP

Projektowanie optymalnej technologii wytwarzania wlewków w procesie ciągłego odlewania jest problemem złożonym i wieloetapowym. Jednym z najistotniejszych czynników wpływających na jakość wlewka ciągłego i mogących służyć za kryterium oceny projektowanej technologii jest pole temperatury w objętości krzepnącego metalu. Technologom i projektantom urządzeń do ciągłego odlewania znane są kryteria dotyczące postulowanego rozkładu temperatury. Można zatem sformułować zagadnienie, w którym poszukuje sie warunków brzegowych, tj.warunków wymiany ciepła na powierzchni wlewka czyniących zadość ograniczeniom nakładanym na pole temperatury. Zazwyczaj poszukuje się współczynnika wymiany ciepła a [W/m^]

lub strumienia ciepła q [W/m2], ponieważ w literaturze, np. [3.1, 3.3, 3.8, 3.10], można znaleźć proste analityczne lub graficzne zależności między tymi wielkościami a parametrami, które umożliwiają sterowanie procesem. Na przykład wg Mullera i Jeschara [3.2] współczynnik wymiany ciepła zależy Jedynie od jednostkowego natężenia wody chłodzącej V [ra3/ra2s] i prędkości wypływu wody z dysz u [m/s]. Zależność określająca współczynnik a ma postać:

a = 1 0u+1000(107 +0.6 8 8 u )V (3.1)

przy czym u < K 1 1 V32> a V €<0.003, 0.009>.

Wyznaczone w poszczególnych strefach chłodzenia współczynniki wymiany ciepła a lub strumienie ciepła q pozwalają konstruktorom urządzeń do ciągłego odlewania zaprojektować systemy chłodzenia o parametrach zapewniających właściwe stygnięcie wlewka.

- 15 -

Rozważa się płaskie i okrągłe wlewki, wytwarzane na urządzeniu pionowym do ciągłego odlewania, odlewane ze stałą prędkością wyciągania w=constans. Wyniki obliczeń można przenosić zresztą na wlewki radialne o odpowiednio dużym promieniu krzywizny, a więc wielkogabarytowe radialne wlewki ciągłe. Przyjmuje się, że produkowane wlewki mają grubość lub średnice rdwną 2R i wytwarzane są z metalu krzepnącego w przedziale temperatur <T„,TL> {^-temperatura solidusu, ^- t e m p e r a t u r a likwidusu),

* M

który wlewany do krystalizatora ma temperaturę T , T >TL>Ts . 3.2. Z A Ł O H E N I A

Rys.3.i. Modelowane obiekty; a) wlewek okrągły, b) wlewek plaski P i g .3.1. Modelled objects; a) cylindrical ingot, b) slab

Obszar wlewka traktuje sio Jako obiekt dwuwymiarowy, dla którego pseudoustalone pole temperatury z uwagi na znikomą ilość ciepła przewodzonego w kierunku ruchu wlewka [3.11] można opisać jednowymiarowym, parabolicznym równaniem różniczkowym cząstkowym. W równaniu tym jedna ze współrzędnych przestrzennych (zgodna z kierunkiem wyciągania) spełnia rolę czasu. Orientacja w przestrzeni modelowanego obiektu pokazana jest na r y s . 3.1.

Przyjęto wyżej założenia oraz warunek, Ze znane sa przybliżone przebiegi temperatury f =f (z) w przemieszczających się wraz z wlewkiem

r» n

punktach kontrolnych r=rn* rn€<0,R>, n=l,N , pozwalają sformułować matematyczny opis zagadnienia w postaci:

d T=*r_ma (rm<? T) , 0<r<R , 0<z<Z (3.2)

z r r

T=T* , z=0 , 0<r<R (3.3)

d T=0 , r=0 , 0<z<Z (3.4) r

T (r ,z)£f (z) f z «£0 (z )c(0, Z> , n=l7N (3.5)

n n n

gdzie r i z oznaczają współrzędne przestrzenne, T=T(r,z) jest temperaturą, indeks m określa geometrię wlewka (ra=0 wlewek płaski, ra=l wlewek okrągły), Z oznacza całkowitą długość drogi chłodzenia podzielonej na K stref odpowiadających odcinkom <zki,zfc>, k=l,K, przy czym z0=o »2k=^» a Dn ^z ^ oznacza dziedzinę funkcji f^(z). Natomiast występująca w równaniu (3.2) funkcja x=x(T) oznacza:

» = M v r c * r ‘ (3.6)

gdzie X Jest współczynnikiem przewodzenia c i e p ł a,y gęstością masy, a tt *

c =c (T) zastępczą pojemnością cieplną [3.7] uwzględniającą ciepło przemiany fazowej.

Podane przebiegi temperatury, jak już wspomniano wcześniej, wynikać mogą z wymogów technologicznych, pomiarów lub innych uwarunkowań nakładanych na modelowany obiekt. Funkcje f^ zadawane są zazwyczaj w postaci tabeli wartości (dziedzina i przociwdziedzina są zbiorami liczb), w postaci wykresu lub wzorem w postaci jawnej.

Z założeń i opisu matematycznego wynika, że nieznane są warunki wymiany ciepła na powierzchni wlewka r=R. Okazuje się jednak, że dodatkowe informacje pozwalają wyznaczyć pole temperatury w całym przekroju wlewka oraz odtworzyć brakujący warunek brzegowy.

- 17 -

Poszukiwać będziemy takich współczynników wymiany ciepła c i=c ł(z) lub strumieni ciepła q=q (z ) , aby uzyskać możliwie bliskie postulowanym przebiegi temperatury w punktach kontrolnych. Parametry te powiązane są z polem temperatury wlewka warunkiem brzegowym o postaci:

przy czym funkcje p=p(z,to) i f?=^(T,u>) w zależności od dwuwar Lościowe j funkcji o>=u>(z) określającej rodzaj warunku brzegowego (u>=0 warunek drugiego rodzaju, c*>=l warunek trzeciego rodzaju) zdefiniowane są następująco:

gdzie T jest temperaturą otoczenia, co

3.3. METODA ROZWIĄZANIA

Konstrukcje systemów chłodzenia w urządzeniach do ciągłego odlewania pozwalają założyć, że w poszczególnych strefach chłodzenia współczynniki wymiany ciepła lub strumienie ciepła są stałe i wynoszą odpowiednio lub q fc, k = 1,K . Można więc przyjąć, że funkcja p=p(z,o>) ma postać:

-Xć^T=pf9 , r=R , 0<z<Z (3.7)

(3.8)

natomiast

(3.9)

K

(3.10)

gdzie w zależności od funkcji t*> zgodnie z (3.8)

(3.11)

natomiast

(•)S0

Mi - )=-( (3.12)

)>0 j o .

l 1 . <■>:

Przy powyższym założeniu postawione zadanie można rozwiązać kilkoma sposobami. Zdecydowano sio wykorzystać metodo optymalizacyjną. Do opisu matematycznego dołączono funkcjonał jakości:

N

F ( p ) = ^ T c nJ ( T n- f n)2d z (3.13)

b ę d ą c y K r y t e r i u m w y b o r u n i e z n a n y c h p a r a m e t r ó w p ^ , k = l , K , a w ł a ś c i w i e w

zależności od funkcji w parametrów lub q^, który minimalizowano.

Wyżej T =T(r ,z;p), n=l,N oznaczają odpowiednie rozwiązania n n

zagadnienia brzegowego opisanego związkami (3.2)+(3.4) i (3.7) przy zadanej postaci funkcji p r a stałe C oznaczają parametry wagowo determinujące

n

w i a r y g o d n o ś ć i n f o r m a c j i o p r z e b i e g a c h t e m p e r a t u r y f w p u n k t a c h

n

k o n t r o l n y c h r .

n

Minimum funkcjonału (3.13) poszukiwać bodziemy metodami gradientowymi (patrz np.[3.2]). W tyra celu należy wyznaczyć pochodne funkcjonału względem poszukiwanych parametrów:

V,F=2'\' C | ( T - f ) H dz , k = l , K (3.14)

1 / n I n n k, n

ft— 1 D

gdzie T F- d F(p), natomiast H ^=a T(r ,z;p).

k P lf n p n

k k

>o wyznaczenia funkcji H =<? T wykorzystuje sio wyjściowe zagadnienie

k pk

brzegowe. Różniczkując ( 3 . 2 ) « M 3 . 4 ) i ( 3 . 7 ) po p^ otrzymuje sio K prostych zagadnień brzegowych o postaci:

a H = « r"“ d ( r * d H ) + « ’r""1 [a ( r " a i ) ] B , 0 < r < R , 0 < z < Z

x k r jr k r r k (3.15)

- 19 -

H =0, z=0, 0<r<R

k ^{0 . 1 6 )}

aH=fl, r=0, 0<z<2 r k (3.17)

-\a^Hk=p p^piplz-z^ i)-/j(z-zfc) ] , r=R , 0<z<2 (3.18)

gdzie « ’=d»/dT natomiast <p.=d p, co oznacza, że w zalalności od funkcji u k pk

zgodnie z (3.8) otrzymamy:

f Hk , w=l

\ 0 , co=0

(3.19)

Poszukując minimum funkcjonału, zakłada sie początkowo przybliżenie dla składowych wektora parametrów, t z n . :

a kolejne przybliżenie znajduje sio zgodnie z ogOlnym algorytmem przyjętej metody minimalizacji (patrz Dodatek).

3.4. WYBÓR PRZYBLIŻENIA. POCZĄTKOWEGO

Wybór przybliżenia początkowego dla składowych wektora parametrów

też w praktycznej realizacji procesu obliczeniowego wybór przybliżenia początkowego nie był przypadkowy, a polegał na tym, że poszukiwano go i to oddzielnie dla każdej ze składowych, poczynając od k=l kolejno aż do k=K.

Wykorzystano w tym celu te samą opisaną wyżej metode optymalizacyjną, ale zawążoną tylko do obszaru odpowiadającego jednej strefie chłodzenia. W tyra przypadku układ równań odpowiadający układowi (3.2)+(3.5) odniesiony do strefy o numerze k bgdzie miał postaó:

(3.20)

Pk°\ k=l,K ma istotny wpływ na liczbę iteracji i czas obliczeń. Dlatego

O T=*r m {9 (r"j Tl ] , 0<r<R , z <z<z

z r r fc-i * ^(3.21)

T=Tkl , z=z , 0<r<R (3.22) k-1

a

r k-1 k

T(r ,z)=f (z) , z«^>k(z)c(z ,z > , n=l,N (3.24)

n n n k-1 k K

gdzie Tk ,=T (r , z^ ^) , przy czym T°=T*, a pozostałe Tk *, k>l wynikają z uzyskiwanych kolejno rozwiązań. Natomiast Dk (z) i N^ oznaczają dziedzinę funkcji f (z) odniesioną do przedziału (zt .zt> i liczbę punktów

n k-1 L

kontrolnych w tyra przedziale odpowiednio.

To samo dotyczy i warunku (3.7); mamy teraz:

-XarT=pkp , r =R , zk j<z<zk (3.25)

Jednocześnie w miejsce funkcjonału (3.13) będziemy rozważać funkcjonał:

N

(3.26)

Niewielkiej modyfikacji ulegną równie! formuły na pochodne funkcjonału:

K

7 F = 2 ) C |(T-f )H dz (3.27)

k / nj ^ » k, n

oraz układy równan opisujące występujące w (3.27) funkcje H^. Mają one postać:

3 H =*r~"d lr"SH )-*-*,r'"[d (r”a T) ]H , 0<r<R . z <z< z (3.28)

2 k r r l r r k k-1 k

H k=0 , 2 = z k_. . 0<r<R (3.29)

a fH k= 0 , r = 0 , z M < z i z fc (3.30)

- 21 -

- X 0 r Hk= p kV * > » r = R » Z k - l< Z “ 2k (3.31)

Całe dalsze postępowanie jest analogiczne jak w przypadku globalnym, a wyznaczone dla kolejnych stref p^, k=l,K minimalizujące odpowiedni d l a danej strefy funkcjonał (3.26) brano jako początkowe przybliżenie d l a składowych wektora parametrów, kładąc p^°*=p,_, k=l,K.

3.5. MODEL NUMERYCZNY

Do rozwiązywania występujących w procesie iteracyjnym zagadnień brzegowych można wykorzystać jedną ze znanych metod numerycznych rozwiązywania nieliniowych równań przewodnictwa (3.6). W pracy wykorzystano metodę różnicową. Istotą tej metody jest aproksymacja równania różniczkowego i warunków brzegowych odpowiadającymi im równaniami, różnicowymi, co sprowadza problem do układu równań algebraicznych (tzw.

układ rozwiązujący).

Na modelowany obszar nakłada się równomierną siatkę różnicową f o stałych krokach Ar=R/I i Az=Z/J, t z n . :

Równania różnicowe wyprowadzono opierając się na niejawnej gwieździe

to(z)=0, ze(0,Z),tzn. gdy poszukiwać będziemy strumienia ciepła q, z układu (3.2)+(3.4) i (3.7) otrzymuje się J układów nieliniowych równań algebraicznych kolejno dla j=l,J, o postaci:

czteropunktowej zgodnie z zasadami podanymi w pracy (3.4). Przykładowo, dla

i=0,1 (3.32)

g d z i e :

(« >a

* Az _ . l.J _ 2i-l .m

><*>=____

Ar

‘I—

)

1 = 1 , I

2 (ra+1 ) « J s Jaz b->=i

0,J ArA 2

2*| ‘ ’Az

bH * +

2 * ' s JA2

'•J A r 2 [ 21

J

2 (ra+1 ) ** s *Az c =-(s)

i,J . 2

Ar

< = >A

c (=).—— . 1 = 1 , 1 - 1

‘■J A r 2

<s>

C.,j=°

d (s)=T , 1=0,1-1 l.J l.J-1

[ 4I j '■

2Az

■J-! *(s) Pj

S. j T

który rozwiązuje się iteracyjnio za pomocą algorytmu Thomasa [3.9]. W ostatnich równaniach s jest numerem iteracji, natomiast =T(iAr,jAz),

(T^ j) , podobnie C*(‘)=C*(Tł t) , przy czym T. (= (T'0*+Tjtł^1/>/2, a pj=p(jAz). Rozpoczynając proces iteracyjny jako początkowe przybliżenie przyjmuje się, że 1=0,1. Proces iteracyjny kończymy, jeżeli:

- 23 -

gdzie c jest zadaną dostatecznie małą liczbą, a jako wynik końcowy można przyjąć wynik uzyskany w ostatniej iteracji.

Podobne układy równań uzyskamy i dla zagadnień brzegowych (3.15)+(3.18). Przykładowo przy o>(z)=0 otrzymamy kolejno dla k = i #K J układów równań o postaci:

a Hk +b Hk +C Hk =d , i=0,I (3.34) l.J 1-l.J i.J U i.J l+l.J I.J

gdzie j=l,J, natomiast:

ao / °

* Az _ *. . _ *

i ,j r 2i-i i _____

i = - + , i=i,i-i

,J Ar2 [ 21 J

2*. , Az

2(ra+l)*o jAz 2 (m+1)*’ ^Az

{ w )

b0 j 1 +

2x Az

b . = 1 + --- - ^ 1 1 ---1 T . - 2 T + 1 ---1 T l , 1 = 1 , 1 - 1

Ar “ Ar

bjJ= ¿ P ~ ~

m '

2(m+l)*0 jAz

°-J A r 2

°l.J=

*1 , j*2 f 2i+l "j"

Ar2 [ 2 i

J

C..J=°

d k = l.J

41-1 41

2Az

c ; . j '

W ostatnich równaniach H(kj=Ht (iAr , j Az ) , » ’ (T^), C* =C*(T ), a r? =p(jAz-z )-/j(jAz-z ). Pozostałe parametry mają takie

i.J i.J J «*i k

samo znaczenie jak w układach równań (3.32).

Dyskretyzacje funkcjonału (3.13) przeprowadza sie zazwyczaj na tej samej siatce różnicowej i zadanie to jest proste, jeżeli węzły siatki pokrywają sie z położeniem punktów kontrolnych, tzn. gdy dla każdego n = l ,N istnieją ie { i : i=0,I} takie, że i Ar=r i jednocześnie istnieją takie

n n n

j €:{ j • j =0, J} , źe j AzeD (z); wówczas można przyjąć, żo:

) ’

n * 1 J

v - 2 t e i*— 4 *— • c J [ f . . j - T. .n n n-/ j K .n nj * k = l t K ( 3 ' 36)

n ■ 1 J

g d z i e f . =f (j A z ) , T , = T ( i Ar.j Az) a _{n j n n ł J n n} H k , =Ht ( i A r . j Az ) , _{i ,J k n n}

n n n n n

- 25 -

Zadani© komplikuj© sio, gdy położeni© punktów kontrolnych nie pokrywa sio z wozłarai przyjotej siatki równomiernej. Wówczas do wyboru są dwie drogi. Pierwsza polega na zmianie siatki, t j . przejściu z siatki równomiernej na siatko o zmiennych krokach. Natomiast w drugim przypadku wystopujące w związkach (3.35) i (3.36) i H* wyznacza sio drogą

n ’ n n* n

interpolacji na podstawie wartości funkcji w węzłach sąsiadujących z danym węzłem kontrolnym.

Podobne układy rozwiązujące otrzymuje sio również dla lokalnych zagadnień (3.21) + (3.23) , (3.25) i drugiego zagadnienia globalnego, gdy to(z)=l, tzn. gdy poszukiwać bodziemy współczynnika wymiany ciepła a.

3.6. PRZYKŁADY OBLICZEŃ

Ilustracją funkcjonowania metody mogą być rozwiązania, jakie uzyskano dla wlewka płaskiego i wlewka okrągłego. W obu przypadkach obliczenia przeprowadzono w ten sposób, że najpierw rozwiązano proste zagadnienie brzegowe przy z góry ustalonych parametrach (parametry termofizyczne, temperatura zalewania, prodkość wyciągania, współczynniki wymiany ciepła - dla wlewka płaskiego, strumienie ciepła — dla wlewka okrągłego), a uzyskane tą drogą pole temperatury wykorzystano przy konstruowaniu postulowanych przebiegów temperatury w punktach kontrolnych dla zagadnienia odwrotnego.

Przykład 3.1

Przyjęto, że wlewek o środnicy 2R=0.16(m] odlewany z prędkością w=l[ra/min] wykonany jest ze stali miękkiej. Parametry termofizyczne dobrano na podstawie danych zawartych w [3.5]. Rozważaną drogę chłodzenia wlewka 2 podzielono na cztery sektory o długościach: I—1 (krystalizator), II— 1.5, III-2.5, 1V-10 [m] oraz przyjęto przy tym, że strumienie ciepła w poszczególnych sektorach wynoszą odpowiednio: 1—360000, 11—330000, III—

-225000, IV—180000 [W/m2 ]. Ponadto założono, że metal wlewany do krystalizatora ma temperaturo T*=1540 [°C].

( X ) p o ło żen ie punktów k o n tro ln y c h

Rys.3.2. Siatka różnicowa. Położenie punktów kontrolnych Fig.3.2. Differential mesh. The position of control points

Na modelowany obiekt nałożono siatkę różnicową zawierającą 11 węzłów na promieniu wlewka oraz 240 węzłów w kierunku przesuwu wlewka, (rys.3.2).

Na rysunku tym pokazano również położenie punktów kontrolnych. Rozwiązanie zagadnienia prostego przedstawiono na rys.3.3, założone przebiegi temperatury w punktach kontrolnych na rys.3.4, a rozwiązania zagadnienia odwrotnego na rysunkach 3.5 i 3.6.

- 27 -

Rys.3.3. Rozwiązanie zagadnienia prostego Fig.3.3. Solution of direct problem

R y s . 3 . 4 . P o s t u l o w a n e p r z e b i e g i t e m p e r a t u r y w p u n k t a c h k o n t r o l n y c h

F i g . 3 . 4 . R e q u i r e d t e m p e r a t u r e c o u r s e s a t c o n t r o l p o i n t s

q., = 360 [kW/mZ]

q2=330

q CkW/mZJ q3 = 240

360

Ocou

2

330

© - s e k to ry c h ło d z e n ia , is 1.4

. ——• - rzeczyw isiy s i rumień ciep<a

300.

2 7 0 .

— — — odtworzony strum ień ciep ło

= 361 619 Cw/m2J

. q2 = 323 566

740

210.

' q3 = 252 369 q4 = 172 703 1R0

-

© © ©

—i---1 © odległość DtG

Rys.3.5. Rozwiązanie zagadnienia odwrotnego - strumienie ciepła Fig.3.5. Solution of inverse problem - heat fluxes

Rys.3.6. Rozwiązanie zagadnienia odwrotnego - przebiegi temperatur Fig.3.6. Solution of inverse problem - temperaturo courses

- 2 9 -

Przykład 3.2

Przyjęto, Ze wlewek o grubości 0.21 [m] odlewany z prędkością w=l.l [m/min] wykonany jest ze stali bardzo miękkiej. Parametry termofizyczne materiału dobrano na podstawie danych zawartych w [3.5] oraz załoiono. Ze warunki brzegowe odpowiadają technologii wytwarzania wlewków ciągłych na urządzeniu łukowym USINOK w Dunkierce (pomijając znikomy wpływ promienia krzywizny urządzenia łukowego).

Rozważaną drogę chłodzenia wlewka Z podzielono na siedem sektorów o długościach: 1-0.7 (krystalizator), II-2.6, .111-1.5, IV-3.5, V-3.95, VI—4.25, VII-4 [ra] oraz przyjęto zgodnie z danymi zawartymi w [3.10], Ze współczynniki wymiany ciepła w poszczególnych sektorach wynoszą odpowiednio: 1-1200, 11-950, 1X1-550, IV-430, V-250, VI-225 i VII-200

2

[W/m K]. Ponadto założono, Ze metal wlewany do fcrystalizatora ma temperaturę T*=1551 [°C],

0.105 CnJ | e

Rys.3.7. Siatka różnicowa. Położenie punktów kontrolnych Fig.3.7. Differential mesh. The position of control points

o d l e g ł o ś ć C m j

■ i ■ "I i i i *"

10 12 14 16 18 z

Rys.3.9. Postulowane przebiegi temperatury w punktach kontrolnych Fig.3.9. Reguired teporature courses at control points

1400_

temperaturo CK3

D , = <0.1 7>

D2 * < 2 . 8 > v < 1 2 .17>

D 3 - <5.17>

ł 3

r W V U

1200_

Na modelowany obiekt nałożono siatkę różnicową zawierającą 10 węzłów na połowie grubości wlewka oraz 500 węzłów w kierunku przesuwu wlewka

(rys.3.7). Na rysunku tym pokazano również położenie punktów kontrolnych.

1200.

o d leg ło ść Cnij

Rys.3.8. Rozwiązanie zagadnienia prostego Fig.3.8. Solution of direct problem temperaturo CK3

- 31 -

Rozwiązanie zagadnienia prostego przedstawiono na rys.3.8, założono przebiegi temperatury w punktach kontrolnych na rys.3.9, a uzyskane rozwiązania zagadnienia odwrotnego na rysunkach 3.10 i 3.11.

Rys.3.10. Rozwiązanie zagadnienia odwrotnego - współczynniki wymiany ciepła Fig.3.10. Solution of inverse problem - heat exchange coefficients

odległość Cm]

-L 1---- f— —~~i 1 i---- r----1---- 1 ! .I

0 2 4 6 8 10 12 14 16 18 z

18 z odległość Cm}

1200-

W s p ó łc z y n n ik wym iany

0 -

rzeczyw isty współczynnik a odtworzony w spółczynnik <5, aŁi = 1 2 0 0 a , = 1141

a .2 ' 9 5 0 a2= 9 6 2

tX3 “ 5 5 0 043« 537

044= * 3 0 0 4 = « 2 c t5 = 2 5 0 045= 241 a 6 = 2 2 5 a 6= 2 2 9

o i7 = 047= 207

te m p e ra tu ra CK3 przebieg rzeczyw isty

przebieg odtworzony

Rys.3.11. Rozwiązanie zagadnienia odwrotnego — przebiegi Łemperatuf Fig.3.11. Solution of inverse problem - temperature courses

2 przeprowadzonych obliczeń wynika, źe mimo znacznych zakłóceń przy wprowadzaniu danych wejściowych (w zamieszczonych przykładach przebiegi temperatury zaburzano losowo w punktach kontrolnych do ±50 [°C]) odtworzone wartości strumienia ciepła (przykład 3.1) i współczynnika wymiany ciepła (przykład 3.2) różnią, się nieznacznie od założonych (patrz rys.3.5 i 3.10).

Błąd nie przekracza 8% w przykładzie 3.1 i 5% w przykładzie 3.2. Wygładzone zostały również pola temperatury, a różnica temperatur w punktach kontrolnych była znacznie mniejsza od umyślnie wprowadzonych zakłóceń

(rys.3.6 i 3.11).

3.7. PODSUMOWAŃIE

W rozdziale tyra przedstawiono dogodny do realizacji na EMC algorytm rekonstrukcji warunków brzegowych dla krzepnącego wlewka odlewanego sposobem ciągłym. W modelu wykorzystuje się informację o przybliżonych przebiegach temperatury w ustalonych przemieszczających się wraz z wlewkiem punktach kontrolnych. Dowolność w wyborze położenia punktów kontrolnych i dziedzin funkcji f^ stwarza dodatkową możliwość wykorzystania algorytmu do rekonstrukcji warunków brzegowych również i na podstawie znajomości przybliżonego położenia frontu krzepnięcia.

Prezentowany algorytm oparto na rozwiązaniu odwrotnego nieliniowego zagadnienia brzegowego dla równania przewodnictwa uzyskanego wraźliwościową metodą optymalizacyjną.

Rozwiązanie ww. zagadnienia pozwala dla z góry ustalonych przybliżonych przebiegów temperatury w punktach kontrolnych lub przybliżonego położenia frontu krzepnięcia wyznaczyć warunki chłodzenia na powierzchni wlewka. Wyznaczone na podstawie modelu realne strumienie ciepła lub współczynniki wymiany ciepła mogą być wykorzystane przy projektowaniu systemu chłodzenia urządzeń do ciągłego odlewania.

4. WYZNACZANIE WARUNKÓW WYMIANY CIEPŁA NA POWIERZCHNI NA PODSTAWIE POŁOŻENIA FRONTU KRZEPNIĘCIA

4.1. WSTĘP

W poprzednim rozdziale pokazano, że możliwe jest przybliżone odtworzenie warunków wymiany ciepła na powierzchni wlewka ciągłego na podstawie znajomości przebiegów temperatury w punktach kontrolnych położonych w obszarze wlewka. Podobnym a istotnym dla analizy procesu ciągłego odlewania zagadnieniem jest odwrotne zagadnienie krzepnięcia (tzn.

odwrotne zagadnienie Stefana) polegające na wyznaczeniu warunków wymiany ciepła na powierzchni wlewka na podstawie znajomości (np. z założeń projektowych) położenia frontu krzepnięcia (granicy rozdziału faz) w zadanych punktach kontrolnych, położonych na drodze chłodzenia wlewka.

Próby rozwiązań tego typu zagadnień przy licznych założeniach uprasz­

czających, jak np. pominięcie entalpii przegrzania ciekłego metalu pokazano w pracach [4.1,4.2,4.4,4.5,4.6,4.11,4.14,4.15] Do rozwiązania powyższego zagadnienia można przy dodatkowych założeniach również wykorzystać metodę zaproponowaną w rozdziale poprzednim. Należy założyć, w jakiej temperaturze T z przedziału <T ,T > (T - temperatura solidusu, T - temperatura

s L S L

likwidusu) następuje" rozgraniczenie faz; obrać taki zbiór punktów kontrolnych, by pokrywały się one z położeniem frontu krzepnięcia oraz przyjąć, że w punktach tych temperatura musi przyjmować wartość równą T.

Przybliżenie to będzie jednak zbyt grube, a tyra samym i uzyskane wyniki mogą znacznie odbiegać od rzeczywistych.

Niżej zostanie przedstawiona odmienna od spotykanych w literaturze metoda rozwiązywania tego typu zagadnień, przy czym podobnie jak w rozdziale poprzednim poszukuje się współczynników wymiany ciepła a lub strumieni ciepła q.

Rozważa się odlewane ze stałą prędkością wyciągania w=constans płaskie i okrągłe wlewki wytwarzane na urządzeniu pionowym do ciągłego odlewania, w którym całkowita długość drogi chłodzenia wynosi Z. Przyjmuje się, że produkowane wlewki mają grubość lub średnicę równą 2R i wytwarzane są z metalu krzepnącego w stałej temperaturze T (T - temperatura krystalizacji),

* * — który wlewany do krystalizatora ma temperaturę T , T >T.

Modelowane wlewki traktuje się jako obiekty dwuwymiarowe, dla których pseudoustalone pole temperatury z uwagi na znikomą przewodność ciepła w kierunku wyciągania wlewka [4.8] można opisać jednowymiarowym parabolicznym równaniem różniczkowym cząstkowym. W równaniu tym jedna ze współrzędnych przestrzennych (zgodna z kierunkiem wyciągania) spełnia rolę czasu.

R y s .4.1. Modelowany obiekt Fig.4.1. Modelled object

O r i e n t a c j a w p r z e s t r z e n i m o d e l o w a n e g o o b i e k t u j a k n a r y s . 4. 1 , p o w y ż s z e

z a ł o ż e n i a o r a z w a r u n e k , ż e z n a n e s ą p r z y b l i ż o n e p o ł o ż e n i a f r o n t u

k r z e p n i ę c i a £ w p u n k t a c h k o n t r o l n y c h £ e < 0 , z * > c < 0 , 2 > , n = l ,N p o z w a l a j ą

n n

s f o r m u ł o w a ć m a t e m a t y c z n y o p i s z a g a d n i e n i a w p o s t a c i :

- 35 -

v9 T =a r “'^(r1" # ^ ) , 0<r<a , 0<z<z* (4.1)

"mar(rmd T2) , o<r<R , 0<z<z* (4.2)

v °zr2=!Lir " dr(rm3 T2) , 0<r<R , z*<z<2 (4.3)

Tt=T*. z=0, 0<r<R (4.4)

® T = 0 , r=0 , 0<z<z* (4.5)

a T = 0 , r=0 , z*<z<Z (4.6) r 2

T#=T , v -a , 0<z<z* , e=l ,2 (4.7)

->-iarT1+X2arTJ=wr2Q a ’ , r=<7 , 0<z<z* (4.8)

?nS a ( C n), n=l7Ń (4.9)

gdzie T^=T^(r,n), a^,X^,y^ oznaczają temperatury, dyfuzyjność cieplną, współczynik przewodzenia ciepła i gęstość masy odpowiednio fazy ciekłej e=l oraz fazy stałej e=2, a r i z są współrzędnymi przestrzennymi. Indeks m określa geometrie wlewka (m=0 wlewek płaski, ra=l wlewek okrągły), Q oznacza ciepło krystalizacji, funkcja o=o(z) opisuje położenie granicy rozdziału faz, a a* jest jej pochodną (a*(z)- d a(z )/dz), z* oznacza maksymalną głębokość zalegania ciekłego metalu, a Z ćałkowitą długość drogi chłodzenia podzielonej na K stref odpowiadających odcinkom ć z ^ . z ^ , k=i,K, przy czym z =0 a z,=Z.

O k

Z założeń i opisu matematycznego wynika, że nieznane są warunki wymiany ciepła na powierzchni wlewka r=R, determinujące proces formowania sie naskórka. Poszukiwać wiec będziemy takich współczynników wymiany ciepła a=oi(z) lub strumieni ciepła q=q(z), aby uzyskać możliwie bliskie postulowanym położenia frontu krzepnięcia w punktach kontrolnych. Parametry te powiązane są z polem temperatury wlewka warunkiem brzegowym o postaci:

-X2tfT2=p-*> t r =R , 0<z<Z (4.10)

przy czym funkcje p=p(z,o>) i p=<p(T,co) w zależności od dwuwar tościowe j funkcji u>=co(z) określającej rodzaj warunku brzegowego (t*>=0 warunek drugiego rodzaju, u>=l warunek trzeciego rodzaju) zdefiniowane są następująco:

i

natomiast:

T-T , (*>=1 OD 1 , (*>=o

(4.12)

gdzie jest temperaturą otoczenia.

4.3. METODA ROZWIĄZANIA

Konstrukcje systemów chłodzenia w urządzeniach do ciągłego odlewania pozwalają założyć, że w poszczególnych strefach chłodzenia współczynniki wymiany ciepła lub strumienie ciepła są stałe i wynoszą odpowiednio lub q^, k=l,K. Można więc przyjąć, źe funkcja p=p(z,<*>) ma postać:

P =X ! [ łJlz-2k-i)~p(z_Zk)] Pk (4.13)

gdzie zgodnie z (4.11) składowe wektora parametrów {p} ={ p^, p^, . . . , p^}

oznaczaja:

, w=l

(4.14)

!O=0

natomiast:

0 , ( ■ ) < 0

p ( • )= -^ . (4.15)

1 . ( • ) > O

- 37 -

Przy powyższych założeniach postawione zadanie rozwiązywano wykorzystując metody optymalizacyjne. Do opisu matematycznego dołączono funkcję jakości:

.f

F ( P > - V c (o -ę )2 (4.16) / . n " "

n = 1

kryterium wyboru nieznanych parametrów p^, k=l,K, a właściwie w zależności od funkcji u> parametrów lub q^, którą minimalizowano.

Wyżej o =o(C ;p) oznacza lokalne położenie granicy rozdziału faz w

n n

punktach kontrolnych wyznaczone w wyniku rozwiązania zagadnienia brzegowego (4.1)+(4.8) z warunkiem (4.10) przy zadanej postaci funkcji p, a stałe C oznaczają parametry wagowe determinujące wiarygodność informacji o położeniu frontu £ w punkcie kontrolnym C •

n n

Minimum funkcji (4.13) poszukuje się metodami gradientowymi [4.3,4.16,4.19], W tyra celu należy wyznaczyć pochodne danej funkcji względem poszukiwanych parametrów. Mamy więc:

7 F = 2 V C (o

k n n

n = 1

k=l ,K (4.17)

gdzie V F = d F(p) a 6 = d cr(f ;p). W wyrażeniach (4.17) występuję po N k p k kji p kŁ n

wartości & , n=l,N odpowiadających punktom kontrolnym funkcji kji

O - d o(z;p). Do ich wyznaczenia wykorzystuje się układ równań (4.1)+(4.8) k pk

wraz z warunkiem (4.10). Różniczkując równania (4.1)+(4.8) i (4.10) po p^

otrzymuje się K prostych zagadnień brzegowych w postaci:

w d U = a r ' " d ( r md U ) , 0<r<o , 0<z<z (4.18)

x ljt 1 r r IJc

wd U =a r ' “ a ( r Ba U ) , 0<r<E , 0<z<z*

z 2k 2 r r U

(4.19)

wtfU = a r " m< M r md U u ) , 0<r<R , z <z<Z (4.20)

z 2Jc 2 r p 2,k

U =0 , z=0 t 0<r<R (4.21) lJc

d U =0 , r=0 , 0<z<z (4.22) r tJc

d U =0 , r =0 , z*<z<Z (4.23) r 2X

U +0 a T =0 , r = a , 0<z<z , e=l,2 (4.24) ejc k r &

- x ( a u + o e 2 t )+x f a u , + a a 2 t,) = 1 r l . k k r p 1 2 r 2 , k k r r 2

= w r Q O ’ , r =o , 0<z<z*

2 k

(4.25)

-X d U =p*?'U +£> lfjlz-z )-p(z-z ) ] , r=R , 0<z<Z (4.26)

2 p 21 21 i.-1 k

gdzie U =a T (r,z) , =d*> (T)/dT a d*=dO (z)/dz.

• Je ■ k k

Minimum funkcji F można poszukiwać jedną z metod opisanych w dodatku, przy czym wybór przybliżenia początkowego ({p}*)*°^ dla wektora parametrów {p}^ nie powinien być przypadkowy. Do jego określenia można wykorzystać metodę zaproponowaną w rozdziale trzecim.

4.4. MODEL NUMERYCZNY

W celu stworzenia praktycznej możliwości realizacji obliczeń przekształcono układ równań i warunków (4.1)+(4.10) wprowadzając nową bezwymiarową zmienną dla każdej z faz w następujący sposób:

, 0<r<o , 0<z<z a

r -a R - a

, cr<r<R , 0<z<z

(4.27)

, 0<r<R . z <z<Z

• 39 -

W przyjętym nowym układzie współrzędnych układ równań i warunków (4.1)+(4.10) przyjmie postać:

w(a T - y a ’a !a T ) =a o Zy ”d (y"a T ) , 0<y<l , 0<z<z* (4.28)

z I y l 1 y y 1

w[a T -(1-y)(R-a) 'o'a T l =

1 (4.29)

a (R-a) [y(R-a)+a] "a {[y(R-a)+alma T } , 0<y<l , 0<z<z*

2 y y 2

w 3 T = a R ' V " d ( y " d T ) , 0<y<l , z*<z<Z (4.30)

z 2 2 y y 2

T i=T , z=0 , 0<y<l (4.31)

d T =0 , y=0 , 0<z<z

y t (4.32)

ayT2=0 , y=0 , z <z<Z (4.33)

Tt (1 f z ) =T2 (0 , z ) =T , 0<zSz* (4.34)

-Xia'1ayT i(l,z)+X2(R-a)‘tdyTJ(0,z)=wy2Q o ’, 0<z<z* (4.35)

? S a(C ) . n = l ,N (4.36)

n n

-X2 (R-a) a yT 2 (l,z)=pp , 0<z<z - X 2R ' 1avT 2 (l,z)=pf> , z*<z<Z

(4.37)

Jak łatwo zauważyć, w nowym układzie współrzędnych obszary wlewka odpowiadające fazie ciekłej i fazie stałej są prostokątne. „Unieruchomiona została również granica rozdziału faz, a funkcja o = o {z) opisująca jej położenie stała się współczynnikiem równań.

Przekształceniu ulegną również i równania wraź 1iwościowe (4.18+4.26).

W nowym układzie współrzędnych mają one postać:

(4.38)

w ( a U , - ( 1 - y ) z 2 ,k ( R - o ) * V a v 2,kU ) =

=a2(R-a) 2[y(R*-a)+o) m<?y{ [y(R~a) + o ] m d ^ U 2 , 0<y<l , 0<z<z*

wali =a R \ md (y™a U ) , 0<y<l , z*<z<Z

z 21 2 V y 2Jt (4.40)

z=0 , 0<y<l (4.41)

a U =0 , y=0 , 0<z<z

» u

* (4.42)

a U =0 , y=0 , z*<z<Z

y 2,k (4.43)

UiA+eko'ła T =0 , y=l , 0<z<z* (4.44)

U +<J (R-o)"'aT=0 , y=0 , 0<z<z*

2 X k y 2 (4.45)

(4.46)

=vr3Q , 0<z<z*

- 41 -

oraz

-X2(R-cj) 1ayU2k=pł3’U2k+f> [(j(z-zk i)-fj(z-zt) ] , y=i , 0<z<Z (4.47)

Do rozwiązywania występujących w procesie iteracyjnyra zagadnień brzegowych można wykorzystać metody różnicowe [4.9,4.17]. Dla układu równań (4.28)+(4.35),(4.37) układ rozwiązujący wyprowadza się opierając się na równaniach (4.28)+(4.34) oraz warunku (4.37), natomiast warunek (4.35) wykorzystuje się do budowy procesów iteracyjnych, w których wyznacza się wartości funkcji a określające na każdym poziomie siatki różnicowej w starym układzie współrzędnych położenie frontu krzepnięcia.

Podobnie w przypadku równań wraż 1iwościowych (4.38) + (4.47) układ rpzwiązujący tworzy się na bazie równań (4.38)+(4.45) i warunku (4.47), natomiast warunek (4.46), podobnie jak wyżej, wykorzystuje się do skonstruowania procesu iteracyjnego.

Na przekształcone obszary fazy ciekłej i fazy stałej nakłada się siatki różnicowe oj i O 2 o krokach Az=Z/J oraz Ay =1/1 i Ay =1/1

1 1 "2 2 odpowiednio, tzn. :

Równania różnicowe wyprowadza się na podstawie niejawnej gwiazdy cztero- punktowej, zgodnie z zasadami podanymi w [4.10,4.12,4.13,4.18].

Układ rozwiązujący dla zagadnienia Stefana (4.28)+(4.35),(4.37) tworzy

°i j = {.yi:yi=1Ayr j=o,j}

= { y t : y t= i A y 2 . i - o . i ^ i z :z =jAz, j = o , j }

j układów równań kolejno dla J=1,J o postaci:

i =0,Ij ( 4 . 4 8 )

( 4 . 4 9 )

t.J ^OJ

2.0 ^ . O + 2.0 ^ . 0 + l.J t - t . J l,J l.J

.2,0 -2,0 ' l . J 1*1.J = d2,(a)

"1.J i=1.1, (4.50)

( « ) f-.1 , ( O ) ,

°J-1 S 1.jL I,'2 -J

g(*>r4TJ,<*1-3T2,<

° 2 . j ( . l . J o . j 1 ,(

(4.51) i ) .J, (s)-)

l.J J

gdzie j* oznacza ten numer warstwy siatki rdZnicowej (podlegający wyzna­

czeniu), na ktdrej kortczy sie proces krzepnięcia, a ^ = a ^ s\j A z ) , T*’jS>=

T (,>(iAy , JAz), e = l ,2, jest numerem iteracji. współczynnik i odpowiednio wynoszą:

< r = °

u lOj"*)wAy|

2 i - l

2i 2 o <* )

2(m+l)a Az

„«■»- 1 4

O.J . ( a ) , 2 . 2

‘ w(Oj ) A y (

2a Az

C - 1 + , ¿o , 2 -at -

(c ) wAy^

1 . 2 a Az c ..u>--- L _

0.J , !

w(cr^ ) A y f

'Id

a ^ Az ‘ 2i+l '

(o‘* ’wAy* 2i 2 d (‘ >

J

i=l,It-2

- 43 -

= tî ^ ■ i=i,-2

t‘ +

a Az

1 ’ 2 1 r 1 1

*

. ( * ) . 2 , 2

( a i ) w A y t 2 1 1 - 2

g(«> = - ai.J

XjAz

w r2 Q A y t o * “ )

(<>

32.j

X^Az w r 2Q A y 2 (R-a‘° )

a.U >=

i.J

a 2Az ( i - i ) A y2( R — a * * } )+ ° j * >

w A y2( R - a * * > ) 2 i A y2 (R-£Jj* *

*

( l - i A y 2 ) (ôj* J- J-l

2(R- O j ' ł)Ay2

2.U)

A y 2 (R-Oj* *)

.

a.U)

A y 2 (R-Oj* ’)

<« )

a2A2 ‘ ( i - i ) A y2( R - o ‘* > ) + o ' * > ■ J w A y2( R - O j“ >) 2 i A y2 (R— c / s ’l + o ' ” }

a 2Az ' { i + f ) A y2(R-cr'J'> ) + o<0 'm"

w A y2( R - ö j * ’ ) 2 i A y ^ R - a' * 5 i + o ' * ’

2 ,<«>_

" U

a 2Az ( i + i ) A y2< R - a < * )+ a < * )

w A y2( R - a * s ’ ) 2 i A y2( R - a ^ * >) + a ‘‘ )

( l - i A y2 ) (a** >- o j _ , ]

2(R- a \ m*Ay2

J 2

, i = i , r 2 - i

t.J-1

a 2Az

w A y2 ( R - C T j ' ' ) 2

i A y 2 ( R - a < “ > ) + * < * >

A y2 ( R - O j * ' i + O j ' 5

( 1 — A y 2 ) (o'" '-Oj. , )

2 ( R -0‘' )) A y2

,2,(*)

d = p «?

i3»j pr j

przy czym py=p (jAz f c*>) f natomiast:

T2» ^“>-t

12 , j 00 , €0=1

, €0=0

Rozpoczynając proces iteracyjny, przyjmujemy, ±e o^ =0^ . Proces itera- cyjny kortczyray, jeżeli:

| Oj*"1' | < c ( 4 . 5 2 )

gdzie c jest zadana dostatecznie mała liczba» a jako wynik końcowy a^ można przyjąć wynik uzyskany w ostatniej iteracji.

JeZeli kontynuując obliczenia otrzymamy dla pewnego j', Ze o*<e*, gdzie c* jest dostatecznie mała stała dodatnia, wówczas możemy przyjąć, Ze zakończyło się krzepnięcie wlewka oraz Ze głębokoSó zalegania ciekłego metalu wynosi z*=j*Az. Dalej będzie odbywało się już tylko stygnięcie wlewka i układ rozwiązujący tworzy J—J* układów równań kolejno dla j=j*+l,J postaci:

bo.j = 1H

2(m + 1 )a Az

2

w R 2A y 2

2a Az

2

— , 1=1, I -1 w R 2A y 2 2

RAy,

c0.j

2 (m+1 )a^Az wR2A y 2

1. J w R 2A y 2

2i+ l

2i

, i=l,I2-l

C .2.J = °

< 1 = > 1=0>I2-1

" ■ 2 . J = P J * J

Podobne układy rozwiązujące otrzymamy również dla równań wrażliwościowych.

W szczególności dla j = l tJ* i k=i71c mamy:

IJc.taty.U.ł tj l-lj i,j i,j

tk.(.) ltk.U) .tjc.<«> — r CM UI*1J * i=0r I,

(4.54)

- 47 -

2,k,<«> 2,k,<«) 2,k,<«> 2,k, (a )

l . J *-t.J l.J l.J

(4.55)

jk'<^> =ok-<“> +g k'(s)fu1’t'(“>+3U1-k'<*>-4U1-k’(*

i

l

J

k.is) f 2J..U) 2 « B>1 k.(S)

2.J [ 1.J

0,1 2,1

(4 .5 6 )

i odpowiednio dla j=j +1,J oraz k=l,K:

U i Z +bU ° U +CM C = • 1=0'I2 <4 -5 7 >

gdzie U*JJ'(!’,=U*S* (iAy^, j Az ) , e=l,2 , ©k,<*>=6**>( j A z ) , s jest numerem iteracji, a współczynniki odpowiednio wynoszą:

a * * ™ = 0 o.j

- a Az 2Ł-1 m i(o -o )..

h H — b H ■

a - . , = 0

2 (m + 1 ) a Az b>*.(.>= 1+ L _

0,J wOjAy2

2a t Az

’ 1=1-1 ,-1 w o jA y i

1 1

OJ 2 2

w o jA y i

r V , a i „

I",

I 2 A 2

L a w A y 1i 2 i J l a ,

J

d! T m>

e k,(e) Q 1

V ---

X A z W «)

i j --- a J&y tvf'2o

X A z

,k.(«)

2 . J ---

(R-oj)Ay2w r 2Q

^ X 2A z ó ^ < 0 t

( R - a ()2wy

^ U>= °

- 49 -

2Jś,(s)_

a l.J

2,k,0=>

a Az

2 ■ ( i - i j A y ^ R - O j J + a j ' ( l - i A y 2 ) [ o } - o ¡ t )' w A y2( R - ö j )

2

j

2,k,(s>

A y J ( R - o j )

b0 j - 1

a Az

1.J " 1+- w A y

2 - r ( i 4 )Äy2(R- ö j )+aji" +

“ (R-cr i ) 2 L lAy^R-e-jł+Oj J

wAy

L2Ai r ii+7 ) A y 2 ( R _ a j , + a jl’,1 ____

|(R-Ö j ) 2 [ l A y ^ R - O j i + O j J I

-1

2,k.(«>_

A y 2 (R-Oj) - p p ł

vr s

_ŁM*>=

'O.j

a 2Az (i+i-)Ay2( R - o J ) + o J (1— iA y2)(°j-°j-j) w A y2( R - ö j) 2 i A y2( R - ö j )+ oj Z { R - a ^ y \

,(* > D IT 2 ,2i,(*> J ° ' J

d # J -

R-CTj

d,k

2

2Jc.<«>_ u 2,k .-— -- di.J - ’ 1=1

J(S>= PJp ( j A z - z k_1)-fj(JAz-zk)j

d la J < j \ przy czym dV j=< ^ T (iAy#I j A z ), 0*1^=* T <iAy#t J A * > , e=l,2

*>’ - d*c>(T2 )/dT oraz:

J 12J 2

a 2,k= O o,J

a 2Az 2i-l *

w R 2 A y 2 [ 2i J

2,k

RAy,

2(m+l)a Az

„ « = 1 + L _

°'J w R 2A y 2

2a^Az bf* * l+-r--

4»J wR Ay

" 1 7 7 ~ pjt;

- SI -

.24c -o.j

2 ( m + 1 ) a 2 Az

vR2Ay2

a 2Az 2i+l »

°ij =

_{,J wR2Ay‘}

_{l 2 i J} 1 *

d!.J =

C <

= PJp ( j A z - z k i)-p(jAz-zk)

dla j>j'

4.5* PRZYKŁADY OBLICZEŃ

W celu zbadania skuteczności metody przeprowadzono szereg obliczeń testujących. W pierwszej kolejności testowano zaproponowaną metodę przybliżonego rozwiązywania zagadnień prostych Stefana, gdyż tylko poprawne funkcjonowanie metody rozwiązywania zagadnień prostych pozwalało oczekiwać wartościowych wyników przy rozwiązywaniu zagadnień odwrotnych. Polegało to na porównywaniu rozwiązań przybliżonych z rozwiązaniami analitycznymi

(przykład 4.1). Dalsze obliczenia testujące dotyczyły już zagadnień odwrotnych (przykłady 4.2 i 4.3), gdzie szczególną uwagę zwrócono na wrażliwość uzyskiwanych rozwiązań na zaburzenia w danych wejściowych.

Przykład 4.1

Jako jeden z przykładów rozważano obiekt o wymiarach R=Z=i, w którym położenie granicy rozdziału faz opisuje funkcja:

a — r 1-0.52

natomiast pola temperatur w obszarach odpowiednio funkcje

T = 300(r3+O.5rz-r )+500 or a z ;

T2 = 250 (r3+0.5rz—r )+500

Tak dobrane funkcje spełniają, równania (4.1) + (4.3) i warunki Stefana (4.7) + (4.8), gdy ra=0, w=i, T*=500, Xt=200, X2=100, a ^ a ^ l / 1 2 , r2= l , a 0=2800001^ (i~0.5z)3 . Nietrudno zauważyć, że funkcje te spełniają również warunki brzegowe:

T t=T** . r=0 , 0<z<Z

-X 3 T = a(T —T ) , r=R , 0<z<Z 2 r 2 2 03

oraz warunek początkowy:

1^=300 (r3-r) +500 , 0<r<R

gdy T**=500, a=-100 a T =0.

CD

Na przekształcone zgodnie z (4.27) prostokątne obszary odpowiadające poszczególnym fazom nakładano siatki różnicowe o różnej liczbie węzłów. Dla każdej pary siatek poszukiwano rozwiązania przybliżonego (pole temperatury, położenie granicy rozdziału faz) i porównywano z rozwiązaniem analitycznym.

Na rysunku 4.2 przedstawiono rozwiązanie rzeczywiste i wyniki, jakie uzyskano dla trzech różnych par siatek różnicowych, w tym dla ^510 *

°5.to’ n \ox 1 oraz " k « 1 nk « ’ przy czym rozwiązanie to dotyczy położenia granicy rozdziału faz.