基于遗传算法及温度预报模型参数优化

基于遗传算法的温度预报模型参数优化
1问题描述
近年来，随着纯净钢生产技术的进步和连铸技术的发展，炉外精炼工艺与设备迅速普及。

其中，LF以其优异的综合性能，在实际生产中得到了广泛应用。

而点测钢水温度的高成本，使精炼炉温度预报成为了极具实际意义的研究。

因此，钢水温度预报模型的建立显得至关重要。

目前，温度预报模型的建立主要采用3种方法，即：机理模型、“黑箱模型”和“灰箱模型”。

机理模型是指用尽可能准确的数学方程来描述过程机理而建立的模型，而“黑箱模型”则采用一些数学方法（智能算法、回归算法等），结合实验数据或实际生产中的数据，生成一种只有输入和输出的模型，而完全不考虑过程机理。

然而，机理模型需要考虑的因素太多，且这些因素具有相当的不确定性，因此模型中的许多参数很难得到，严重影响了温度预测的精度；“黑箱模型”则完全建立在数据的基础上，如果生产环境和条件改变导致了数据改变，原先模型的准确性就得不到保障，即模型的可移植性差。

因此，本文主要研究“灰箱模型”的建立，即采用机理建模与数据建模相结合的方式，应用智能优化算法对机理模型中较难获得的参数进行辨识与确定，其流程如下图所示：
图1 温度预报模型建立的方案流程图
而参数的辨识过程，即是基于智能优化算法，根据输入输出条件，对机理模型不易获得的参数进行寻优的过程。

2 理论基础
2.1 机理模型
本研究所基于的机理模型，是根据LF 炉精炼过程中的传热基本方程、能量守恒方程和质量守恒方程等来建立的。

并运用有限差分法和有限元法等，通过控制初始条件和边界条件来对模型进行求解，得到钢包内的温度情况。

2.1.1 热平衡分析
3.1.1 吸热与散热
精炼过程中，钢水热量的来源与去向大致如图2所示。

图2 热量的来源与去向
将钢水和炉渣作为一个系统，来推导吸热和散热与其温度变化的关系。

由于系统在加热与散热的过程中，其内能变化都体现在温度的变化上，所以系统实际吸收的热量为通过电弧加热所吸收的热量与冶炼过程中散去的热量之差，即：
()st st sl sl e sa ch ls c m c m T Q Q Q Q +∆=--- （1）
式中，st c 、sl c 为钢水和渣的比热容，st m 、sl m 为钢水和渣的重量，T ∆为系统总的温度变化，e Q 、sa Q 、ls Q 、ch Q 分别代表通过电弧加热所吸收的热量，和通过渣面及炉盖散热、炉壁及炉底散热、加渣料及加合金影响所损失的热量。

为最终求出温度随时间变化的曲线，可对式（1）两端对时间取微分。

由此，对每个温度的影响因素进行分析与建模，分别求解它们的热量变化速率，再将模型综合后求积分，即可得出温度随时间变化规律的曲线，温度的计算公式为：
()()e sa ch ls st st sl sl dQ dQ dQ dQ T c m c m dt dt dt dt dt
=---+⎰ （2）
2.1.2 机理模型描述
（1）电弧加热模型
系统所吸收的热量应与总的耗电量相关，二者间是一个乘系数的关系，即：
e arc all Q Q η= （3）
其中all Q 为该次冶炼过程中的总耗电量，而arc η是系统的能量吸收系数，即电弧效率。

（2）渣面及炉盖散热模型
LF 精炼过程中，通过渣面损失的热量sa Q 包括通过渣面的辐射和对流的热损失sr Q ，还包括吹氩带走的热量g Q ，即：
sa sr g Q Q Q =+ （4）
渣面及炉盖散热的近似计算方法，是将渣面散热和炉盖散热视为一个整体，通过计算水冷炉盖中冷却水带走的热量，来确定渣面及炉盖的散热量，即有：
()sr w w w all wout win Q c G t T T ρ=- （5）
式中，w c 、w ρ分别为冷却水的比热容和密度，w G 为冷却水流量，其单位为3/m h ；all t 为冶炼时间，单位为h ；wout T 、win T 分别为冷却水的出口、入口温度。

而计算高温气体带走热量的经验计算方法为：对小于120t 吨位的钢包炉，气体带走热量约占有功功率的5%-6%；而对于大于120t 吨位的钢包炉，气体带走的热量大约占有功功率的7%-9%。

因此有：
g g Q P η= （6）
其中g η为气体带走热量占有功功率的百分数，P 为有功功率。

将sr Q 和g Q 综合考虑，并考虑精炼时间和吹氩时间，以冷却水流量、冷却水出入口温度，以及有功功率作为输入，以g η为模型待确定参数，即可搭建该部分模型。

（3）炉壁及炉底散热模型
炉壁散热可当作一维非稳态导热问题建模，利用其初始条件和边界条件来求解。

为简化问题，将炉壁的倾斜度忽略不计，将其视为圆柱体；并且将炉壁材料的导热系数和比热容取平均值，视为常量。

计算时，可将炉壁视为无限长圆筒壁，在柱坐标下求解。

就炉壁而言，其导热微分方程式为：
22()c c c p p p T T T c t r r r
ρλ∂∂∂=+∂∂∂ （7） 上式中，p ρ、p c 、p λ分别为钢包材料的密度、比热容和导热系数；c T 为炉壁径向上某点的温度，r 为该点距钢包中轴线的距离（即半径）；t 为时间，单位是s 。

式（7）的初始条件（0t =）为：
020021ln(/)ln(/)
c ls st ls T T r r T T r r -=- （8） 式中，0st T 、0ls T 分别为首次测温时钢水温度和钢包外壁的温度；1r 、2r 分别为钢包的内径和外径。

式（8）的边界条件是，在1r r =时有：
c st T T = （9）
在2r r =有：
2|()c p r r ac c a T h T T r
λ=∂-=-∂ （10） 以上两式中，st T 、a T 为钢水和环境的实时温度；ac h 为钢包侧壁的综合对流换热系数，其单位为2/()W m C ⋅。

运用有限差分的方法，可对该微分方程进行求解。

同理，对于炉底散热问题，也可用同样的方法进行讨论，这里不再赘述。

（4）加渣料与合金的影响
加料是LF 的精炼流程中的一个重要环节，不同渣料的熔化热不同，因此加渣料对于钢水温度的影响取决于不同渣料的热物性参数和加入量的多少；与加渣料类似，不同合金的熔化热、熔解热以及化学反应吸收或放出的热量都不同。

根据大量的实际数据统计，以下给出一些常见渣料、合金的温降系数值，如表1所示。

其中，加渣料、合金对钢水的平均温降系数i k ，单位为/100C kg 。

表1 部分常见渣料与合金的温降系数
需要注意的是，上表中的负值表示降温。

于是，加渣料与合金的热效应可由式（11）计算。

100
ch i i
melt
dQ k m
dt t
=∑（11）
其中
i
m为加某种渣料或合金的重量（kg），
melt
t为熔化时间。

2.1.3 机理模型仿真
依据各部分机理模型，以其所有输出与钢水、渣的重量及比热容作为总模型输入，最终通过仿真输出钢水温度。

该模型的Simulink框图及仿真结果如下所示。

图3 总机理模型Simulink仿真框图
图4 总机理模型Simulink仿真结果
上图中，横坐标表示时间（s），纵坐标表示钢水温度（℃），该图显示了一个精炼周期内（约45min）钢水温度的变化情况。

显然，在电弧加热的时间段内，钢水温度显著上升；在其它时段，由于散热，钢水温度有缓慢下降的趋势。

而且在一个精炼周期内，钢水大约升温70C，基本符合实际情况。

然而，由于电弧效率、钢包材料的导热系数与比热容等参数很难准确得知，机理模型中只能大致在一定范围内取值，因此机理模型仅能准确反映温度变化趋势，难以达到精确预测温度的效果。

因此，本研究将利用遗传算法来对这些难以准确获得的参数进行辨识与优化。

2.2 遗传算法
2.2.1 算法简介
遗传算法（GA）是基于Darwin进化论和Mendal遗传学说的一种优化搜索方法，从20世纪60年代开始兴起。

它是基于自然选择和基因遗传学原理的搜索方法，将“优胜劣汰，适者生存”的生物进化原理引入待优化参数形成的编码串群体中，按照一定的适配值函数及一系列遗传操作对各个个体进行筛选，从而使适配值高的个体被保留下来，组成新的群体。

新群体包含上一代的大量信息，并且引入了新的优于上一代的个体。

这样周而复始，群体中各个体适应度不断提高，直至满足一定条件。

此时，群体中适配值最高的个体即为待优化参数的最优解。

遗传算法主要有以下特点：
（1）遗传算法的处理对象不是参变量本身，而是参变量编码后的称为人工染色体的位串，使得遗传算法可直接对集合、队列、矩阵、图表等结构对象进行操作。

这一点使得遗传算法的应用范围很广。

（2）遗传算法是多点搜索，而不是单点，从而避免陷入局部最优，逐步逼近全局最优解。

也正是其固有的并行性，是遗传算法优于其他算法的关键。

（3）与自然界类似，遗传算法对求解问题的本身可以一无所知，如对搜索空间没有特殊要求（如连续、凸性等），对目标函数也不要求诸如连续性、导数存在和单峰等假设，它所需要的仅是对算法所产生的每个染色体进行评价，并基于适应值排列等级来选择染色体，使适应值好的染色体比适应值差的染色体有更多的繁殖机会。

遗传算法利用简单的编码技术和繁殖机制来表现复杂的现象，从而解决非常困难的问题。

（4）遗传算法是一种自适应随机搜索方法。

遗传算法是以概率原则指导搜索，适应值高的个体，在进化过程中将被赋予更高的选择概率，在下一代中有更大繁殖机会。

遗传算法可以更有效地利用已有的信息来搜寻那些有希望改善解质量的串，比一般的随机搜索有更高的搜索效率。

因此，利用遗传算法的上述优点与广泛适用性，对机理模型中难以获得的参数加以辨识后，再应用到机理模型中，是一种能进一步发展和完善机理模型的切实可靠的方法。

2.2.2 编码方法
编码是应用遗传算法时要解决的首要问题，也是设计遗传算法时的一个关键步骤。

编码除了决定个体的染色体排列形式之外，还决定了个体从搜索空间的基因型变化到解空间的表现型时的解码方法，编码方法也影响到交叉算子、变异算子等遗传算子的运算方法。

迄今为止，最常用的是二进制编码方法和浮点数编码方法。

本例中采用浮点数编码方法，即将个体的每一个基因值用某一范围内的一个浮点数来表示，个体编码长度等于其决策变量的个数。

因为这种编码方法使用的是决策变量的真实值，所以浮点数编码方法也称真值编码方法或实参编码方法。

浮点数编码方法大致具有以下几个优点：
（1）适合于在遗传算法中表示范围较大的数。

（2）适合于精度要求较高的遗传算法。

（3）便于较大空间的遗传搜索。

（4）改善了遗传算法的计算复杂性，提高了运算效率。

（5）便于设计针对问题的专门知识的知识型遗传算子。

根据之前建立起的机理模型，选择系统中不易获得的参数进行辨识，它们是：电弧效率arc η、钢包材料导热系数p λ、钢包材料比热容p c 。

将这些参数作为遗传算法中组成种群中个体的基因，并以向量(,,)T arc p p c ηλ来表示。

采用浮点数编码方法，即用这些参数的真值作为其编码，并且限定这些参数的范围。

根据经验，给出电弧效率、钢包材料导热系数以及比热容的参考范围，即：
[0.5,0.8]arc η∈，[1.3,6.5]p λ∈，[1000,1380]p c ∈ （12）
2.2.3 适应度函数的建立
个体适应度函数的值不仅代表了该个体本身的好坏，还决定了在下一代遗传中，该个体存留下来的概率。

本研究利用模型计算的终点温度值与实际测量的终点温度值间的误差来构造适应度函数，即：
2
1
1
(')i n j j j f T T ==-∑ （13）
式中，i f 表示某个体的适应度，n 为输入输出数据的组数，而'j j T T -为第j 组数据的计算温度值与实际温度值的误差值。

为方便计算，将图3所示的Simulink 机理模型编写为一个m 函数，供遗传算法所调用。

该函数以总耗电量、加热时间、有功功率、钢包尺寸、加料量、钢水初始温度等信息为输入，用(,,)T arc p p c ηλ作为待优化参数，以钢水终点温度为函数计算的输出，以钢水的终点测温值作为比较。

而输入信息和钢水的终点测温值，均取自多炉次的现场操作记录与数据。

2.2.4 基本操作
（1）选择操作
选择是从一个旧种群中选择生命力强的个体位串，产生新种群的过程。

本例中选择操作采用适度比例法，它根据某个体的适应度与该代全部个体适应度之和的比值，来决定其子孙遗留的可能性，即在第n 代中，某个体i 被选取的概率为：
i si N
f P f =∑ （14） 其中，i f 为某个体的适应度，N f ∑为该代全部个体的适应度之和，可根据式（13）来进行计算。

确定概率后，可用轮盘赌的选择方法来实现选择操作。

例如，用计算机生成0～1之间均匀分布的随机数，若0.5si P =，则当产生的随机数在0～0.5之间时，该串被复制，否则被淘汰。

（2）交叉操作
交叉运算是指两个相互配对的染色体按某种方式相互交换某部分基因，从而形成两个新的个体。

在遗传算法中，将种群中的N 个个体以随机的方式组成/2N 对配对个体组，交叉操作在这些配对个体组中的两个个体间进行的。

考察某配对个体组中的两个个体a 、b ，交叉操作采用一定方式将它们变为两个新的个体'a 、'b 。

在遗传算法中，交叉操作过程需要满足：
''a b a b +=+ （15）
基于上式，浮点数编码的交叉操作采用如下方式来实现：
'(1)'(1)a a b b b a αββα=-+⎧⎨=-+⎩
（16） 其中，α、β为区间(0,)r 上均匀分布的随机数，且有1r ≤，调整r 的大小即可控制交叉操作的变化范围。

（3）变异操作
使用变异算子来调整个体编码串中的部分基因值，可以从局部的角度出发使个体更加逼近最优解，提高遗传算法的局部搜索能力。

变异操作将某个个体的参
数c ，操作变为域内的另一个值'c 。

浮点数编码的遗传算法采用下式进行：
'(,)c N c σ= （17）
其中，(,)N c σ表示平均值为c ，方差为σ的正态分布随机数。

可以看出，变异操作即是以当前值为中心，主要在一个小范围内进行随机扰动的变化。

3 算法设计
根据以上介绍，就遗传算法本身而言，其步骤主要如下：
（1）确定编码方式，选取一定大小的初始种群；
（2）计算所选种群中每一个个体的适应度函数值及复制概率；
（3）用轮盘赌方法，淘汰相应个数的函数值较小的个体，替换为相应个数的函数值较大的个体；
（4）对个体随机两两配对，按指定概率c P 进行交叉操作；
（5）对每一个体中的每一参数，按指定概率m P 进行变异操作；
（6）若满足收敛条件则输出最优解并退出，否则返回执行步骤（2），如此往复。

而根据本机理模型具体实例而言，现将具体的算法流程作以介绍。

步骤一：产生待辨识参数的初始种群
根据式（12）中所示的待辨识参数上下限范围，采用随机的方式，生成3个参数的初始种群，种群大小为100N =，并设置最大进化代数100G =。

若
(,,)T arc p p c ηλ为种群中某个体的向量表示，
则用MATLAB 随机生成的种群可用一个3100⨯的矩阵来表示。

步骤二：计算个体的适应度
选取50炉钢的操作记录（即数据组数50n =）作为输入数据及温度比较标准，
对于选中的每个个体，通过n组输入数据的计算结果，根据式（13）计算适应度函数的值，并判断当前的进化代数是否达到设定值。

若达到设定值，则结束遗传算法寻优过程，此时最大适应度值对应的个体即为机理模型参数的辨识结果；若进化代数未达到设定值，则继续执行步骤三。

步骤三：进行遗传算法操作
根据步骤二中计算出的适应度值与式（14）计算出的选择概率以及预先设定的交叉和变异概率，进行遗传算法的选择、交叉和变异操作，产生下一代群体，而后将遗传代数G增加1再返回步骤二。

根据以上的三个步骤，绘制算法的流程图如下所示：
图5 算法流程图
4仿真实验与结果分析
综合前文所述，参与参数优化的输入数据有：
总耗电量，冷却水流量，冷却水入口、出口温度，有功功率，精炼周期，加渣、加合金重量，钢包内外径，钢包高度，包底厚度，钢水初温，钢水进站与出站重量等。

利用MATLAB中英国谢菲尔德大学遗传算法工具箱，很容易利用工具箱内
的函数进行生成种群、选择、交叉、变异操作。

经MATLAB 运行仿真，最终遗传代数到达100代时得到种群中的最优个体如下：
(,,)(0.7352,4.132,1223)T T arc p p c ηλ= （18）
而随着遗传代数的增加，种群中的个体逐步优化，计算温度与实际温度的误差减小，即个体适应度增大，图6为每代最优个体适应度值的变化趋势。

图6 误差趋势图
可见，该趋势符合实际情况。

得到参数后，将最终得到优化的参数代入到钢水温度预报的机理模型中，即可得到更加精确的温度预报曲线。

为检验模型的精确性，另取20组现场数据，以式（18）所给出的参数计算终点温度，并与实测温度比较，结果如下表所示。

表2 模型验证结果
为使结果更加直观，用折线图的方式绘制温度计算值与实测值的偏差如下：
图7 预报温度与实际温度误差折线图
可见，在20组数据中，有11组的预报结果误差在5℃以内，有16组的预报结果误差在10℃以内，而全部预报结果的误差均在15℃以内，平均预报温度误差为6.2℃。

可以看出，精炼过程中钢水温度的预报值与实测值比较吻合，说明该模型建立即参数优化比较合理。

5总结与展望
本文在钢水温度预报机理模型的基础上，应用遗传算法对机理模型中不易获得的参数进行辨识与寻优，取得了初步的成功。

然而，在温度预报方面，虽然本文中的机理模型可以正确预测温度变化趋势，混合模型可以利用智能方法辨识出部分机理模型的参数，使模型更加准确，但由于冶炼现场的复杂性，任何确定性或非确定性因素都有可能影响温度变化的走势。

因此，若想建立更加精确的温度预报模型，还需要对现场的实际情况做更加深入的考察与分析。

此外，应用遗传算法时，待辨识参数的选择也尤为关键。

若选择较多的参数进行辨识，有可能使辨识后模型的结果更贴近于实际情况，但计算相对复杂；而选择较少的参数辨识，可使计算较为简化，但模型精度变差。

因此，结合现场经验，合理把握辨识参数的程度亦十分关键。

而在执行遗传算法的过程中，若种群中个体的个数为N，数据的组数为n，则每代种群需进行n N
次计算，效率相对较低，所以研究并开发其它相对高效的智能算法做钢水温度预报研究，将成为我们今后研究方向之一。

附录 MATLAB源程序代码
clc
clear all
data=xlsread('C:\wrydata.xls');
P=data(:,1:end-1);
T=data(:,end);
N=3;
threshold=[-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1 1;-1
1];%******************************************************************************* ******************
%%ÒÅ´«Ëã·¨²ÎÊý
NIND=100;%ÖÖȺ´óС
MAXGEN=100;
PRECI=3;%************************************************************************* ************************¸öÌ峤¶È
GGAP=0.95;
px=0.7;
pm=0.01;
trace=zeros(N+1,MAXGEN);
Chrom=crtrp(NIND,[0.5 1.3 1000;0.8 6.5 1380]);
%%ÓÅ»¯
gen=0;
for k=1:MAXGEN
Obj=zeros(size(P,1),NIND);
tin=25;R1=1.85;R2=2;Tevn=30;Height=4;Z1=0.25;
for j=1:NIND
for i=1:size(P,1)
qall=P(i,1);
iter=Chrom(j,1);
flow=P(i,2);
tout=P(i,3);
pactive=P(i,4);
RealTime=P(i,5);
Cp=Chrom(j,3);
lamad=Chrom(j,2);
zha=P(i,6);
cao=P(i,7);
casi=P(i,8);
sife=P(i,9);
Tst0=P(i,11);
Tls0=P(i,12);
mstin=P(i,13);
mstout=P(i,14);
Obj(i,j)=mechanicalmodel(qall,iter,flow,tout,tin,pactive,RealTime,Cp,lamad,zha,cao,casi,sife,mnfe,R1, R2,Tst0,Tls0,Tevn,Height,Z1,mstin,mstout);
end
end
ObjV=zeros(1,NIND);
fprintf('%d\n',gen);
for j=1:NIND
ObjV(j)=sum((Obj(:,j)-T).^2);
end
ObjV=ObjV';
FitnV=ranking(ObjV);
SelCh=select('sus',Chrom,FitnV,GGAP);
SelCh=recombin('reclin',SelCh,px);
ObjSel=zeros(size(P,1),size(SelCh,1));
for j=1:size(SelCh,1)
for i=1:size(P,1)
qall=P(i,1);
iter=SelCh(j,1);
flow=P(i,2);
tout=P(i,3);
pactive=P(i,4);
RealTime=P(i,5);
Cp=SelCh(j,3);
lamad=SelCh(j,2);
zha=P(i,6);
cao=P(i,7);
casi=P(i,8);
sife=P(i,9);
mnfe=P(i,10);
Tst0=P(i,11);
Tls0=P(i,12);
mstin=P(i,13);
ObjSel(i,j)=mechanicalmodel(qall,iter,flow,tout,tin,pactive,RealTime,Cp,lamad,zha,cao,casi,sife,mnfe, R1,R2,Tst0,Tls0,Tevn,Height,Z1,mstin,mstout);
end
end
ObjSelV=zeros(1,size(SelCh,1));
for j=1:size(SelCh,1)
ObjSelV(j)=sum((ObjSel(:,j)-T).^2);
end
ObjSelV=ObjSelV';
[Chrom,ObjV]=reins(Chrom,SelCh,1,1,ObjV,ObjSelV);
gen=gen+1;
[Y,I]=min(ObjV);
trace(1:N,gen)=Chrom(I,:);
trace(end,gen)=Y;
end
%%½ø»¯Í¼
outcome=trace;
figure(1);
plot(1:MAXGEN,trace(end,:));
grid on;
xlabel('ÒÅ´«´úÊý')
ylabel('Îó²î±ä»¯')
title('½ø»¯¹ý³Ì')
bestX=trace(1:end-1,end);
bestErr=trace(end,end);
function
y=mechanicalmodel(qall,iter,flow,tout,tin,pactive,RealTime,Cp,lamad,zha,cao,casi,sife,mnfe,R1,R2,Ts t0,Tls0,Tevn,Height,Z1,mstin,mstout)
y1=absorbedheat(qall,iter);
y2=surfacedissipation(flow,tout,tin,pactive,RealTime);
y3=materials(zha,cao,casi,sife,mnfe);
y4=sidedissipation(R1,R2,Tst0,Tls0,Tevn,Height,RealTime,Cp,lamad);
y5=bottomdissipation(Z1,Tst0,Tls0,Tevn,R1,RealTime,Cp,lamad);
y=(y1-y2-y4-y5)/(460*mstin+1050*(mstout-mstin))-y3+Tst0;
function y=absorbedheat(qall,iter)
y=qall*iter*3.6*1000000;
function y=surfacedissipation(flow,tout,tin,pactive,RealTime)
y=RealTime/3600*flow*1000*4200*(tout-tin)+pactive*0.065*RealTime/3600*3.6*1000000;
function y=materials(zha,cao,casi,sife,mnfe)
y=(zha+cao+casi*1.05-sife*0.9+mnfe*0.9)*0.01;
function y=sidedissipation(R1,R2,Tst0,Tls0,Tevn,Height,RealTime,Cp,lamad)
deltaR=(R2-R1)/9;
for i=1:1:10
R(i)=R1+deltaR*(i-1);
end
T0=zeros(1,10);
for i=1:1:10
T0(i)=Tls0+(Tst0-Tls0)*((log(R2/R(i)))/(log(R2/R1)));
end
lamad1=lamad;
deltaTau=10;hac=10;rho=2000;B=hac*deltaR/lamad1;
Fo=(lamad1*deltaTau)/(rho*Cp*deltaR*deltaR);
T=zeros(270,10);
for j=1:1:270
for i=1:1:10
switch i
case 1
T(j,i)=Tst0;
case 10
T(j,i)=(1-2*Fo*(1-deltaR/(2*R(i)))-2*Fo*B)*T0(i)+2*Fo*(1-deltaR/(2*R(i)))*T0(i-1)+2*Fo*B*Tevn;
otherwise
T(j,i)=Fo*(1-deltaR/(2*R(i)))*T0(i-1)+(1-2*Fo)*T0(i)+Fo*(1+deltaR/(2*R(i)))*T0(i+1);
end
end
T0=T(j,:);
end
for j=1:1:270
Qac(j)=lamad1*(T(j,1)-T(j,2))/deltaR*2*pi*R1*Height;
end
x=RealTime/10;
if (x==0)
Result=0;
else
Result=Qac(fix(x));
end
y=Result*RealTime;
function y=bottomdissipation(Z1,Tst0,Tls0,Tevn,R1,RealTime,Cp,lamad) deltaZ=Z1/9;
for i=1:1:10
Z(i)=deltaZ*(i-1);
end
T0=zeros(1,10);
for i=1:1:10
T0(i)=Tst0-(Tst0-Tls0)*Z(i)/Z1;
end
lamad2=lamad;
deltaTau=10;hab=5;rho=2000;B=hab*deltaZ/lamad2;
Fo=(lamad2*deltaTau)/(rho*Cp*deltaZ*deltaZ);
T=zeros(270,10);
for j=1:1:270
for i=1:1:10
switch i
case 1
T(j,i)=Tst0;
case 10
T(j,i)=(1-2*Fo*(1+B))*T0(i)+2*Fo*T0(i-1)+2*Fo*B*Tevn;
otherwise
T(j,i)=Fo*T0(i-1)+(1-2*Fo)*T0(i)+Fo*T0(i+1);
end
end
T0=T(j,:);
end
for j=1:1:270
Qab(j)=lamad2*(T(j,1)-T(j,2))/deltaZ*2*pi*R1*R1;
end
x=RealTime/10;
if (x==0)
Result=0;
else
Result=Qab(fix(x));
end
y=Result*RealTime;。