Решив эту систему, получили уравнение зависимости С2(t)=3,85 – 0,0025t.
После соответствующих преобразований и подстановки в соотношение (7) определим темп выделения твёрдой фазы путем дифференцирования уравнения (7) для последующего использования этого выражения в математической модели (6):
, тогдаВведём обозначение dm/dt = A, тогда из (6) получим следующее дифференциальное уравнение:
a (t –tep )F dτ = –VСp dt + VLp Adt, откуда получаем
(10) Условие окончания второго этапа: (11)III Для третьего участка кристаллизации (перитектическое превращение при постоянной температуре tр) уравнение баланса энергии принимает вид:
a(t –tep )F dτ = VLpdm, откуда (12)
(13)
Условие окончания третьего этапа: (14)
где mф – количество твёрдой фазы, которая остаётся после завершения перитектической реакции. m(1) легко определяемое по правилу отрезков на диаграмме состояния (ниже линии tр).IV. Для четвертого участка охлаждения твердой фазы в интервале температур от tр до tк уравнение баланса энергии по виду аналогично первому участку:
Условие окончания четвёртого этапа: (16)
После определения всех уравнений модели переходим к выбору метода и составлению алгоритма численного решения задачи.
Выбор метода и разработка алгоритма решения задачи:
Преобразуем дифференциальные уравнения модели в разностную форму методом Эйлера:
(17)(18)
(19) (20)=(17)Программирование задачи:
Оформленная в системе MATLAB программа по описанному выше итерационному алгоритму:
%Programma of coling Fe-0.16C alloy
clear;
c0=0.16;C=444;R=7000;L=277000;a=126.5;F=0.025;V=0.00025;
t(1)=1800;tp=1499;tk=800;tsr=20;dtau=1;m(1)=0.854;n=10000;
tliq=1539-78.43*c0;
%step_1
for i=1:n; if t>tliq;
t(i+1)=t(i)-a*F*dtau*(t(i)-tsr)./(V*C*R);
end;end; s1=length(t);
%step_2
for i=s1:n; if t>tp;
t(i+1)=t(i)+a*F*dtau*(t(i)-tsr)./(-V*C*R-0.0016*V*L*R/(15.77-0.01025*t(i)).^2);
end;end; s2=length(t);
%step_3
m(s2)=0.854; for i=s2:n; if m<1;
m(i+1)=m(i)+a*F*dtau*(t(i)-tsr)./(V*L*R);
t(i+1)=t(i);
end;end; s3=length(t);
%step4
for i=s3:n; if t>tk;
t(i+1)=t(i)-a*F*dtau*(t(i)-tsr)./(V*C*R);
end;end; sk=length(t);
h=0:dtau:dtau*(length(t)-1);
plot(h,t); grid
title('temperature as function of time')
ylabel('temperature,C')
xlabel('time,sec');
Полученная в результате запуска программы математическая модель процесса кристаллизации сплава системы Fe – 0,16 % C
Проверка адекватности модели:
Проверка адекватности программной модели состоит в сопоставлении расчетных ур и экспериментальных уэ значений, которые выдает преподаватель, и в оценке величины коэффициента парной корреляции rУрУз по формуле:
и коэффициента регрессии a=tgφ
где yрi и yэi - пары соответственных значений расчетных и экспериментальных данных;
и — средние значения ури уэпо совокупности всех п сопоставляемых величин.Таблица 1
Анализ согласованности расчётных и экспериментальных данных.
время τ, сек | 50 | 100 | 150 | 200 | 250 | 300 |
t расчётная, °С | 1520 | 1510 | 1490 | 1215 | 1000 | 810 |
t эксперим., °С | 1600 | 1550 | 1450 | 1100 | 850 | 730 |
Рис. 5 Проверка адекватности численной модели с помощью графического образа.
rкр=0,707 (число степеней свободы n=6, уровень доверительной вероятности р=0,95)
Из графического и численного анализа согласования полученных расчетных ури экспериментальных уэданных (рис. 3) можно заключить о принципиальной адекватности полученной модели на основе тесной зависимости между расчетными и экспериментальными данными (так как точки на графике располагаются вблизи биссектрисы угла, образованного координатными осями и r> rкр(0,95;6) и а≈1).
Исследование процесса с помощью модели.
Поскольку результат проверки адекватности модели положителен, проведём серию расчётов по данной модели, варьируя ряд исходных данных.
Изменим коэффициент теплоотдачи α со 126,5 до 1000 Вт/(м2·К):
Увеличим температуру среды до 200 °С (tsr начальная = 20°С):
В 2 раза увеличим размеры слитка с 0,05×0,05×0,1 м3 до 0,1×0,1×0,2 м3
Вывод:
Мы ознакомились с методикой построения детерминированной модели и разработали численную модель металлургического процесса для системы с сосредоточенными параметрами. С помощью коэффициента корреляции мы доказали, что построенная модель является адекватной и отражает реальную картину охлаждения и кристаллизации заданного расплава.
При изменении некоторых параметров (коэффициента теплоотдачи α=1000 Вт/(м2К), температура среды tsr=200°C и размеры слитка 0,1×0,1×0,2 м3) налицо изменения самой модели.
Получили, что при увеличении коэффициента теплоотдачи α примерно в 7,9 раз охлаждение расплава протекает быстрее во столько же раз. А также площадка, соответствующая перитектической реакции, стала ровной, что в большей степени отвечает реальному процессу, поскольку при нонвариантном превращении температура сплава остается неизменной.
При увеличении температуры среды III и IV стадии протекают дольше. Несмотря на то, что tsr является параметром, который присутствует во всех стадиях, её изменение не оказало видимого влияния на замедление скорости протекания I и II стадий.
При увеличении размеров слитка значительно замедлился процесс кристаллизации сплава на всех стадиях. Здесь столь сильное влияние объясняется тем, что увеличив размеры, мы изменили сразу два параметра: поверхность теплообмена F и объем V, которые фигурируют на каждом шаге (на втором шаге V фигурирует дважды ).