CC BY
36
УДК 544.3 + 519.63 + 517.956.4
DOI: 10.14529/met160406
МАТЕМАТИЧЕСКАЯ МОДЕЛЬ КИНЕТИКИ ФАЗОВОГО ПЕРЕХОДА ПРИ НАГРЕВЕ ПОВЕРХНОСТИ ЦИЛИНДРИЧЕСКОГО ТЕЛА
11 2 А.Д. Дрозин, Н.М. Япарова, В.Я. Гольдштейн2
1 Южно-Уральский государственный университет, г. Челябинск, 2ЗАО «Ферросплав», г. Челябинск
В случае если в теле при его нагреве или охлаждении происходит фазовый переход, процесс распространения в нем тепла претерпевает значительные изменения. При переходе через температуру ликвидус энтальпия, плотность и коэффициент теплопроводности как функции температуры претерпевают разрыв. При нагреве поверхности твердого металлического цилиндрического тела выше температуры его плавления у поверхности возникает жидкая фаза и соответствующая граница раздела фаз, перемещающаяся к оси цилиндра. В статье приведен вывод уравнений распространения тепла в теле, учитывающий, помимо фазового перехода, изменение размеров участков тела вследствие зависимости плотности от температуры. К полученной системе уравнений был применен метод выпрямления фронтов - были введены такие координаты, относительно которых поверхность раздела фаз является неподвижной. Полученную систему дифференциальных уравнений сводили к конечно-разностным уравнениям. Для решения полученной системы разностных уравнений была разработана компьютерная программа. В статье приведены результаты одного из таких расчетов. Разработанная методика позволяет рассчитывать скорость движения границы раздела фаз, а также температуру тела в любой его точке и в любой момент времени. Результаты работы могут быть интересны не только металлургам, но и метрологам при разработке теории самотестирующегося датчика температуры.
Ключевые слова: теплопроводность; фазовый переход; модель.
Введение
Расчет скорости прогрева образца (деталей при термообработке, кусков шихты в начальный период ее плавления) значительно осложняется наличием фазового перехода. На поверхности раздела фаз терпят разрыв все термические коэффициенты, выделяется теплота фазового перехода [1-3], граница раздела фаз движется. В данной работе приведен метод расчета кинетики прогрева первоначально твердого образца цилиндрической формы.
1. Вывод уравнения теплопроводности
Вначале рассмотрим линейный случай. Пусть среда представляет собой некоторую трубу с постоянной площадью поперечного сечения F и с теплоизолированной боковой поверхностью и
пусть ось Ох направлена вдоль трубы (рис. 1).
Выделим некоторый участок [х;х + Дх] и составим уравнение баланса энтальпии в нем за время х ? + А?]. В момент времени / выделенный объем
имеет энтальпию ДxFр(t, X) h (?, X) , где X - некоторая
точка (х < X < х + Дх) . В момент времени ? + Д? выде-
ДxFр(t + Д?,х)И(? + Д?,х) .
х+Дх
Рис. 1. К выводу уравнения баланса энтальпии
ленный объем имеет энтальпию
где х - некоторая точка
(
х < х < х +
Дх ). (Мы применили теорему о среднем значении функции). Изменение энтальпии
выделенного участка составит
ДЯ =
р (t + hi, х) h (t + hi, х) - р (t, х) h (t, х)
F Дх.
Через площадку с координатой х за время [?, ? + Д? ] вследствие теплопроводности поступит энтальпия q (7, х) FД?, где q (?, х) - плотность теплового потока (имеет размерность Дж/ м с),
Дрозин А.Д., Япарова Н.М., Математическая модель кинетики фазового перехода
Гольдштейн В.Я. при нагреве поверхности цилиндрического тела
где / - некоторая точка (t < Г < t + Д/) . Через площадку с координатой х + Дх за это время уйдет энтальпия д (Г, х + Дх) FДt, где / - некоторая точка (/ < t < t + Дt). Изменение энтальпии за счет теплопроводности составит
ДН ,=
д (Г, х)- д (Г, х + Дх )
Е Д/ .
Вследствие изменения температуры изменяется плотность системы, а следовательно, и ее объем, в нашем случае, длина. Материальные точки системы движутся. Пусть w{t, х) - скорость материальной точки с координатой х в момент времени t. За время Д/ через левую границу выделенного объема зайдет новая масса, равная т1 = Ер(Г,х)w(t ,х)Д/, где / - некоторая точка
(/ </ < / + Д/). Эта масса принесет энтальпию Ер(/,х)w(х)h(/,х)д/, где / - некоторая точка
(/ < / < / + Д/). За это же время через правую границу выделенного объема уйдет масса, равная т2 = Ер(/ ,х + Дх)w(t ,х + Дх)Д/, где / - некоторая точка (/ </ < / + Д/). Эта масса унесет энтальпию Ер(/ ,х + Дх)w(/,х + Дх)h(/,х + Дх)д/, где / - некоторая точка (/ < / < / + Д/). Изменение энтальпии за счет массопереноса составит
ДНт = Е
Е Д/.
(/ , х | w (/ , х | И (/, х )-р(/, х + Дх | w (/' , х + Дх | И (/' , х + Дх ) Составим уравнение баланса энтальпии: ДН = ДН^ + ДНт . Получим
Е Д/ +
р(/ + Д/,х)И(/ + Д/,х)-р(/,х)И(/,х) ЕДх = д(Г,х)-д(Г,х + Дх) +Е |р (/ , х ) w (/, х ) И (/ , х ) - р (/, х + Дх ) w (/ , х + Дх ) И (/' , х + Дх )
р(/ , х) w (/, х) И (Г , х )-р(/ , х + Дх) w (/ , х + лх | И (/ , х + лх | Е Д/ или
р (/ + Д/, х) И (/ + Д/, х) - р (/, х) И (/, х) д (Г, х + Дх) - д (Г, х) Д/ Дх
р (/ , х + Дх) w (/ , х + Дх ) И (/ , х + Дх) - р (/ , х) w (/, х ) И (/, х )
Дх
Переходя к пределам при Д/ ^ 0, Дх ^ 0 , получим
Я/(р(/,х)И(/,х)) = -£(д(/,х))-£-(р(/,х)w(/,х)И(/,х)) . (1.1)
Для преобразования этого уравнения составим дополнительно уравнение баланса массы выделенного участка.
В момент времени / выделенный объем имеет массу ДхЕр (/, хс) . В момент времени / + Д/ выделенный объем имеет массу ДхЕр (/ + Д/, хс) . Изменение массы выделенного участка составит
Дт =
р( / + Д/, х) - р (/, х)
Е Дх.
За время Д/ через левую границу выделенного объема поступит масса, равная Ер(/,х)w(/,х)Д/, а через правую границу уйдет масса, равная Ер(/ ,х + Дх)w(t,х + Дх)Д/. Изменение массы составит
Дт = Е
р(/, х) w (/, х) - р(/ , х + Дх) w (/ , х + Дх )
ЕД/.
Составив уравнение баланса массы и перейдя к пределам при Д/ ^ 0, Дх ^ 0 , получим
я я
Ъ(р(х)) = -&(р(х)w(/,х)) . а.2)
Мы получили систему уравнений (1.1), (1.2): Переписав уравнение (1.1) в виде
dh , р— + h Ct
Sp + C(pw)
Ct СХ
\
Cq Ch
pw- , Cx Cx
и используя уравнение (1.2), получим
дh дq дh
Р^7 = -я—Р^. (13)
дt дх дх
Так как к = к (Т) является функцией только температуры и состава системы, то к зависит от
I, х только через температуру. Поэтому
дк dк дТ дк dк дТ „.
— =--, — =--. (1.4)
дt dT дt дх dT дх
^ dк / ч dк
Производная — является теплоемкостью (Ср ) : сР = —, которая также зависит только от dT dT
температуры и состава системы.
Полагая, в соответствии с первым законом Фурье,
дТ
q ^, х )=-х(Т )&' (1-5)
где X - коэффициент теплопроводности, с учетом (1.4), (1.5) получим
дТ д ( дТ \ дТ
дt дх ^ дх) дх
Был рассмотрен линейный случай. Общий, трехмерный случай имеет вид
рсР дТ = div ( XgradT ) - р жр gradT . (1.7)
2. Распространение тепла на участке с фазовым переходом
Найдем условия на движущейся границе раздела фаз. Пусть в момент времени t граница была на плоскости с координатой ^ (t), а в момент времени t + Дt - на плоскости с координатой
+ (рис.2). Заметим, что Т7^,^)) = -тем-
пература фазового перехода. Составим баланс энтальпии на участке ;^ + Д)] за время t + Д ].
В момент времени t этот участок был занят твердой фазой, и его энтальпия была равна
Н ^ ) = 5 )(т (^ х( 5))) к(5 )(т (t, х( 5))). В момент
времени t + Дt этот участок будет занят жидкой фазой и его энтальпия будет равна Н ^ + Дt) = ^Д^р ' (т (t + Дt, х(ь))) к(ь) (т ^ + Дt, х( ь))). Здесь
ш) ш+до
Рис. 2. К выводу уравнения баланса энтальпии на участке раздела фаз
%(t)< Х(L)< % (t + At), A% = % (t + At) - % (t) . Изменение энтальпии составит
AH = F A%
(L )(r (t + A, Х( L))) h( L )(r (t + A, Х( L )))-p( 5 )(r (t, Х( 5))) h( S )(r (t, Х( S)))
Через площадку с координатой х = %(t) за время [t, t + Л ] поступит энтальпия
q(L)(f, %(t)) FAt. Через площадку с координатой x = %(t + At) за время [t, t + At] уйдет энтальпия
Л 5)
(t, %(t + At))
) SAt. Здесь t < t < t + At, t < t < t + At.
P
Дрозин А.Д., Япарова Н.М., Математическая модель кинетики фазового перехода
Гольдштейн В.Я. при нагреве поверхности цилиндрического тела
Изменение энтальпии за счет теплопроводности составит
ДН ,=
д(")(I ^))-д( 5 )(Г, ^ + Д/))
Е Д/.
За время Д/ через левую границу выделенного объема зайдет новая масса, равная т1 = Ер((/,))) w(, ))Д/, где / - некоторая точка (/ < / < / + Д/). Эта масса принесет
энтальпию Ер(Ь)(г (/,£(/))) w(Ь)(/, £(/))И(Ь)(г (/,£(/))) Д/, где Г - некоторая точка (/ < Г < / + Д/). За это же время через правую границу выделенного объема уйдет масса, равная т2 = Ер(5^ (т (/,+ Д/))) w(5)(/ , + Д/))Д/, где / - некоторая точка (/ < / < / + Д/) . Эта масса
унесет энтальпию Ер(5)(т(/,+ Д/)))w(5)(/ ,+ Д/))И(5^(т(/ ,/ + Д/)))Д/, где / - некоторая
точка (/ < / < / + Д/) . Уравнение баланса энтальпии в выделенном объеме имеет вид:
Е Д^Гр( ^(т (/ + Д/, х( *))) И( * )(т (/ + Д/, х( * )))-р( 5 )(т (/, х( 5))) И( 5 )(т (/, х( 5 )))" = = д( ")(Г, ^))-д( 5 )(Г, ^ + Д/))] е д/ + е [р( ^(т (Г, $(/))) w(*)(/, $(/)) И(* )(т (/, $(/)))--р( 5 )(т (/£(/+Д/))) w( 5)(/£(/+Д/)) И( 5 )(т (/£(/+Д/)))]д/.
Перейдем в этом равенстве к пределам при Д/ ^ 0 . Учитывая, что при этом + Д) ),
), х(5) ^ ), Д^ , Г ^ /, Г ^ /, / ^ /, Г ^ /, / ^ /, / ^ /, получим "р( ^)(т(5 И( £)(т(5 ^ ))-р( 5 )(т(5 И( 5 )(т(5 ^^ = д( , £(/)))-д( 5)(/, £(/)) +
+р( * )(т(5 ^)) w( % £(/)) И( ^(т(5 ^)-р( 5 )(т(5 ^)) w( 5)(/, £(/)) И( 5 )(т(5 -Учитывая, что
И( ^(т5 ^ )-И( 5 )(т5 ^ ) = ДИ( 5 - ^;
р№ (т5^ ) - р(5 (т5^ ) = Др(5^), (2.1)
где ДИ(5 ^ - удельная теплота фазового перехода Др(5 - изменение объема при плавлении, из предыдущего уравнения получим
И( 5) (т5 ^ ) Др( 5 ^ + ДИ( 5 - Ь) (т5 ^ )] ^ = д(Ь) ((/, £(/))) - д( 5) (/, £( /)) +
+р( * )(т(5 ^)) w( ь)(/, £(/)) И( ^(т(5 ^)-р( 5 )(т(5 ^)) w( 5)(/, £(/)) И( 5 )(т(5 -Учитывая первый закон Фурье (1.5) и опуская аргументы, найдем
(И(5)Др(5^) + ДИ(5^щ^ = 5) ^ - ^) ^ + р(^w(i)И(Ь) - р(5)w(5)И(5). (2.2)
V > Ж Ях Ях
В момент времени / этот участок был занят твердой фазой, и его масса была равна т (/) = ЕД^р(5)(т (/, х(5))). В момент времени / + Д/ этот участок будет занят жидкой фазой, и
его масса будет равна т(/ + Д/) = ЕД^р(^(т(/ + Д/,х(^))И. Здесь £(/)<х(ь),х{ь)<^(/ + Д/),
Д£ = £(/ + Д/) -^(/) . Изменение массы составит Дт = 5Д^ р((т(/ + Д/,¿с())- р(5) (т(/,х(5))) . За время Д/ через левую границу выделенного объема зайдет новая масса, равная
Физическая химия и физика металлургических систем_
т = Fр(i)(т (/, ^(t))) ^(t , где t - некоторая точка (t < t < t + Дt). За это же время
через правую границу выделенного объема уйдет масса, равная т2 = Fр(5)(Т(7,^(t + Дt)))х
х w(S) (t , %(t + At))At, где t - некоторая точка (t < t < t + At) . Изменение массы составит
Am = F
(L )(Г (t, %(t))) w( L)(t, %(t ))-p( S )(T (t, %(t + At))) w(S)(?, %(t + At))'
At.
Уравнение баланса массы в выделенном объеме имеет вид:
F
= F
р( L )(T (t + At, Х( L )))-p( S )(T (t, Х(S))) p( L) (t (t, %(t))) w( L )(t, %(t ))-p(S) (t (t, %(t + At))) w(S )(f, %(t + At))
At
или
p( L)(t (t + At, Х(L)))-p( S )(T (t, X(S)))
At
= p( L)(T (t, %(t))) w( L)(t, %(t ))-p( S )(T (f, %(t + At))) w( S )(f, %(t + At)).
Перейдем в этом равенстве к пределам при At ^ 0 и получим p(
(L)(t(S-L))-p(S)(t(S-L)) = p(L)(t(S-L))w(L)(t,%(t))-p(S)(t(S-L))w(S)(t,%(t)).
pv
S - L'
(t,%(t)). Так как
p(L (ts l ) - p(S (ts l ) = Ap(S L), где Ap(S L^ - изменение объема при плавлении, из предыдуще-
го уравнения найдем
= p(
Ap( S -L) = p( L )(T(S -L)) w( L )(t, %(t ))-p(S )(T(S -L)) w( S )(t, %(t))
или, опуская аргументы,
Ap( S - L) = p( L) w( L )-p( S) w( S). F dt У У
Подставляя (2.3) в (2.2), получим
Ah(S-L)p(L) = S)CT - X(L) CT + p(L)w(L) (h(L) - h(S)) .
dt cx Cx v '
(2.3)
или
Ah( S -L )p(L) = S) — - X(L) — + p(L) w(L)Ah( S -L). dt Cx Cx
(2.4)
3. Цилиндрическая задача для системы с движущейся границей раздела фаз 3.1. Постановка задачи для системы с движущейся границей раздела фаз
Рассмотрим цилиндрическое тело бесконечной длины, нагревающееся по произвольному закону с боковой поверхности (рис. 3). Упростим, в первом приближении задачу. Не будем учитывать движение точек среды.
В цилиндрических координатах при осевой симметрии gradT = —, ДТ = ——| г — | [4]. Переписав, с уче-дг г дг \ дг )
том этого, уравнения (1.7), (2.4), получим следующую систему уравнений:
Рис. 3. Схема прогрева цилиндрического образца
^f=И И" f)0 < - <«); (3.1)
(LCL)CT=1 Afx(LV—1, %(t)
P dt r Cr У Cr) w
|<r <R;
P
Дрозин А.Д., Япарова Н.М., Математическая модель кинетики фазового перехода
Гольдштейн В.Я. при нагреве поверхности цилиндрического тела
т (/, 5( / )) = т(5 -Ь ]; (3.3)
ДИ( 5 -Ь )р(Ь) Ъ = 5 ) Ят (t, 5^)-0) - ¿(Ь ) Ят (t, 5^) + 0) ;
ят (/,0)
р( Ь^ = х( 5)-^^->--Х(Ь)-^ ^ ' + 1 '; (3.4)
& Яг Яг
= 0, т(/,Я) = ц(/), / > 0; (3.5)
Яг
5(0) = Я, т(0,Г) = ф(Г), 0 <г <Я. (3.6)
3.2. Применение метода выпрямления фронтов для системы
с движущейся границей раздела фаз
Граница раздела фаз Ь и 5 является движущейся. Поэтому будем строить разностную схему математической модели методом выпрямления фронтов [5]. Для этого произведем замену переменной г на новую независимую переменную г г-5( /)
х = —— при 0 <г <5(/) , х =-^-4 +1 при 5(/)<г <Я . (3.7)
5(/) Я-5(/)
В результате подвижная граница г = 5 (/) областей действия уравнений математической модели заменится неподвижной границей х = 1, а тоже подвижная (из-за изменения плотности) граница г = Я (/) заменится неподвижной границей х = 2 .
Краевая задача (3.1)-(3.6) сведется к виду
1. Уравнение теплопроводности твердой фазы (0 < х < 1) :
Ят ' АГхЮхЯПДхвт.
Я/ р(5)с{5)х52 Ях^ ~ ЯххЯх ' (3'8)
2. Уравнение теплопроводности жидкой фазы (1 < х < 2) :
Ят 1 1 ЯГ (ги , ч/ ччЯт ^ Я'х + 5'(1 - х) ят
— =-7-Цт7-ч- , Л , Г--—I ¿(Ь)(5 + (Я-5)(х-1))— 1 +--—. (3.9)
Я/ 5+(Я-5)(х-1) р(Ь)с(Ь)(Я-5)2 Ях^ ^ ^ Ях) я-5 Ях
3. Первое условие на границе раздела фаз (х = 1):
Г (/,1) = Г(5 -Ь). (3.10)
4. Второе условие на границе раздела фаз (х = 1):
= ¿(5) 1 ят (/,1 - 0) ¿(Ь) 1 ят (/,1+0) (311)
Ж/ "дй(5-Ь)р(Ь) 5 Ях ДИ(5-Ь)р(Ь)(Я -5) Ях • ( - )
5. Условие на левой границе твердой фазы (х = 0):
ят (/,0)
—= 0. (3.12)
Ях
6. Условие на правой границе жидкой фазы ( х = 2) :
Г(/,2) = ц(/), />0. (3.13)
7. Начальное условие движущейся границы раздела фаз (/ = 0):
5(0) = Я (Я0 <я). (3.14)
8. Начальное распределение температур (/ = 0):
т(0,х) = ^Ф(Я»х)'°£х £1' (3.15)
1ф(Я +(Я - Я )(х -1)), 1 < х < 2.
3.3. Переход к разностным уравнениям для системы
с движущейся границей раздела фаз
Перейдем от дифференциальных уравнений к разностным [6]. Зададимся некоторым натуральным числом nc и разобьем отрезки 0 < x < 1 и 1 < x < 2 каждый на nc равных частей точками x0 = 0, x1,...,x^c = 1, x c+1,..., x^c = 2 с длиной каждого отрезка l = \/nc. Таким образом,
xj = jl, j = 0,...,2nc.Ось времени t > 0 разобьем на отрезки точками t0 = 0,t1,t2,... с шагом т , который может быть переменным.
Образованная сетка состоит из узлов (x j, tk ), j = 1,2,...,2nc; к = 0,1,2,...
Решать разностную задачу будем с помощью неявной четырёхточечной схемы. Процесс решения состоит в переходе от к -1 -го временного слоя (соответствующего t = tk_1), на котором все функции уже определены, к новому к -му временному слою. На каждом этапе перехода будем обозначать функцию y(t,x) в узле (x;-,tk-1) символом , а в узле (x;-,tk ) - символом ,
не выписывая явно временной индекс.
Аппроксимируем частные производные следующими соотношениями:
дТ (
{fk, xj)'
dx (
x dx l dx
tj — tj .
? X
Tj+1 — TJ
2l
J. x=xj,
X ,-j.i /о x
Tj +1 " Tj
— X ,_,/n,x
T — T—11
j+1/2 j+1/2 l A"j—1/2xj—1/2 I
(3.16)
(3.17)
(3.18)
t=tk
где
rj—12
ri —1 + rj
Cj+V2
rj + rj +1
X
j—12
X (T(^, rj—1 )) + X (T(tk, rj)) X (T(tk, rj )) + X (T(tk, rj+1))
X _
2 ' j+1/2 2 Для односторонних производных в (3.17) используем аппроксимации
at(s) T(S} - T(S) яг(L) T
—(tk,1 — 0) ——-n!—l, —(t„1 + 0)~-П
(tk ,1 + 0)'
(L) — T(L) nc +1 nc
дх I дх I
Выпишем полную систему разностных уравнений математической модели. 1. Основные уравнения твердой фазы:
а) to =
1
0_ гр (s Ше212+1^
(S) (S р0 CP,0s
4X j )х
T1 +
p0S ^s212
4X j )x
Г p0S )cPSS0e212
t0
4Xj )x
+1
б) j xj—12 — 2 P((S j2**') tj—1 —
f P(S)C(S)l2 >
X (S) + X (S) , P j CP, jl s2
X j —12 xj —12 + X j+12 xj+12 + _ xj S
(3.19)
(3.20)
(3.21)
(3.22)
tj +
1 Л P (S )C(S )l 2
1 p(S)C(DS)lx2SS'V.,, = — Pj Cp,j xjS2TTj, j = 1,2,...nC — 1.
+ jXj+12 + 2pj)Cp;^SS' |j
(3.23)
1
2
2
Дрозин А.Д., Япарова Н.М., Математическая модель кинетики фазового перехода
Гольдштейн В.Я. при нагреве поверхности цилиндрического тела
2. Основные уравнения жидкой фазы:
а) (* + (R - *)(х}-12 -1)) -2(* + (Л - *)(ху -1))(Л - *)' х
х(К'х] +^'(1 - ху)))Ту- -
-(Х(-!/2 (* + (Л - *)(х;-/2 - 1)) + х("+!/2 (* + (Л - *)(х>+12 - 1)) +
Т +
р(^)С(Ь)( Л ?)2 '2
ру у' (5+(Л-«(ху -1))
У
^ (* + (Л - *) (х;+12 - 1)) + 2р^ (* + (Л - ху -1)) (Л Ч)' х
б) Т , =
(3.30)
(^С^) (Л -£)2 '2
х(Л'ху +^'(1 -ху )))ту+1 =-ру Ср,У ^ ^ (^ + (ЛЧ)(ху - 1))ТУ, у = пС + 1,...,2пс -1; (3.24)
б) Т2пс = ), t > 0. (3.25)
3. Основные уравнения поверхности раздела фаз:
а) Т (5 )= Т (^ = Т(5^); (3.26)
X5) Х( ^)
б) ?' =_пС_|т(5-^)-Т I__^_(т - • (3 27)
б) м^р^ ^'' пС ^ Дк(.-.)р(.) ( л'I пС+1 > (. )
в) $ = | + ^'т. (3.28)
4. Начальные условия:
а) | = ^о (^о < Л); (3.29)
ф( Ло ху ) ,у = 1,2,...пС;
ф(Л + (Л-Л0)(ху -1)), у = пС +1,пС + 1,...,2пС
5. Дополнительные соотношения:
а) р(/}=р(5)(Ту), р^)(Ту); (3.31)
б) р)= ср)(Ту), ср)= ср)(Ту); (3.32)
в) Х(у5) = Х(5} (Ту), = Х(1) (т у). (3.33)
3.4. Алгоритм расчета по конечно-разностной модели
3.4.1. Разработка алгоритма
Для решения уравнений (3.22)-(3.33) воспользуемся алгоритмом, предложенным в работе [5]. Применительно к нашей задаче, согласно этому алгоритму, переход от k -1 -го временного слоя к k -му осуществляется с помощью следующих последовательных операций. Задавшись каким-либо значением , решаем систему уравнений (3.22)-(3.26), а затем по уравнениям (3.27), (3.28) находим новое значение . Повторяем вычисления со всё новыми значениями до тех пор, пока новое значение практически не перестанет отличаться от старого.
Рассмотрим этот процесс подробнее. Условимся обозначать величину Т, относящуюся к временному слою k и итерации с номером р символом Т(р). Пусть для k -го слоя уже опреде-
е(т) )(т) I ■ п с\ ^(¿)(т) / . с ~ с\
лены величины, ^ ', Ту ' ', I у = 0,...,п I, Ту Л ', I у = п ,...,2п I, относящиеся к т -й итерации. Теперь по формулам (3.27), (3.28) найдем т+1) и т+1).
Уравнения (3.42), (3.43) можно представить в виде:
Т0(5 )=^т— 5 ) + Л; (3.34)
уз - СуТ5) + ВуТ(5+1=-Еу, (у = 1,..., пс -1), (3.35)
где Ау, Су-, В у, Fу - коэффициенты при Т^) в уравнениях (3.22), (3.23). С учетом номеров итераций, в соответствии с [5], запишем уравнения (3.34), (3.35) в виде:
Т(5)(т+1) = ^(т)т(5)(т+1) т); (3 36)
А(т)т(3( т+1) - су.т)т)5)(т+1) + в(т)т)+1)( т+1) = - у), (у = 1,..., пс -1). (3.37)
В этих уравнениях все коэффициенты вычисляются при уже известных значениях
(5)(т) (5)(т) л (5)(т) е,(т+1) е(т+1) -ту
ру- , сРу' , Xу , , ^ , и, следовательно, могут быть определены. Для решения этих
уравнений применим метод правой прогонки [7]. Последовательно определим прогоночные коэффициенты а у, р у, у = 1,..., пс по формулам:
«1 ау+1 = уУ(<у)-ауА(т)), у = пс -1;
&j=(4%+j)/(Cm)-jjm)), j=l,..., -1
(3.38)
двигаясь по сетке слева направо. Если значение Т(5)( определено, то все остальные значения
п
у)(*+1) можно найти по формуле [6]
т(5—(т+1) =а)(т+1) + Ру, у = пс, пс -1,..., 2,1, (3.39)
двигаясь по сетке справа налево.
Аналогично вышеизложенному представим уравнения (3.24)-(3.25) для перехода к т + 1-й итерации в виде
т+1)=Ц ^ ) (3.40)
А(ь)(т)т(ь)(т+1) - С(Ь)(т)Т(Ь)(т+1) + В(Ь)(тТ(ь)(т+1) = -F(L)(т) у = п +1 2пс -1 (3 41)
где Ау, Су, В у, Fj - коэффициенты при Тр) в уравнениях (3.24).
Решаем уравнения (3.24)-(3.25) методом левой прогонки [7]. Последовательно определяем прогоночные коэффициенты а у, р у, у (у = пс + 1,...,2пх -1) по формулам:
«2„< = 0, р2пс = ^), ау = А(.ту(с(.т)-ау+1В(.т));
' (3.42)
ру = (в(т)Ру+1 + Fjm))|(Cm)-ау+1Вут)), у = 2пс - пс - —
двигаясь по сетке справа налево. Если значение Т^)( т+1) определено, то все остальные значения
п
т+—) находятся по формулам [7]:
у)(т+1) =ау(т+1) +ру, у = пс +1,..., 2пс -1 (3.43)
при движении слева направо по узлам сетки.
После того, как найдены все значения температур на т -й итерации, по формулам (3.27),
(3.28) находим (т+2) и т+2). Если разность т+2) -£,(т+1) мала, переходят к расчету следующего слоя. В противном случае, продолжают итерации.
Дрозин А.Д., Япарова Н.М., Гольдштейн В.Я.
Математическая модель кинетики фазового перехода при нагреве поверхности цилиндрического тела
4. Пример расчета
В качестве примера приведены результаты расчета прогрева цилиндра из олова диаметром 5 мм, первоначально имеющего температуру 20 °С, температура поверхности которого с момента I = 0 стала поддерживаться равной 300 °С (рис. 4, 5).
Рис. 4. Изменение температуры в центре образца со временем
в) г)
Рис. 5. Распределение температуры по сечению образца через 0,01 с (а), 0,08 с (б), 0,36 с (в), 1 с (г)
после начала прогрева
Параметры расчета
Параметр Sn (тв.) Sn (ж.)
Удельная теплоемкость (Дж/(кг К) cp) = 243 cp) = 240
Удельная плотность (кг/м3) р(S) = 7280 р (L) = 6850
Коэф. теплопроводности (Вт/(м К)) тв. S) = 65,8 L) = 31,6
Удельная теплота плавления (°С) 231 ,93
Удельная теплота плавления (Дж/кг) 59 000
Подобная постановка задачи имеет отношение к термообработке деталей, к расплавлению кусков шихты. Она также интересна метрологам для обоснования теории самотестирующегося датчика [8-11].
Исходные данные [12-13] приведены в таблице.
Заключение
Разработана методика, позволяющая рассчитывать скорости процессов, происходящих при нагреве или охлаждении образца с учетом происходящих при этом фазовых переходах.
Литература
1. Любое, Б.Я. Теория кристаллизации в больших объемах /Б.Я. Любое. - М.: Наука, 1975. -256 с.
2. Скрипов, В.П. Спонтанная кристаллизация переохлажденных жидкостей / В.П. Скрипов, В .П. Коверда. - М.: Наука, 1984. - 232 с.
3. Авдонин, Н.А. Математическое описание процессов кристаллизации / Н.А. Авдонин. -Рига: Зинатне, 1980. - 180 с.
4. Тихонов, А.Н. Уравнения математической физики / А.Н. Тихонов, А.А. Самарский. - М. : Наука, 1972. - 736 с.
5. Будак, Б.М. Разностная схема с выпрямлением фронтов для решения многофронтовых задач типа Стефана / Б.М. Будак, Н.Л. Гольдман, А.Б. Успенский // Докл. АН СССР. - 1966. -Т. 167, № 4. - С. 735-738.
6. Самарский, А.А. Теория разностных схем /А.А. Самарский. - М. : Наука, 1977. - 656 с.
7. Самарский, А.А. Методы решения сеточных уравнений /А.А. Самарский, Е.С. Николаев. -М. : Наука, 1978. - 592 с.
8. Шестаков, А.Л. Модель самодиагностирующегося датчика параметра с нелинейной функцией преобразования /А.Л. Шестаков, А.С. Семенов //Измерение. Мониторинг. Управление. Контроль. - 2015.- № 1 (11). - С. 17-22.
9. Шестаков, А.Л. Методы теории автоматического управления в динамических измерениях / А.Л. Шестаков. - Челябинск: ЮУрГУ, 2013.
10. Самокалибрующийся термометр на основе точек плавления. Конструкция и алгоритмы работы / М.Д.Белоусов, В.В. Дьячук, Д.А. Мирзаев, А.Л. Шестаков // Вестник ЮУрГУ. Серия «Компьютерные технологии, управление, радиоэлектроника». - 2013. - Т. 13, № 1. - С. 26-33.
11. Белоусов, М.Д. Оценка собственного состояния термометров сопротивлений / М.Д. Белоусов, А.Л. Шестаков, Н.М. Япарова // Вестник ЮУрГУ. Серия «Компьютерные технологии, управление, радиоэлектроника». - 2012. - № 35. - С. 105-109.
12. Физические величины. Справочник / под ред. И.С. Григорьева, Е.З. Мейлихова. - М: Энер-гоатомиздат, 1991. - 1250 с.
13. Рабинович, В.А. Краткий химический справочник / В.А. Рабинович, З.Я. Хавин. - Л: Химия, 1978. - 392 с.
Дрозин А.Д., Япарова Н.М., Гольдштейн В.Я.
Математическая модель кинетики фазового перехода при нагреве поверхности цилиндрического тела
Дрозин Александр Дмитриевич, д-р техн. наук, профессор кафедры материаловедения и физико-химии материалов, Южно-Уральский государственный университет, г. Челябинск; drozi-nad@susu.ru.
Япарова Наталья Михайловна, канд. физ.-мат. наук, доцент кафедры прикладной математики и программирования, Южно-Уральский государственный университет, г. Челябинск; уара-rovanm@susu.ru.
Гольдштейн Владимир Яковлевич, д-р техн. наук, профессор, ведущий научный сотрудник научно-внедренческого центра, ЗАО «Ферросплав», г. Челябинск.
Поступила в редакцию 7 октября 2016 г.
DOI: 10.14529/met160406
MATHEMATICAL MODEL OF THE KINETICS OF PHASE TRANSITIONS AT THE HEATING OF THE SURFACE OF A CYLINDRICAL OBJECT
A.D. Drozin1, drozinad@susu.ru, N.M. Yaparova1, yaparovanm@susu.ru, V.Ya. Gol'dshteyn2
1 South Ural State University, Chelyabinsk, Russian Federation,
2 JSC "Ferrosplav", Chelyabinsk, Russian Federation
If due to heating or cooling of a single-component object, a phase transition occurs, the propagation of heat significantly changes. When passing through the liquidus temperature, enthalpy, density and thermal conductivity as a function of temperature undergo discontinuity. By heating the surface of the solid metal cylindrical object above its melting temperature, the liquid phase occurs at the surface, and the corresponding phase boundary moves to the cylinder axis. The article presents the derivation of heat propagation equations in the body, taking into account, in addition to the phase transition, the change in the size of parts of the body due to density dependence on temperature. To the resulting system of equations, straightening fronts method was applied: such coordinates have been entered that for which interface surface is stationary. The resulting system of differential equations is reduced to finite-difference equations. A computer program was developed to solve the resulting system of difference equations. Results of one of such calculations are given in the article. The developed method allows to calculate the velocity of the phase boundary, as well as the object temperature at any point and at any time. The results may be interesting not only for metallurgists, but also in metrology to develop self-testing temperature sensor theory.
Keywords: heat conduction; phase transition; model.
References
1. Lyubov B.Ya. Teoriya kristallizatsii v bol'shikh ob"emakh [Theory of Crystallization in Large Volumes]. Moscow, Nauka Publ., 1975. 256 p.
2. Skripov V.P., Koverda V.P. Spontannaya kristallizatsiya pereokhlazhdennykh zhidkostey [Spontaneous Crystallization of Undercooled Liquids]. Moscow, Nauka Publ., 1984. 232 p.
3. Avdonin N.A. Matematicheskoe opisanie protsessov kristallizatsii [Mathematical Description of Crystallization Processes]. Riga, Zinatne Publ., 1980. 180 p.
4. Tikhonov A.N., Samarskiy A.A. Uravneniya matematicheskoy fiziki [Equations of Mathematical Physics]. Moscow, Nauka Publ., 1972. 736 p.
5. Budak B.M., Gol'dman N.L., Uspenskiy A.B. [Difference Scheme with Straightening of Fronts for Solving Stefan Type Multi-Front Problems]. Doklady ANSSSR, 1966, vol. 167, no. 4, pp. 735-738. (in Russ.)
6. Samarskiy A.A. Teoriya raznostnykh skhem [Theory of Difference Schemes]. Moscow, Nauka Publ., 1977. 656 p.
7. Samarskiy A.A., Nikolaev E.S. Metody resheniya setochnykh uravneniy [Methods of Solving Grid Equations]. Moscow, Nauka Publ., 1978. -592 s.
8. Shestakov A.L., Semenov A.S. [Model of a Self-Diagnosing Parameter Sensor with Nonlinear Transformation Function]. Izmerenie. Monitoring. Upravlenie. Kontrol', 2015, no. 1 (11), pp. 17-22. (in Russ.)
9. Shestakov A.L Metody teorii avtomaticheskogo upravleniya v dinamicheskikh izmereniyakh [Methods of Automatic Control Theory in Dynamic Measurements]. Chelyabinsk, South Ural St. Univ. Publ., 2013.
10. Belousov M.D., D'yachuk V.V., Mirzaev D.A., Shestakov A.L. [Self-Calibrating Thermometer Based on Melting Points. Construction and Operation Algorythms]. Bulletin of the South Ural State University. Ser. Computer Technologies, Automatic Control & Radioelectronics, 2013, vol. 13, no. 1, pp. 26-33. (in Russ.)
11. Belousov M.D., Shestakov A.L., Yaparova N.M. [Estimation of Eigenstate of Resistive Thermometers]. Bulletin of the South Ural State University. Ser. Computer Technologies, Automatic Control & Radioelectronics, 2012, no. 35 (294), issue 17, pp. 105-109. (in Russ.)
12. Fizicheskie velichiny. Spravochnik [Physical Values. Reference Book]. Grigor'ev I.S. and Mey-likhov E.Z. (Eds). Moscow, Energoatomizdat Publ., 1991. 1250 p.
13. Rabinovich V.A., Khavin Z.Ya. Kratkiy khimicheskiy spravochnik [Brief Chemical Reference Book]. Leningrad, Khimiya Publ., 1978. 392 p.
Received 7 October 2016
ОБРАЗЕЦ ЦИТИРОВАНИЯ
FOR CITATION
Дрозин, А.Д. Математическая модель кинетики фазового перехода при нагреве поверхности цилиндрического тела / А.Д. Дрозин, Н.М. Япарова, В.Я. Гольдштейн // Вестник ЮУрГУ. Серия «Металлургия». - 2016. - Т. 16, № 4. - С. 54-66. DOI: 10.14529/те060406
Drozin A.D., Yaparova N.M., Gol'dshteyn V.Ya. Mathematical Model of the Kinetics of Phase Transitions at the Heating of the Surface of a Cylindrical Object. Bulletin of the South Ural State University. Ser. Metallurgy, 2016, vol. 16, no. 4, pp. 54-66. (in Russ.) DOI: 10.14529/met160406
