29/44
Solidification of Metals and Alloys, Year 2000, Volume 2, Book No. 44 Krzepnięcie Metali i Stopów, Rok 2000, Rocznik 2, Nr 44 PAN – Katowice PL ISSN 0208-9386MODELOWANIE ZADAŃ Z OSTRYM FRONTEM KRZEPNIĘCIA Z WYKORZYSTANIEM II SCHEMATU MEB
J. MENDAKIEWICZ1, A. PIASECKA BELKHAYAT2, R. SZOPA3
1 Katedra Wytrzymałości Materiałów i Metod Komputerowych Mechaniki Politechnika Śląska
2 Instytut Matematyki i Informatyki, Politechnika Częstochowska
STRESZCZENIE
W pracy przedstawiono sposób modelowania procesu krzepnięcia zachodzącego w stałej temperaturze (problem Stefana), przy czym rozpatrywano zadanie 1D. Wykorzy- stano II schemat metody elementów brzegowych. Omówiono algorytm rozwiązania oraz pokazano przykład obliczeń numerycznych.
1. SFORMUŁOWANIE ZADANIA
Jednym z najważniejszych modeli matematycznych opisujących proces krzepnięcia i stygnięcia metalu jest sformułowany ponad 100 lat temu problem Stefana. W zadaniu tym rozpatruje się półprzestrzeń ograniczoną płaszczyzną o temperaturze Tb, w której wyróżnia się dwa zmienne w czasie podobszary, a mianowicie podobszar fazy stałej (x∈(0, η)) oraz podobszar cieczy (x∈( η, ∞)). Temperatura powierzchni kontaktu x = η jest równa Tkr (temperatura krzepnięcia). Dla czasu t= 0 temperatura w całym obszarze wynosi T(x, 0) = T0. Również dla x → ∞ temperatura cieczy jest równa T0. Problem ten posiada rozwiązanie analityczne [1], które będzie wykorzystane do przetestowania dokładności i efektywności metody elementów brzegowych w tzw. zadaniach z ruchomymi brzegami, przy czym zastosowany zostanie II schemat MEB przy założeniu, że obszar odlewu jest obszarem ograniczonym (x∈( η, G)) i na granicy x=G obowiązuje warunek brzegowy qL (G, t)=0.
Nieustalone pole temperatury w rozważanym obszarze opisuje układ dwóch równań różniczkowych w postaci
1 Dr inż., e-mail: mendoz@rmt4.kmt.polsl.gliwice.pl,
2 Dr inż. e-mail: alicja@rmt4.kmt.polsl.gliwice.pl
3 Dr hab. inż., prof. P.Cz.,
224
2 2
2 2
) , ( )
, : (
) , ( )
, : (
0
x t x a T
t t x G T
x
x t x a T
t t x x T
L L L
s S S
∂
= ∂
∂
< ∂
<
η
∂
= ∂
∂ η ∂
<
<
(1)
Dla x = η przyjmuje się warunek brzegowy
⎪⎩
⎪⎨
⎧
=
=
= η
∂ λ ∂
∂ − λ ∂ η
=
kr L
S
L V S L
S
T t x T t x T
L t x
t x T x
t x T x
) , ( ) , (
d ) d , ( )
, (
: (2)
lub
(3)
⎩⎨
⎧
= η
= η
= η
− η
kr L
S
V S
L t T t T
T
wL t q t q
) , ( ) , (
) , ( ) , (
Ponadto dla x=0 i x=G
0 ) , ( , ) , 0
( t =T q G t =
TS b L (4)
oraz dla t=0
(5)
kr
L x T T
T ( ,0)= 0 ≥
W równaniach (1), (2) as, aL oznaczają współczynniki dyfuzji ciepła fazy stałej i ciekłej, λS, λL - współczynniki przewodzenia ciepła, LV - utajone ciepła krzepnięcia, natomiast w=dη/dt jest prędkością krzepnięcia.
2. METODA WYZNACZANIA POŁOŻENIA FRONTU KRZEPNIĘCIA
Obszar płyty [0, G] dzielimy na n elementów liniowych, co oznacza, że na odcinku [xj-1 , xj ], j=1, 2, ..., n, pole temperatury przybliżamy funkcją liniową. Tak więc siatkę geometryczną tworzą punkty
G x x
x x
x < < < j < j < < n =
= 0 1 Κ −1 Κ
0 (6)
Załóżmy, że w chwili czasu t=t f-1 współrzędna geometryczna frontu krzepnięcia jest równa xj-1 oraz, że po czasie ∆t front krzepnięcia „przemieści się” do punktu xj . Czas ∆t odpowiadający temu przesunięciu (xj-1→ xj ) nie jest znany. W chwili t + ∆t=t f węzły x0, x1, ..., xj-1 należą do obszaru fazy stałej, natomiast xj+1, xj+2, ..., xn - do fazy ciekłej.
Kolejne dwa położenia frontu krzepnięcia pokazano na rysunku 1. Na początku obliczeń dotyczących przejścia od chwili t f-1 do t f zakładamy pewien krok czasu ∆t
225
(wynika z niego bezpośrednio chwilowa prędkość krzepnięcia w = ∆x/∆t) i przyjmując dla x= xj : TS (xj , t f )= TL (xj , t f ) =Tkr sprowadzamy problem do dwóch oddzielnych (rozprzężonych) zadań z warunkiem I rodzaju na granicy rozdziału faz, to znaczy dla x=xj
Rozwiązania tych zadań zawierają informacje o brzegowych strumieniach ciepła na froncie krzepnięcia, co pozwala skorygować przyjęty pierwotnie krok czasu.
Wykorzystujemy tutaj pierwszą część warunku Stefana (3), którą zapiszemy w postaci
) , ( ) ,
( j f S j f
L
V
t x q t x q
x t L
−
= ∆
∆ (7)
W ten sposób metodą kolejnych przybliżeń (procedura iteracyjna jest szybko zbieżna) wyznaczamy czas przejścia frontu krzepnięcia od węzła xj-1 do węzła xj . Jako punkt startowy procesu iteracyjnego można przyjąć wartość ∆t z przejścia poprzedniego.
Rys. 1. Przemieszczający się front krzepnięcia.
Fig. 1. Moving solidification front.
3. II SCHEMAT MEB
Istotę algorytmu nazywanego II schematem metody elementów brzegowych omówimy na przykładzie obszaru jednorodnego o stałych granicach [0, L]. Temperaturę w tym obszarze oznaczymy przez T (x, t ), natomiast współczynnik dyfuzji ciepła i współczynnik przewodzenia przez a i λ. Taki opis II schematu jest całkowicie wystarczający, ponieważ zaproponowana metoda obliczeń sprowadza się do
„rozprzężenia” odlewu na jednorodne podobszary cieczy i ciała stałego, przy czym na ich granicach zadana jest temperatura lub brzegowy strumień ciepła.
Do rozważań wprowadzamy siatkę czasu
226
(8)
∞
<
<
<
<
<
<
<
<
<
<
<
=t0 t1 t2 Κ ts−1 ts Κ tf−1 tf Κ tF 0
ze stałym krokiem ∆t.
W przypadku II schematu MEB dla przejścia t f-1 → t f otrzymuje się następujące równanie całkowe [2]
∫
∫
∫
ξ + ξ
⎥ +
⎥
⎦
⎤
⎢⎢
⎣
⎡ ξ
= ρ
⎥ =
⎥
⎦
⎤
⎢⎢
⎣
⎡ ξ
+ ρ ξ
∗
−
=
=
∗
=
=
∗
−
−
L
f f
L x
x t
t
f
L x
x t
t
f f
x t x T t t x T t
Z t
t x T t t x c q
t t x q t t x c T
t , T
f
f
f
f
0
0 0 1
0 0
d ) , ( ) , , , ( ) , ( d
) , ( ) , , , 1 (
d ) , ( ) , , , 1 (
) (
1 1
(9)
przy czym Z(ξ, t f-1 )=0 dla f=1, natomiast dla f=2, 3, ..., F:
∑ ∫
∑ ∫
−
=
=
=
∗
−
=
=
=
∗
−
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡ ρ ξ
−
⎥ +
⎥
⎦
⎤
⎢⎢
⎣
⎡ ρ ξ
= ξ
−
−
1
1 0
1
1 0
1
1 1
d ) , ( ) , , , 1 (
d ) , ( ) , , , 1 (
) , (
f
s
L x
x t
t
f f
s
L x
x t
t
f f
s
s s
s
t t x q t t x c T
t t x T t t x c q
t Z
(10)
W równaniach (9), (10) ξ∈(0, L) oznacza punkt obserwacji, natomiast T * jest tzw.
rozwiązaniem fundamentalnym (podstawowym) i dla omawianego zadania ma postać
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡
− ξ
− −
− π
=
∗ ξ
) ( 4
) exp (
) ( 2 ) 1 , , , (
2
t t a
x t
t a t
t x
T f f
f (11)
Strumień ciepła wynikający z rozwiązania fundamentalnego definiuje się następująco
x t t x t T
t x q
f f
∂ ξ λ∂
−
=
ξ ∗
∗ ( , , , )
) , , ,
( (12)
Rozważać będziemy tzw. stałe elementy po czasie, czyli przyjmiemy założenie, że dla t∈[t f-1, t f ]: T(x, t)=T(x, t f ), q(x, t)=q(x, t f ), oraz dla t∈[t s-1, t s ]: T(x, t)=T(x, t s ) i q(x, t)=q(x, t s ).
Wprowadzamy następujące oznaczenia
227
∫
− ξ= ρ
ξ ∗
f
f
t
t
f t t
t x c T
x g
1
d ) , , , 1 (
) ,
( (13)
∫
− ξ= ρ
ξ ∗
f
f
t
t
f t t
t x c q
x h
1
d ) , , , 1 (
) ,
( (14)
∫
− ξ= ρ
ξ ∗
s
s
t
t
f
s T x t t t
x c g
1
d ) , , , 1 (
) ,
( (15)
∫
− ξ= ρ
ξ ∗
s
s
t
t
f
s q xt t t
x c h
1
d ) , , , 1 (
) ,
( (16)
Powyższe całki wyznacza się analitycznie [2]. Zmierzamy z punktem obserwacji ξ do brzegów obszaru (ξ → 0+ oraz ξ → L¯ ) i otrzymujemy układ dwóch równań
(17)
⎥⎥
⎦
⎤
⎢⎢
⎣ +⎡
⎥⎥
⎦
⎤
⎢⎢
⎣ +⎡
⎥+
⎥⎦
⎤
⎢⎢
⎣
⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡
−
−
−
= −
⎥⎥
⎦
⎤
⎢⎢
⎣
⎥⎡
⎦
⎢ ⎤
⎣
⎡
−
−
−
−
−
−
+ +
) , (
) , 0 ( ) , (
) , 0 (
) , (
) , 0 ( 1 ) , ( ) 0 , (
) , 0 ( 1 ) 0 , 0 ( )
, (
) , 0 ( ) , ( ) 0 , (
) , 0 ( ) 0 , 0 (
1 1 f f f
f
f f f
f
t L Z
t Z t
L p
t p
t L T
t T L L h L
h
L h h
t L q
t q L L g L g
L g g
przy czym dla f =1: Z(0, t 0 ) = 0, Z(L, t 0 ) = 0, natomiast dla f = 2, 3, ..., F:
(18)
∑
∑
−
=
−
= − −
+ +
−
−
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡
−
− −
⎥+
⎥⎦
⎤
⎢⎢
⎣
⎡
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡
−
= −
⎥⎥
⎦
⎤
⎢⎢
⎣
⎡
1
1 1
1 1 1
) , (
) , 0 ( ) , ( ) 0 , (
) , 0 ( ) 0 , 0 (
) , (
) , 0 ( ) , ( ) 0 , (
) , 0 ( ) 0 , 0 ( )
, (
) , 0 (
f
s s
s s
s
s s
f
s s
s s
s
s s
f f
t L q
t q L L g L g
L g g
t L T
t T L L h L h
L h h
t L Z
t Z
oraz
(19)
∫
ξ=
ξ ∗
L
f
f T xt t T xt x
t P
0
0 0) ( , )d ,
, , ( ) , (
Po rozwiązaniu układu (17), znane są wartości brzegowe temperatur i strumieni ciepła i na tej podstawie wyznacza się temperatury w zbiorze punktów wewnętrznych.
228
4. PRZYKŁAD OBLICZEŃ
Analizowano krzepnięcie płyty o grubości 0.04 [m] wykonanej z miedzi. Parametry termofizyczne przyjęto zgodnie z [1]. Założono Tb =1000o C, Tkr =1083o C, T0 =1100o C.
Na rysunku 2 przedstawiono rozwiązanie uzyskane za pomocą II schematu metody elementów brzegowych dla czasów 1 - 0.31 [s], 2 - 1.09 [s], 3 - 2.29 [s], 4 - 5.87 [s] oraz 5 - 11.13 [s]. Należy podkreślić, że zgodność otrzymanych wyników z rozwiązaniem analitycznym dotyczącym półprzestrzeni jest w pełni zadowalająca.
Rys. 2. Rozwiązanie numeryczne.
Fig. 2. Numerical solution.
LITERATURA
[1] W. Longa, Krzepnięcie odlewów w formach piaskowych, Śląsk, Katowice, 1974 [2] E. Majchrzak, Zastosowanie metody elementów brzegowych w termodynamice
procesów odlewniczych, Wyd. Pol. Śl., Mechanika, 102, Gliwice, 1991 MODELING OF THE STEFAN PROBLEM USING THE BEM SUMMARY
The numerical model of 1D Stefan problem is solved using the 2nd scheme of the BEM.
The theoretical background and also the example of numerical simulation are presented.
Reviewed by prof. Stanisław Jura
