CC BY
12УДК 539.3, 621.74
Г. М. Севастьянов, В. И. Одиноков, И. Г. Сапченко
ОБ ОДНОМ АЛГОРИТМЕ РЕШЕНИЯ УРАВНЕНИЯ ТЕПЛОПРОВОДНОСТИ ДЛЯ КЕРАМИЧЕСКОЙ ОБОЛОЧКОВОЙ ФОРМЫ В ЛИТЬЕ ПО ВЫПЛАВЛЯЕМЫМ МОДЕЛЯМ
Институт машиноведения и металлургии ДВО РАН
Приводится оригинальный алгоритм решения задачи распространения теплоты в керамической форме при заливке и затвердевании расплава стали и расчёта фронта кристаллизации. Алгоритм предназначен для использования в программных продуктах научно-исследовательской и прикладной направленности для литейного производства.
Ключевые слова: литейное производство, керамические оболочковые формы, уравнение межфазового перехода, уравнение теплопроводности, численные методы.
Введение. Рассмотрим процесс заливки стали температуры в керамическую оболочковую форму, разогретую до температуры 90 . Ограничимся рассмотрением стандартного стояка формы, представляющего собой осесимметричное тело, состоящее из сферической и цилиндрической частей. Геометрические параметры стояка: Я - радиус сферической части, £ - высота цилиндрической части, Н - толщина формы. Ввиду осевой симметрии и однородности материала формы и заливаемого металла будем рассматривать сегмент стояка с произвольным углом развертки ю. Зададим систему криволинейных координат и разобьём рассматриваемую область ортогональными плоскостями, как показано на рис. 1. Число элементов разбиения по координате а: в форме - , в металле - Ы1М ; по координате а2: в сферической части стояка - N°, в цилиндрической части стояка - Ж^1. Таким образом, общее число элементов равно + Ы1м + ЫН ).
Рис. 1. Область расчёта с обозначением поверхностей, разбитая на ортогональные элементы
© Севастьянов Г. М., Одиноков В. И., Сапченко И. Г., 2010.
Формализация процесса. Моделируемый процесс распространения теплоты состоит из двух этапов: этап заливки расплава, на котором непрерывно меняется граница раздела металл - форма и этап последующего застывания металла в форме.
Время заполнения формы рассчитаем при постоянной скорости поступления расплава:
Т1 = Т ф+ Т Я:
где т - время заполнения сферической части; тн - время заполнения цилиндрической части;
_ 2лЯ3 _ пЯ2 И
тФ ~ , , ти - ,
Эу V
V - скорость поступления расплава.
Общее время, в течение которого будем рассматривать процесс, складывается из т и времени т2 - продолжительности второго этапа.
Примем ряд предположений:
• в процессе заливки не происходит фазовый переход расплава; это предположение в реальных условиях выполнено практически в точности из-за незначительности времени заполнения формы и процессов перемешивания металла под напором струи; по тем же причинам принимается, что температура расплава остаётся неизменной до начала второго этапа процесса;
• в течение некоторого времени после окончания заливки в жидком металле сохраняется область, температура которой равна начальной температуре заливаемого металла; граничная поверхность этой области на рис. 1 обозначена ; фактически, представленное предположение адекватно, если поверхность задана на достаточном удалении от внутренней поверхности формы;
• температура внешней грани формы во время заливки постоянна, а в дальнейшем изменяется по заданному закону 9| - 0(т), где т - время с начала процесса;
• считаются заданными температурные зависимости теплофизических величин для стали и керамики формы: теплопроводность Хм (9) и Хр (9), теплоёмкость См(9) и СД9), плотность Р м (9) и Р F (9).
Зададим разбиение времени процесса по шагам: ш1 - число шагов при заполнении сферической части стояка; ш1н - число шагов при заполнении цилиндрической части стояка, ш11 - число шагов на втором этапе процесса (этап затвердевания полностью залитого металла).
Тогда начальные условия на первом шаге процесса включают в себя поле температур в элементах формы и расплава 9^, где /,] - индексы, определяющие положение элемента по
сетке; I - номер столбца (нумерация от оси симметрии); ] - номер элемента в столбце (нумерация от внешней грани формы):
9 [9М, при ^ < ] < ^ + КМ [9F, при ] < N1 .
На последующих шагах начальными условиями по температуре являются рассчитанные на предыдущем шаге поля температур.
Граничные условия на этапе 1 9^ - 9^ , 9|5„ - 90, 9|5 - 90, или в матричной форме:
_/9М. пРи ' < п(Т)
190, при I > п(т)'
а = )°м > пРи 1 ^ n(T) СП
9гД - 9F, / -1,(Кф1 + К°). (2)
Здесь п(т) - число дуг на внутренней поверхности формы, составляющих &2, иными словами, это номер последнего элемента «залитого» металлом. Оно определяется следую-
щим образом. Возьмём произвольное т < т . Объём расплава, поступивший в форму за это время, равен ту . Этот расплав занимает в стояке элемент объёма высотой к, равный
пj^R2 " (R - x)2 dx (рис. 2), тогда tv =1 nh2(3R - h).
3
Это кубическое уравнение решается после соответствующих преобразований относительно к по формулам Кар-дано или же методом половинного деления на полуинтервале к е [О,Я) , где оно имеет один действительный корень. Да-
1 Н
лее, учитывая соэд = 1 — —, выразим
п(т) = round
Рис. 2. Элемент объёма в сферической части формы
2Nj • arccosil - h(T)
л
при
т < т„
где
round (x) - округленное к ближайшему целому значение х. При произвольном тф < т < тф + тн n(т) = Ny + round
V
Граничные условия на этапе 2 б|5 = , 9|5 = 9(т) или
^(т - Тф )NH Л nR 2 H
9,(n1+NL) = 9M , 9а = 9(т), i = U< + NH)
(3)
Алгоритм расчёта температурных полей. Согласно методу, описанному в [1, 2], предпримем следующую последовательность шагов.
1. Производится подсчёт дуг ортогональных элементов:
• массив дуг £21}., отложенных по координате а, ^ = 1, (N1 + N° +1), ] = 1, (+ ),
S2,., =
H
Nf d
—, если j < NlF
N,
—, если j > NF
где d - расстояние между поверхностями S3 и S2;
массив дуг S1i,j, отложенных по координате а2, i = 1, (NF + N"), j = 1, (NlF + NlM +1),
(R + H - K1)-^, еслиi < N
S1,7 = 1
2 N
r11
' ф
0, если j = 1
где K1 = .
^ S2hk, если j *1
к=1
• массив дуг S3i,j,
j = 1,(NF + NM +1)
S
N
п,еслиi>N
H
отложенных по координате а3, i = 1,(N" + N" +1),
0
У
S3i , = (R + H - Ki)sin
f \ -Лт K 2
V2< J
cosro.
10,если 1 = 1г., . , АгП , ,
где К1 = Ь-'х 7 К2 V-1, если' ^ 1
где К1 = { ^^2^, если ] ф 1' К 2 = {жфп, если/ > ^ +1'
=1
2. Задаётся начальное поле температур (на первом шаге используются начальные условия задачи, на последующих - значения, рассчитанные на предыдущем шаге) 0;.,
i = 1,(N" + N"), j = 1,(N +1). Согласно граничным условиям (1)-(2) задаются элементы
массива 0и , / = 1,(N° + Ж"), 1 = 1, (N1 +1) (внешняя поверхность формы и поверхность соприкосновения металла с формой).
3. В соответствии со сформированным полем температур 00. для каждого элемента области задаются массивы теплопроводности Х^, теплоёмкости Сг ; и удельного веса уг- ■.
4. Для внутренних элементов области производится расчёт коэффициентов Шг;, ¿12г;,
t22u, 121 i,j i = 2, (N" + NH -1), j = 2, NF по формулам:
tii = + K-1 )Skj • (S3,j + S3,+i.j )AA (41)
4 CjjYi,jVi,j(S2y + S2,+i,7 + S2.-i + S2,+и-i), (.)
tl2 = (hj + hj+i)Sii.j+i(S 3i.j+i + S3i+i,.j-+i) • A4 (42)
h] CjYjVj (S2i.j + S2.+i.j + S2i.j+i + S2.+и+i), (.)
t22i = (4j + ^+ij)S2i+i.j (S3i+i.j + S3i+i.j+i)AA ^ (4 3)
4 Ci.jYi.jVi.j (Sii.j + Sii.j+i + Sii+i.j + Sii+i.j+i) , .
t21 _ + Кц )S 2u (S3,j + S3y+i)AA
iJ CuYi.jVi.j(Sii.j + Sii.j+i + Sii-i.j + Sii-i.j+i), (.)
где Атк - величина к -го шага по времени; Vj - объём элемента (i,j), вычисленный приближённо по формуле Vl} = (S\Uj + Sll}+1)(S2U + S2г+1у.)(S3hJ + S3U+1 + S3,+w + S3i+1J+i)/16 .
Отметим, что согласно формулам (4.1)-(4.4) указанные коэффициенты не зависят от угла развертки ю сегмента формы. Это подтверждает возможность произвольного его выбора. По формуле
A««- e0j +9JMj +9Л1^ + е:^1 • t22tJ + eje^t21,. j
9, J --,
4 1 + t11tJ + t12tJ +122tJ +121tJ
где 0^ - элементы сформированного в пункте 3 поля температур; iter - номер итерации,
представляющей собой итеративную процедуру, рассчитываются значения температуры во внутренних элементах области.
5. Для элементов первого столбца элементов от оси симметрии (i — 1, j — 2, NlF ) коэффициент 121. j — 0, коэффициенты t12,j и 122,j рассчитываются по формулам (4.2) и (4.3), а коэффициент t11 , исходя из соображений симметрии, по формуле
2 (S3:,j + S 3,+1,J )Д4
tu, =
i.j
Ci.j Yi.jVi.j(S 2i.j + S 2i+i.j + S 2i.j -i + S 2i+u-i)
Итеративная процедура для этих ячеек имеет вид
90 j + 9it;-11t11i,j + 9it;+11t12i, j + 9lt+j 22,
9,- i =-.
l,j 1 + t11, + t12,,j +122,, j
6. Температура в элементах последнего столбца (i = Nipp + N°, j = 2, NF ) рассчитывает-
fl0 (Yter
niter 9i,j + 9i-1,j
ся, исходя из сглаживающего условия 9^- = ——-.
Пункты 4-6 повторяются до сходимости по температуре, то есть до тех пор, пока не выполнится условие max{|9ljr-9l,jr-1 |j<s, i = 1,(N" + Npp), j = 1,(NlF +1), где s - требуемая
точность. Сходимость процедуры доказана в [2].
Описанная процедура повторяется (mp + mlH) раз, что имитирует заполнение формы
целиком. На этом завершается первый этап расчёта. Отметим, что п. 1 выполняется только один раз в начале процедуры, поскольку считается, что термические деформации формы, в силу хрупкости керамики, остаются малы в рассматриваемом процессе и не вносят сильных искажений в сетку (по крайней мере, если не произошло нарушения целостности формы).
На втором этапе процесса начинается процесс кристаллизации металла в форме. Расчёт температурных полей на этом этапе также происходит по описанной процедуре 1 -6, за единственным отличием: область расчёта увеличивается за счёт элементов в жидком и закристаллизовавшемся металле, соответственно все массивы в пп. 2 и 3 (9г-., ., Сг ;, у; .)
имеют размерность i = 1,(N° + N"), j = 1, (Nj, + NlM ) . Соответствующим образом изменяются и границы массивов в пп. 4-6. Начальным полем температур в п. 2 на втором этапе служит поле
9,,=19У
9т+тн, при j < Nlp +1
9M, при j > Nlp +1
где I = 1,« + ] = 1, (N + К1М).
В п. 2 используются граничные условия (3). Число шагов на втором этапе процесса равно ш11. Для каждого шага по времени после расчёта поля температур выполняется расчёт фронта кристаллизации.
Алгоритм расчёта фронта кристаллизации. Расчёт фронта кристаллизации проведём следующим образом [3]. Рассмотрим уравнение межфазового перехода
X, —X 2 Хр, (5)
дп дп дт
где 9 и 92 - температура в твердой и жидкой фазах соответственно; X и Х2 - коэффициенты теплопроводности в соответствующих фазах; А - толщина корочки; Ь - скрытая теплота плавления; р - плотность металла; п - нормаль к границе раздела двух фаз, направление нормали совпадает с направлением координатной оси 1 .
Решение этого уравнения получается при некоторых упрощающих предположениях.
Во-первых, примем градиент температуры в жидкой фазе равным нулю, то есть
092= О. (6)
дп
Во-вторых, примем, что температура в твердой фазе изменяется линейно, то есть
91 = 91 + (9? — 91)^, (7)
где 0| - температура на границе образовавшейся на временном шаге Ах корочки и твердого металла (или формы на первом временном шаге кристаллизации); 02 - температура на границе образовавшейся на временном шаге Ат корочки и жидкого металла; х - координата по толщине корочки. Тогда, учитывая (6) и (7), уравнение межфазового перехода запишется в виде
(8)
А ат
Сгруппируем и проинтегрируем:
^02 — 01 е
01 01VI =[ АаА, (9)
Ь • р
откуда зависимость для определения прироста по толщине закристаллизовавшейся корочки на временном шаге Ат будет
А =
"V
2СА6Л Ат
pL
где A0j - перепад температур в твердой фазе вблизи фронта кристаллизации.
Тогда используем следующую итерационную процедуру для расчёта фронта кристаллизации металла на к -м шаге второго этапа:
1) для каждого номера столбца элементов i определим такой номер элемента в столбце
jк, что 0 0 < 0к, 0 .„ > 0к, где 0к - температура кристаллизации металла;
i,jk i,j0+1
2) итеративно (iter - номер итерации) рассчитывается прирост толщины корочки на к-м временном шаге по формуле
АГ =
2 • j " ^W ■
-Ат^ , где j(k_r - наибольший номер элемента в твёрдой
р(ек)Ь
фазе i -го столбца на (k -1) -м временном шаге, j = N^ ; 3) рассчитывается величина
/Г = round
( к iter \
1
Т NM
V d
где d - расстояние между поверхностями и .
Пункты 2-3 повторяются до тех пор, пока не выполнится равенство j1£er = jer~l, тогда jk = jf. Таким образом, найдено уточнённое значение величины прироста корочки затвердевшего металла для i -го столбца элементов на к-м временном шаге.
После этого уточняются температуры в i -м столбце с учётом фазового перехода:
0i j = 0к при j = jk, j'0 . Описанную процедуру необходимо провести для каждого столбца элементов области. Алгоритм отражает картину протекания процесса кристаллизации металла и позволяет определить границу раздела фаз в данный момент времени.
Совокупность описанных алгоритмов была реализована в среде Compaq Visual FORTRAN 6.0 и используется для численного моделирования напряженно-деформированного состояния керамических оболочковых форм.
Библиографический список
1. Математическое моделирование сложных технологических процессов / В.И. Одиноков [и др]. - М.: Наука, 2008. - 176 с.
2. Одиноков, В. И. Численное исследование процесса деформации материалов бескоординатным методом / В.И. Одиноков. - Владивосток: Дальнаука, 1995. - 168 с.
3. Севастьянов, Г. М. Моделирование процесса заливки металла в керамическую осесиммет-ричную оболочковую форму / Г. М. Севастьянов, В. И. Одиноков, И. Г. Сапченко // Прикладные задачи механики деформируемого твёрдого тела и прогрессивные технологии в машиностроении: сб. ст. - Комсомольск-на-Амуре: ИМиМ ДВО РАН, 2009. Вып. 3. Ч. 2. С. 18-38.
Дата поступления в редакцию 06.04.2010
G. M. Sevastyanov, V. I. Odinokov, I. G. Sapchenko
THE SOLUTION ALGORITHM FOR HEAT CONDUCTION EQUATION FOR CERAMIC SHELL MOLD IN INVESTMENT CASTING
Original solution algorithm for heat conduction problem in ceramic mould during a process of pouring and solidification of a steel melt and crystallization front computation are given. Algorithm is assigned for use in research and applied software products for foundry engineering.
Key words: foundry engineering, ceramic shell moulds, interphase conversion equation, heat conduction equation, numerical methods.
