Курсовая работа: Математическое моделирование тепловой работы вращающейся печи
где а, b – числовые коэффициенты.
Согласно принятой здесь нумерации точек на гранях контрольных объемов, формулы для расчета коэффициентов дискретных уравнений будут представлены следующим образом:
Коэффициент теплопроводности слоя гарнисажа, образованного на внутренней поверхности футеровки застывшим клинкерным расплавом, принимают равным 1 Вт/(м•К).
В цикле обратной прогонки вычисляют в узлах сетки относительные разности значений температуры в текущей и предыдущей итерациях:
и выбирают из них максимальную разность. В конце общего
итерационного цикла производят оценку сходимости итераций, сравнивая абсолютную
величину 
с заданным малым числом. Если |
| больше чем, например, 0,00001, то итерации
повторяются, если же меньше, то итерации завершаются. Чтобы повторить вновь
выбор максимальной относительной разности температур, в начале каждой итерации
величину 
устанавливают равной нулю.
В конце программы предусматривают вывод на экран и на печать исходных данных и результатов расчета. При этом температуру представляют в градусах Цельсия.
Вторая программа для расчета стационарного температурного поля в футеровке должна учитывать теплообмен футеровки с диффузионным факелом и воздухом в печи, согласно граничным условиям третьего рода. С этой целью в предыдущую программу вносятся соответствующие изменения.
Коэффициенты прогонки Рj, Sj на внутренней поверхности футеровки или гарнисажа определяют теперь, предварительно вычислив коэффициенты дискретного уравнения. При этом используются значения степени черноты и поглощательной способности газообразной среды, приведенной степени черноты и коэффициента конвективной теплоотдачи, найденные при расчете параметров диффузионного факела. Температуру на внутренней поверхности футеровки или гарнисажа определяют обратной прогонкой.
Третью программу составляют для расчета нестационарного теплообмена, имея в виду, что при вращении печи температура внутренней поверхности футеровки изменяется. При контакте с технологическим материалом она равна температуре этого материала, что соответствует граничным условиям первого рода, а при нагреве диффузионным факелом зависит от условий радиационной и конвективной теплоотдачи, согласно граничным условиям третьего рода.
Для учета нестационарности в программе организуют цикл по интервалам времени, который является внешним по отношению к итерационному циклу. Цикл по интервалам времени выполняется в течение нескольких оборотов печи, так чтобы расчетное время прогрева футеровки оказалось достаточным для стабилизации изменений температурного поля в следующих друг за другом оборотах.
Чтобы задать начальное температурное поле в футеровке, целесообразно усреднить стационарные распределения температуры, полученные в двух предыдущих программах при граничных условиях первого и третьего рода. По-прежнему допускается рассматривать задачу приближенно как одномерную. Коэффициент aj, Dj дискретных уравнений рассчитывают теперь с учетом нестационарных членов.
Окружность печи делят на 16, расчетных отрезков и определяют интервал времени, необходимый для прохождения каждого отрезка расчетным сечением футеровки при вращении печи
где nоб – частота вращения печи (1/с).
Если расчетный отрезок футеровки находится под слоем технологического материала, в программе используются граничные условия первого рода. Для остальных участков футеровки расчет выполняется при граничных условиях третьего рода.
Чтобы выделить интервалы времени, когда футеровка находится под слоем технологического материала, разделим порядковый номер расчетного интервала времени на число расчетных отрезков по окружности печи. Те расчетные моменты времени, для которых величина получаете целой, соответствуют началу нового оборота печи. На первом обороте печи это начальный (нулевой) момент времени, а на остальных оборотах печи это моменты времени, порядковые номера которых кратны числу nф.
Если расчетным отрезкам по окружности футеровки, находящимися под
слоем технологического материала, присвоить порядковые номера, то
соответствующие им расчетные моменты времени можно определить из условия, что
отношение 
, является целым числом.
Поперечное сечение вращающейся печи:
1-футеровка; 2 – слой гарнисажа; 3 – технологический материал.
Программа №1
Program Kurs;
{Математическое моделирование тепловой работы вращающейся печи}
uses Crt;
const SI = 5.67e-8; {Постоянная Стефана-Больцмана, Вт/(м2*К4)}
SC = 0.85; {Турбулентное число Шмидта}
type tFurnace = object
X: integer;
Dk, Dw, D0, Hw, Hw1, Tm, Ew1, Ew2: real;
X0, RF, FR, LF, AG, ES, EG: real;
BT, G0, R0, U0, H0, QT, GT: real;
GB, GB0, ROB, TB, HB: real;
ALB, NB, VB0, V1, V2, V3, VG0: real;
TF, TK, TP, QF, QP, HF: real;
Cg, Cv, CF, ROF, ROV, RO: real;
P1, P2, P3, PT, PG, PB: real;
Procedure Date;
Procedure Fakel;
Procedure GasPost;
Procedure GasHeat;
Procedure HeatLine;
Function fHdis:real;
Procedure Names;
Procedure Result;
Procedure Nomina;
end;
var Furnace: tFurnace;
Var File1: text;
{============================================================================}
Procedure tFurnace. Date;
{============================================================================}
{Исходные данные}
Begin
DK:=6.2; {Диаметр корпуса печи, м}
HW:=0.265; {Толщина футеровки, м}
HW1:=HW*0.2; {Толщина гарнисажного слоя, м}
DW:=4.6–2*(HW+HW1); {Внутренний диаметр печи, м}
TM:=1465; {Температура обжигаемого материала,°C}
EW1:=0.7; {Степень черноты внутренней поверхности стенки}
QT:=35510; {Низшая теплота сгорания топлива, кДж/кг}
R0:=0.740; {Плотность топлива, кг/м3 (н.у.)}
BT:=7000; {Расход топлива через горелку, м3/ч (н.у.)}
G0:=BT*R0/3600; {Расход топлива через горелку, кг/с}
U0:=160; {Скорость истечения топлива из горелки, м/с}
D0:=Sqrt (4*G0/R0/U0/Pi); {Диаметр горелки, м}
ALB:=1.12; {Коэффициент избытка воздуха}
VB0:=9.45; {Объем теоретически необходимого воздуха, м3/м3}
ROB:=1.293; {Плотность воздуха, кг/м3 (н.у.)}
NB:=ROB*VB0/R0; {Стехиометрический коэффициент, кг/кг}
GB0:=NB*G0*ALB; {Начальный расход горячего воздуха, кг/с}
TB:=550; {Температура горячего воздуха,°C}
HB:=(1.287+0.0001201*TB)*TB; {Энтальпия горячего воздуха, кДж/кг}
{Объем продуктов горения, м3/м3 топлива}
V1:= 1; {Углекислый газ}
V2:= 2.12; {Водяной пар}
V3:= 7.49; {Азот}
VG0:=10.61; {Продукты стехиометрического горения}
End;
{============================================================================}
Procedure tFurnace. Fakel;
{============================================================================}
{Расчет параметров диффузионного факела}
var GP, ROP, HBP, HHF, HHB, ROD, RS, GTB, MG, ZT: real;
XX, XB, XXB, LX, LXB, SC1, SC2: real;
Begin
Date;
ROF:=R0; ROV:=1.3; CF:=1.5;
Names;
While GB*G0<=GB0 do begin {переход к очередному сечению}
{Расстояние от горелки до расчетного сечения, D0*м}
X:=X+1;
{Параметры факела в предыдущем сечении}
GP:=GB; ROP:=ROF; TP:=TF; QP:=QF;
repeat {начало итераций}
{Средняя плотность факела на очередном малом участке, кг/м3}
ROD:=(ROP+ROF)/2;
{Расход воздуха, вовлеченного в струю, G0*кг/c}
GB:=GP+0.322*sqrt (ROD/R0);
{Средняя плотность факела на расчетном участке, кг/м3}
Ro:=R0*sqr (GB/X/0.322);
RS:=sqrt (Ro/R0); {Соотношение плотностей}
{Расстояние до полюса струи, D0*м}
X0:=2.37/RS;
{Расчетная длина факела, D0*м}
LF:=X0*((NB+1)*(2*SC+1) – 1);
{Расстояние от полюса струи до расчетного сечения и конца факела, D0*м}
XX:=X+X0; LX:=LF+X0;
{Радиус турбулентной струи, м}
RF:=0.211*XX*D0;
{Расход несгоревшего топлива через расчётное сечение, G0*кг/с}
SC2:=2*SC; SC1:=SC2+1;
GT:=XX/X0/NB*(LX/XX/SC1+SC2/SC1*exp (ln(XX/LX)/SC2) – 1);
{Средняя энтальпия газов в расчётном сечении факела, кДж/кг}
HHF:=(H0+HB*GB+QT*(1-GT))/(1+GB); {без теплообмена}
HF:=HHF-QF/(1+GB)/G0; {c теплообменом}
{Средняя температура в расчётном сечении факела,°C}
ZT:=TF; TF:=(HF-fHDis)/CF; TK:=TF+273;
{Плотность газов в расчётном сечении факела, кг/м3}
ROF:=ROV*273/TK;
GasPost; Gasheat; HeatLine;
until ABS (1-ZT/TF)<0.0001; {конец итераций}
{Радиус фронта пламени, м}
FR:=XX/7.9*sqrt (exp(ln (LX/XX)/SC2) 1);
{Вывод результатов на печать}
Result;
end;
{–}
{Параметры факела в начальном сечении зоны догорания}
GTB:=GT; GB:=GB0/G0; HHB:=HHF; RF:=DW/2;
XB:=X; XXB:=XB+X0; LXB:=LX;
MG:=(1-exp (ln(XXB/LXB)/SC2))/NB/X0/GTB;
While GT>0.01 do begin {переход к очередному сечению}
{Расстояние от горелки до расчётного сечения, D0*м}
X:=X+1;
{Параметры факела в предыдущем сечении}
TP:=TF; QP:=QF;
repeat {начало итераций}
{Расход несгоревшего топлива через расчетное сечение, кг/с}
GT:=GTB*exp (MG*(XB-X));
{Средняя энтальпия газов в расчетном сечении факела, кДж/кг}
HHF:=HHB+QT*(GTB-GT)/(1+GB); {без теплообмена}
HF:=HHF-QF/(1+GB)/G0; {c теплообменом}
{Средняя температура в расчетном сечении факела, °C}
ZT:=TF; TF:=(HF-fHDis)/CF; TK:=TF+273;
{Плотность газов в расчетном сечении факела, кг/м3}
ROF:=ROV*273/TK;
Gasheat; HeatLine;
until ABS (1-ZT/TF)<0.0001; {конец итераций}
{Вывод результатов на печать}
Result;
end;
Nomina;
End;
{============================================================================}
Procedure tFurnace. GasPost;
{============================================================================}
{Расчет параметров диффузионного факела}
var RO1, RO2, RO3, ROG, VT, VB, VG, VF: real;
Begin
P1:=V1/VG0; P2:=V2/VG0; P3:=V3/VG0;
RO1:=44/22.41;
RO2:=18/22.41;
RO3:=28/22.41;
{Плотность стехиометрических продуктов горения, кг/м3 (н.у.)}
ROG:=RO1*P1+RO2*P2+RO3*P3;
{Объемный расход несгоревшего топлива, м3/с (н.у.)}
VT:=GT*G0/RO;
{Объемный расход воздуха, м3/с (н.у.)}
VB:=(GB – (1-GT)*NB)*G0/ROB;
{Объемный расход продуктов горения, м3/с (н.у.)}
VG:=(1-GT)*(1+NB)*G0/ROG;
{Объемный расход газообразной среды, м3/с (н.у.)}
VF:=VT+VB+VG;
{Объемные доли компонентов в газообразной среде}
PT:=VT/VF; {Топливо}
PB:=VB/VF; {Воздух}
PG:=VG/VF; {Продукты горение}
{Плотность газообразной среды, м3/кг (н.у.)}
ROV:=RO*PT+ROB*PB+ROG*PG;
End;
{============================================================================}
Procedure tFurnace. GasHeat;
{============================================================================}
{Расчет теплоемкости газообразной среды}
var CB, CT, C1, C2, C3: real;
Begin
{Объёмная теплоёмкость природного газа, кДж/м3*К (н.у.)}
IF TF=0 THEN CT:=2.522+0.0005815*TF ELSE
CT:=(-317.9+2.522*TF+0.0005815*sqr(TF)+86840/(TF+273))/TF;
{Объемная теплоемкость воздуха, кДж/м3*К (н.у.)}
CB:=1.287+0.0001201*TF;
{Объёмная теплоёмкость углекислого газа, кДж/(м3*К) (н.у.)}
IF TF=0 THEN C1:=2.081+0.0002017*TF ELSE
C1:=(-139.5+2.081*TF+0.0002017*sqr(TF)+38110/(TF+273))/TF;
{Объёмная теплоёмкость водяного пара, кДж/(м3*К) (н.у.)}
C2:=1.49+0.0002303*TF;
{Объёмная теплоёмкость азота, кДж/(м3*К) (н.у.)}
C3:=1.28+0.0001103*TF;
{Объёмная теплоёмкость продуктов горения, кДж/м3*К (н.у.)}
CG:=C1*P1+C2*P2+C3*P3;
{Объёмная теплоёмкость газов в факеле, кДж/(м3*К) (н.у.)}
CV:=CT*PT+CB*PB+CG*PG;
{Массовая теплоёмкость газов в факеле, кДж/кг*К}
CF:=CV/ROV;
End;
{============================================================================}
Function tFurnace.fHdis:real;
{============================================================================}
{Расчет теплоты диссоциации факела}
var PCO2, PH2O, KCO2, KH2O, ACO2, AH2O, PS, LR, K1, K2: real;
Begin
fHdis:=0;
if TF>1500 then begin
{Объёмные доли трехатомных газов}
PCO2:=P1*PG; {…углекислого газа}
PH2O:=P2*PG; {…водяного пара}
{Константа равновесия для углекислого газа}
KCO2:=exp (ln(10)*(4.47–14700/TF));
{Константа равновесия для водяного пара}
KH2O:=exp (ln(10)*(3.05–13160/TF));
{Степень диссоциации углекислого газа}
ACO2:=exp (ln(2*sqr(KCO2)/PCO2)/3);
{Степень диссоциации водяного пара}
AH2O:=exp (ln(2*sqr(KH2O)/PH2O)/3);
{Теплота диссоциации трехатомных газов, кДж/кг}
fHdis:=(12630*ACO2*PCO2+10800*AH2O*PH2O)/ROV;
end;
End;
{============================================================================}
Procedure tFurnace. HeatLine;
{============================================================================}
{Расчет результирующего излучения факела}
var TG, FF, FW, PC, PH, PS, LR, K1, K2: real;
Begin
{Средняя температура факела на малом участке,°C}
TG:=(TP+TF) /2;
{Расчетная поверхность малого участка факела, м2}
FF:=2*PI*RF*D0;
{Внутренняя поверхность стенки печи, м2}
FW:=PI*DW*D0;
{Объёмные доли излучающих газов}
PC:=P1*PG+PT; {углекислый газ + топливо}
PH:=P2*PG; {водяной пар}
PS:=PC+PH; {трёхатомные газы}
{Эффективная толщина излучающего слоя, м}
LR:=1.8*RF;
{Коэффициент поглощения при температуре газов, 1/м}
K1:=PS*((0.78+1.6*PH)/sqrt (PS*LR) – 0.1)*(1–0.37*(TG+273)/1000);
{Степень черноты факела}
EG:=1-exp (-K1*LR);
{Коэффициент поглощения при температуре стенки, 1/м}
K2:=PS*((0.78+1.6*PH)/sqrt (PS*LR) – 0.1)*(1–0.37*(TM+273)/1000);
{Поглощательная способность факела}
AG:=1-exp (-K2*LR);
{Приведенная степень черноты}
ES:=1/(1/AG+(1/EW1–1)*FF/FW);
{Поток результирующего излучения факела, кВт}
QF:=QP+ES*SI*(EG/AG*sqr (sqr(TG+273)) – sqr (sqr(TM+273)))*FF/1000;
End;
{============================================================================}
Procedure tFurnace. Names;
{============================================================================}
Begin
Assign (File1,'Furnace.txt'); Rewrite(File1);
Writeln(File1); Writeln (file1,' ИСХОДНЫЕ ДАННЫЕ:');
Writeln (file1,' Внутренний диаметр печи DW:=', DW:5:2,' м');
Writeln (file1,' Toлщина футеровки HW:=', HW*1000:4:0,' мм');
Writeln (file1,' Температура обжигаемого материала TM:=', TM:5:0,'°C');
Writeln (file1,' Низшая теплота сгорания топлива QT:=', QT:5:0,' кДж/кг');
Writeln (file1,' Плотность топлива R0:=', R0:3:3,' кг/м3 (н.у)');
Writeln (file1,' Расход топлива через горелку G0:=', G0:5:3,' кг/с');
Writeln (file1,' Скорость истечения топлива из горелки U0:=', U0:5:1,' м/с');
Writeln (file1,' Диаметр горелки D0:=', D0*1000:3:0,' мм');
Writeln (file1,' Коэффициент избытка воздуха ALB:=', ALB:3:2);
Writeln (file1,' Температура горячего воздуха TB:=', TB:4:0,'°C');
Writeln(File1);
Writeln (file1,' X LF GT GB EG AG ES HF TF QF ');
End;
{============================================================================}
Procedure tFurnace. Result;
{============================================================================}
Begin
Append(File1);
if X mod 2 = 0 then
Writeln (file1, x:3,' ', lf:5:1,' ', gt:6:4,' ', gb:6:3,' '
eg:5:3,' ', ag:5:3,' ', es:5:3,' ', hf:6:0,' ', tf:6:0,' ', qf:5:0);
End;
