Построение детерминированной программой математической модели кристаллизации сплава системы Fe-C

Федеральное агентство по образованию

ГОУ «Санкт-Петербургский  государственный политехнический университет»

Факультет технологии и исследования материалов

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























Отчёт

о лабораторной работе №1 «Построение детерминированной программной модели».

























Работу выполнила студентка группы 4064/1                           		Петрова С.С.

Преподаватель                                                                             		Вяххи И.Э.



























Санкт-Петербург

2009 г.

Цель работы – знакомство с методикой построения детерминированной модели и разработка численной модели металлургического процесса для системы с сосредоточенными параметрами.

Система – целесообразная совокупность взаимодействующих элементов, которая ориентирована на выполнение той или иной функции. Системы классифицируются по различным параметрам. По пространственной структуре бывают системы с сосредоточенными и распределёнными параметрами.

Система с сосредоточенными параметрами имеет равномерное распределение выходных параметров по объёму и характеризуется их осреднёнными величинами.

Система с распределёнными параметрами имеет более сложную структуру, в ней выходные параметры неравномерно распределены по объёму.



Модель – объект, находящийся по отношению к натурному объекту в отношении подобия (т.е. взаимно-однозначного соответствия); приближённое описание процессов, происходящих в системе, ориентированное на выполнение определённых функций.

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

Постановка задачи:

Необходимо разработать математическую модель процесса кристаллизации сплава системы Fe – 0,16 % C с учётом равновесного выделения твердой фазы по диаграмме состояния. Изменение температуры данного сплава происходит в интервале от tз=1800°С (температуры заливки в форму) до tк=800°С (конечной температуры охлаждения слитка). Удельная теплоёмкость сплава С=444 Дж/(кг·К), его плотность ρ=7000 кг/м3, коэффициент теплоотдачи α=126,5 Вт/(м2·К), скрытая теплота фазового превращения L=277 кДж/кг.







Поскольку мы имеем дело с системой с сосредоточенными параметрами, то модель, с помощью которой мы будем описывать процессы, происходящие в системе, так же будет с сосредоточенными параметрами, то есть будет характеризоваться осреднёнными параметрами.







Построение физической модели:

На рисунке 2 показан характер изменения температуры во времени при охлаждении сплава заданного состава. Разобьем температурный интервал охлаждения слитка сплава  Fe – 0,16 % C на участки:































	Участок жидкого состояния.

	При заданной температуре заливки tз=1800°С сплав Fe – 0,16 % C находится в жидком состоянии; жидкость охлаждается до пересечения с линией А-В при температуре ликвидуса (tliq) в точке А1.

	 В точке А1 начинается выпадение из расплава кристаллов феррита, 

							  Ж → Ф

	при этом концентрация жидкости меняется по линии А-В, а концентрация феррита – по линии А-С.

	III. В точке С1, когда сплав достигает температуры перитектики (tp=1499°С), выпадение кристаллов феррита заканчивается. При данной температуре происходит перитектическое превращение:

							Ж + Ф → А

	С-С1-В – линия нонвариантного перитектического превращения.

	IV. Ниже точки С1 охлаждение идёт в области диаграммы, соответствующей твёрдому состоянию сплава, до заданной температуры охлаждения слитка  tк=800°С.

		Для дальнейшего построения физической модели примем следующие допущения:                          

		для простоты модели будем рассматривать систему как объём с сосредоточенными параметрами, т. е. не учитывая перепада температур по сечению слитка и принимая его в качестве материальной точки, имеющие постоянные поверхность теплообмена F, объем V и плотность ρ;













На рисунке 3 изображён слиток заданного сплава:



10 см

5 см

5 см

Рис. 3 Нагреваемый слиток металла















(м2), (м3).  

		теплообмен слитка со средой происходит по закону конвективного теплообмена Ньютона с постоянным коэффициентом теплоотдачи α, т.е. пренебрежём лучистым (радиационным) теплообменом;

		основные параметры системы (плотность р, теплоемкость С, теплота кристаллизации L) являются постоянными, не зависящими от температуры и состава выделяющихся фаз, в том числе не учитываются объёмные изменения и физико-химическое взаимодействие между сплавом и окружающей средой; 

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



Формулировка математической модели.



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

, (1)

где α – коэффициент теплоотдачи от поверхности слитка площадью F и объёмом V в окружающую среду; 

tcp - температура окружающей среды; 

t - температура слитка; 

τ - время.

В общем случае изменение внутренней энергии сплава в зависимости от этапа кристаллизации имеет вид

	

Для I и IV участков, где  dm=0,

для II участка, 

для III участкa, где  dt=0,	(2)

	

	где dm – изменение относительного количества твердой фазы ( т ≤ 1).                Рассмотрим участки охлаждения и кристаллизации сплава Fe – 0,16 % C 

I. Для первого участка охлаждения жидкой фазы от температуры заливки (tз) до температуры ликвидуса этого сплава (tn) уравнение баланса энергии имеет вид обыкновенного дифференциального уравнения первого порядка:

		, откуда 	.(3)	

	

	Условие окончания первого этапа охлаждения жидкого металла имеет вид:

		

	

	

	Значение tл(С0) рассчитываем по диаграмме состояния в виде линейного соотношения

		, где С0 – состав сплава.

	Зависимость температуры ликвидуса от состава сплава при допущении о прямолинейности линий диаграммы состояния находим в виде уравнения прямой

	по двум точкам. Как видно из диаграммы состояния (рис. 1), на линии А-В при С=0% tл=1539°С; при С=0,51% tл=1499°С. Составим систему уравнений с использованием данных значений:

												 (4)

	

	

	

	

	Решив эту систему, получили уравнение зависимости температуры ликвидуса от состава сплава: 

               tл =1539–78,43С 			(5)

II. Для второго участка охлаждения в интервале температур от перитектики (tp) уравнение баланса энергии с учетом теплоты фазового превращения при выделении феррита имеет вид:

a (t–tep)Fdτ = –VСpdt + VLpdm,                  (6)

где в правой части содержатся дифференциалы двух взаимозависимых переменных – t и m. Для исключения «лишней» переменной и преобразования этого уравнения к виду

  



необходимо связать дополнительным соотношением t и m , а затем выразить dm через количество твердой фазы – феррита при температуре t в условиях равновесной кристаллизации найти из диаграммы состояния по правилу отрезков:



		,                                       (7)

	

где  С1(t) и С2(t) – концентрация углерода С в жидком сплаве (линия А-В) и в феррите (линия А-С) при температуре t.. 

	Зависимости С1(t) и С2(t) находим, используя допущение о том, что линии диаграммы состояния являются прямыми. Как видно из диаграммы состояния (рис. 1), на линии А-В при t=1539°С С=0%; при t=1499°С С=0,51%. Составим систему уравнений с использованием данных значений:												        (8)

	

	

	а = – 0,01275;  

	

	Решив эту систему, получили уравнение зависимости С1(t)=19,62 – 0,01275t .

		Аналогично находим уравнение для прямой С2:

	На линии А-С при t=1539°С С=0%; при t=1499°С С=0,1%. Составим систему уравнений:													    (9)

	

	 Решив эту систему, получили уравнение зависимости С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к уравнение баланса энергии по виду аналогично первому участку:

		

		,  	 (15)	

	

	Условие окончания четвёртого этапа:	 (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 фигурирует дважды ).











2