Тестирование численных моделей затвердевания металла Текст научной статьи по специальности «Физика»

УДК 519.63

ТЕСТИРОВАНИЕ ЧИСЛЕННЫХ МОДЕЛЕЙ ЗАТВЕРДЕВАНИЯ МЕТАЛЛА З.К. Кабаков, И.Ю. Мазина

В статье изложен способ тестирования численных моделей задачи Стефана с использованием точного решения. Представлены результаты исследования влияния настроечных параметров численного алгоритма на погрешность моделирования

Ключевые слова: математическая модель затвердевания металла, тестирование, задача Стефана, погрешность решения

Для изучения закономерностей затвердевания и охлаждения металлических слитков применяют математические модели. Основу этих моделей составляет численное решение задачи Стефана, включающее определение границы фазового перехода из жидкого состояния в твердое и полей температуры в жидкой и твердой фазе слитка. Для примера приведем математическую модель одномерного симметричного процесса

затвердевания, которая включает сквозное уравнение теплопроводности, общее для твердой и жидкой фаз:

где тл = ткр + АТ > Тс = Ткр - ^ - фиктивные

температуры начала и окончания затвердевания элементарного объема, с(Т) - теплоемкость

металла, Ь - удельная теплота кристаллизации.

При моделировании слитков из двойных сплавов величины Тл, Тс и АТ имеют физический смысл температур ликвидус и солидус и интервала температур кристаллизации.

На рис. 1 показана схема затвердевания

половины сечения плоского слитка.

дТ=51 ъТ

эфР д/ дх ^ дх, интегрируемого в области: 0 < х < £, 0 < / < /к начальное условие: Т| = Т0;

граничные условия:

при х = 0

при х = S

-A — = a(T - T );

дх ' cp)

A — = 0, дх

(1)

(2)

(3)

(4)

где р - плотность металла, а - коэффициент теплоотдачи, £ - половина толщины слитка, Тср -температура среды, Т - начальная температура расплава, X - коэффициент теплопроводности.

При моделировании затвердевания чистого металла источник тепла кристаллизации, действующий на фронте кристаллизации (х = в, Т = Ткр), «размазывают» по фиктивному

интервалу температур АТ = Тл - Тс. При этом выделение тепла кристаллизации в интервале АТ учитывают с помощью эффективной теплоемкости

L"эф

задаваемого выражением:

Сэф =

c(T), c(T) +

T<T, T>T;

с? л’

T < T < T

Кабаков Зотей Константинович - ЧГУ, д-р техн. наук, профессор, тел. 89095998386 .

Мазина Ирина Юрьевна - Череповецкий филиал ВКА им. А.Ф. Можайского, соискатель, e-mail: [email protected].

Рис. 1. Схема расчетной области: 1- твердая фаза, 2 -жидкая фаза, 3- двухфазная зона, А1 - ширина двухфазной зоны, вс, в и вл - координаты изотерм солидуса, температуры кристаллизации и ликвидуса

Система уравнений (1)-(4) в общем случае может быть решена только численным методом. При использовании метода конечных разностей (МКР) значения температуры находят в узлах расчетной области, координаты которых находят по формуле х1 = (I - 0,5)Ах , для дискретных моментов

времени /п = А/ • п , где I = 0, N +1. N - количество узлов внутри расчетной области, 0 и N+1 - номера фиктивных узлов, находящихся за пределами

_ Дх

области на расстоянии —,

Д S

Дх = n ~ расстояние

n = 0,

между узлами.

(n = 0 - начальный

д _

момент

моменты времени времени), Дґ -

расчетный шаг по времени. Для краткости температуру Т (х1, гп) обозначают Т " .

При использовании явной схемы аппроксимации производных по координате температуру в следующий момент времени п +1 в N внутренних узлах определяют по формуле:

Tn+i t

Дґ

‘ С Tn )р{т,п )Дх2

A. 1 T+1 - T )-A 1 T - T-1)

где

i _ 1, N, A +1 _я((т+1 + Tn )), * 2

A 1 _я((т-1 + Tin )).

Температуру в начальный момент времени задают по формуле:

T0 _ T0 для i _ 0, N +1;

температуру в фиктивных узлах определяют по формулам

T _ 10

(1 - ж )т + 2ж Tcp

1+ж

ж_

аДх 2A ;

T _ T 2n+1 2n-

Толщину твердой корки определяют по температуре кристаллизации в поле температуры в цикле по I = 2..^ при условии:

если Т 1 < Ткр < Т , то в = Ах | г - — | +Ах • —-— .

'-1 41 ^ 2 ^ Т1 - Т-1

Численное решение при явной схеме аппроксимации является условно устойчивым. В этом случае алгоритм МКР устойчивый при

выполнении условия Дґ _

Дх2

ky ■ a'

где ky > 2 .

Погрешность численного решения в данном случае будет зависеть от настроечных параметров алгоритма N и ку, а также параметра AT в

математическом описании затвердевания пластины. Необходимо эти параметры выбрать таким образом, чтобы погрешность результатов моделирования не превосходила заданную.

Как правило, в известных работах, в которых используется математическое моделирование процесса затвердевания, этому вопросу не уделяется внимания.

В данной работе предложен способ тестирования численных методов решения задачи Стефана с использованием точного решения этой задачи [1], который позволяет устранить этот пробел. В работе [1] приведено точное решение задачи Стефана при граничном условии I рода, включающее поля температуры: в твердой фазе

( \

erf

Чх,ґ) _тп +T -тп)•

(

(5)

erf

в жидкой фазе

erfc

T2(х,ґ) _ т -( -Тр)•

2yfa2t

(6)

erfc

2yfa2

(7)

и формулу для расчета координаты границы фазового перехода

є = ,

где к - корень трансцендентного уравнения (8)

— Г

4а =р!к^П • (8)

А Т - Тп)

JO! • erf

А (Т - Тр )

2JaO1

Ja2 • ЄФ

2у[°2

В формулах (5)-(8) использованы следующие обозначения: А1 и Л2 - теплопроводность твердой и жидкой фазы, с1 и с2 - теплоемкость твердой и

жидкой фазы, р - плотность, а. =—— -

СР

температуропроводность (г = 1 - твердая фаза, г = 2 - жидкая фаза), Тп - температура поверхности,

T

J- п

начальная температура расплава,

erf (х) _~2= ] e

d#;

\[я 0

erfc(х) _ 1 - erf (х) _

j e # d# .

Тестирование выполним на конкретном примере затвердевания и охлаждения алюминиевой пластины толщиной 2£=0,03 м, охлаждаемой в симметричных условиях. Исходные данные для моделирования и расчета по формулам (5)-(8) приведены в таблице.

Размер расчетной области и граничные условия для модели выберем с учетом того, что точное решение получено для граничных условий первого рода и расчетной области в виде полупространства, а задача затвердевания сформулирована для граничных условий третьего рода и расчетной области конечной толщины.

Тестирование результатов моделирования выполнили путем сравнения с точным решением и оценкой погрешности моделирования. При тестировании исследовали влияние настроечных параметров конечно-разностного решения задачи Стефана на результаты и погрешность моделирования.

Результаты исследования представлены на рис.2-7.

На рис. 2 приведены результаты по росту корки с течением времени для различных значений количества узлов.

2

Исходные данные

№ пп Величина, размерность Значение

в модели в точном решении

1 Половина толщины, м 0,015 <Х)

2 Начальная температура, °С 710 710

3 Температура поверхности, °С - 600

4 Температура среды, °С 600 -

5 кг Плотность жидкой и твердой фаз, — м 2700 2700

6 Вт Коэффициент теплопроводности жидкой и твердой фаз, °С 238 238

7 Теплоемкость жидкой и твердой фаз, Дж кг-°С 920 920

8 Температура кристаллизации, °С 660 660

9 У Дж Удельная теплота кристаллизации, кг 357000 357000

10 Вт Коэффициент теплоотдачи, — м2 -°С 1010 -

11 Конечное время процесса, с 0,8 0,8

2 1.8 1.6 5 1,4

3

1.2

О.

0

і; 1

лз

1

I О.й £ о,б

0.4 0.2 О

О 1 2 3 4 5 6 7 8 9 10 11 12 13 14 16

Время, с10!

Рис. 2. Изменение толщины корки с течением времени:----------------------толщина корки, рассчитанная по формулам (7)-(8);

толщина корки, рассчитанная с помощью модели при к =5, АТ =5: я - N=30,------------ N=50, * - N=100,

- N=150.

Как следует из рис. 2, при N=100 и 150 На рис. 3 приведена зависимость

величина корки практически сходится к точному относительной погрешности моделирования роста решению. корки от времени для N=150, ку =5, АТ =5.

Рис. 3. Относительная погрешность расчета толщины корки в зависимости от времени (N=150, ку =5, АТ =5).

Относительная погрешность рассчитывалась для каждого момента п по формуле:

¿ = 1^—^-100% ,

А/

задать А/ = —.

р к

где к >1 - параметр для повышения

точности результатов моделирования в начале

где еп - результат моделирования роста корки в

*

момент времени п, еп - точное решение в момент времени п.

Как видно на рис. 3, в первые моменты времени / < 0,075 с относительная погрешность превышает 1%. Такое поведение погрешности в первые моменты времени объясняется высокой нестационарностью процесса при резком снижении температуры поверхности от 710°С до 600°С на первом шаге по времени. Для уменьшения погрешности в первые моменты времени необходимо измельчить расчетный шаг по времени А/ . Для этого следует в начальный момент времени

физического процесса. Затем А/р увеличивать по А/

формуле А/р = А/р +---до значения А/. В данной

к

работе при к = 20 достигнуто понижение погрешности расчета корки на начальном этапе приблизительно на 1%.

На рис.3 видно, что с уменьшением

нестационарности процесса (0,1 < / < 0,4) погрешность снижается почти до 0,6%. При / > 0,4 с погрешность начинает возрастать. Объяснение этого факта можно увидеть на рис. 4, на котором приведено поле температуры по толщине слитка в различные моменты времени, полученное моделированием и с использованием точного решения.

Рис. 4. Распределение температуры по толщине слитка для различных моментов времени:------точное решение

(формулы (5)-(8) ),---модель при N=50, ку =2.1, АТ =5.

На рис. 4 видно, что погрешность прогноза толщины корки с использованием модели связана с отклонением поля температуры в жидкой и твердой фазе от точного решения в первые моменты времени. Начиная с момента времени / =0,4 с опять начинается отличие численного решения от точного. Именно с этого момента начинает возрастать относительная погрешность прогноза величины корки. Отличие температурных полей начинается с жидкой фазы и связано с конечной толщиной пластины. Отсюда следует, что тестирование численного решения необходимо выполнять только до / <0,4 с или до толщины корки в = 0,2£ . Таким образом, чтобы протестировать затвердевание половины толщины пластины, необходимо в модели этот размер увеличить в 5 раз.

В данной работе также проведено исследование влияния количества узлов N и параметра ку на среднеквадратичную погрешность моделирования роста корки. Величину среднеквадратичной погрешности определяли по формуле:

8 =

1 1

-К -в) -100%,

где п =

к_

А/

на рис. 5 и 6.

. Результаты исследования приведены

Рис. 5. Зависимость среднеквадратичной погрешности 8 от N (ку =5, АТ =5).

Рис. 6. Зависимость среднеквадратичной погрешности 8 от ку (N=150, АТ =5).

Как видно на рис. 5 и 6 с увеличением N и ку погрешность уменьшается, что связано с

измельчением расчетных шагов по / и х в первом случае (рис. 5) и только по / во втором (рис. 6).

Исследование влияния фиктивного интервала АТ (рис. 7) показало, что имеется оптимальное значение АТ =5°С, при котором численное решение по динамике роста корки имеет наименьшее значение средней относительной погрешности. Уменьшение АТ приводит к тому, что количество узлов, попадающих в фиктивный интервал (А/ на рис. 1), становится недостаточным для

аппроксимации источника тепла кристаллизации. При увеличении АТ >5°С погрешность повышается, так как увеличивается различие характера температурного поля в окрестности координаты фазового перехода для численного и точного решений (рис. 4).

Рис. 7. Зависимость среднеквадратичной погрешности 8 от АТ (параметры модели: N=200, ку =20).

Следует отметить, что такое большое число узлов, обеспечивающих погрешность менее 1%, обусловлено «жестким» граничным условием I рода. При использовании граничных условий III рода следует ожидать существенного уменьшения допустимого количества узлов.

В заключение следует отметить:

1. Одновременное измельчение сетки и соответствующее уменьшение шага по времени согласно условию устойчивости влияет на уменьшение погрешности более эффективно, чем только уменьшение расчетного шага по времени при увеличении параметра ку .

2. В результате тестирования установлено, что для уменьшения средней относительной погрешности до 1% необходимо взять количество узлов сетки не менее 150 узлов. Следует отметить, что для менее «жестких» граничных условий, например III рода, это ограничение может измениться до существенно меньших N.

3. Уменьшению погрешности при «жестком» граничном условии (I рода) способствует введение в начале физического процесса меньшего, чем рассчитанного из условия устойчивости, шага по времени, который затем постепенно увеличивают до А/.

4. При тестировании численного решения задачи полного затвердевания пластины толщиной 28 при симметричном охлаждении размер ее

п=1

следует взять равным 10-8 и сравнение проводить на толщине 8.

5. Величину фиктивного интервала АТ при численном решении задачи Стефана следует подбирать из условия минимума

среднеквадратичной погрешности.

Литература

1. Любов Б.Я. Теория кристаллизации в больших объемах. - М.: Наука, 1975. 256 с.

Череповецкий государственный университет Череповецкий филиал Военно-космической академии им. А.Ф. Можайского

TESTING NUMERICAL MODELS OF METAL SOLIDIFICATION Z.K. Kabakov, I.Yu. Mazina

In article the method for testing numerical models of the Stefan problem using the exact solution is stated. Results of research of influence of setting parameters of numerical algorithm on a modeling error are presented

Key words: mathematical model of metal solidification, testing, Stefan problem, solution error