基于机理和数据的高炉智能建模方法

本发明涉及高炉冶炼，尤其涉及一种基于机理和数据的高炉智能建模方法。

背景技术：

1、高炉炉温理论研究已经有很多年的历史，其主要包括专家系统方法、理化机理方法、数理方法和机理数据结合方法。专家系统方法将专家经验直观化，将操作具体到细节，消除了由于人员不同而产生的操作差异，提高了高炉运行的稳定性，但是专家系统方法基于大量数据与知识，当知识库不足时，专家系统无法做出正确的判断；同时专家系统方法的结果多为表征意义，无法给出直观的温度场分布。机理方法详细的描述了高炉内部的多相流和相间传递，但是高炉内部机理复杂，耦合度高，因此一般的机理方法都做了很大的假设，比如说将密度、传热系数作为常数，完全忽略化学反应等，虽然借助机理在一定程度上保证了温度趋势上的正确性，但是由于假设性太大，其结果整体上偏差可能较大；且机理模型多为局部模型，只对高炉局部进行分析，忽略了高炉的整体性。大部分高炉炉温建立都是通过数据驱动的方法来建立输入和输出的关系，通过特征提取等方法来逼近输入和炉温或其他可以表征炉温的变量之间的定量关系，并通过这种关系来对下时间步的炉温进行预测和表征。该方法跳过了高炉机理的分析，基于高炉已有的历史数据对高炉未来数据进行预测，避开了机理方法中复杂的高炉机理。如果建立的关系足够精确，那么其输出对指导高炉调节和保证高炉稳定运行具有重要意义，但是当数据集的大小和数据不满足时，其准确度呈线性下降，即使数据集满足，当生产条件发生改变时，其准确度也会产生偏差，如果建立的关系本身没有问题，生产条件也符合，其结果也仅仅是通过其他数据来表征炉温，无法得到炉温的细节分布。

2、近些年来，随着计算机硬件和算法的发展，数学与物理相结合已成为一种主流的高炉温度场建立方法，但是目前的主流方法主要存在以下缺陷：对于边界条件设置均比较模糊且误差较大；无法表征详细的炉身温度场分布；当前大多方法均借用fluent等流体仿真软件采用离散法进行迭代计算，计算步骤复杂且计算量较大，计算慢，滞后性大。

技术实现思路

1、针对现有技术的不足，本发明提供一种基于机理和数据的高炉智能建模方法。

2、一种基于机理和数据的高炉智能建模方法，包括以下步骤：

3、步骤1：建立炉壁内表面温度场分布模型；

4、步骤1.1：获取炉壁温度概况；

5、在高炉炉壁的不同高度安置若干热电偶，同一高度的热电偶位于冷却管与高炉内壁之间沿着周向平均分布，通过热电偶的数据反馈，测得炉壁对应位置的温度，确定炉壁的物理参数以及炉壁内部的热电偶温度参数。

6、步骤1.2：建立高炉炉壁的导热微分方程并离散；

7、所述导热微分方程由固体能量守恒与傅里叶定律得到，有：

8、

9、其中λ为导热系数，t为温度，x、y、z为高炉的坐标方向。

10、将式(1)中的直角坐标系的导热微分方程经x＝rcosθ，y＝rsinθ，z＝h变换得到：

11、

12、其中r为高炉炉壁圆环横截面的半径，θ为高炉炉壁圆环横截面的角度，h为高炉的高度。

13、将高炉炉壁分成若干个横截面，形成若干个二维横截面导热模型，针对每个二维横截面导热模型，忽略纵向，得到：

14、

15、其中r0为二维横截面导热模型内环的半径，r2为二维横截面导热模型外环的半径。

16、采用有限差分法将偏微分方程离散成代数方程；由于目标区域规整，采用等距步长将目标区域划分为网格，使用网格上的代数值代替连续的数值域；将目标域在径向方向划分为m步，周向方向上划分为n步，整个目标域内包括内外边界共存在n*(m+1)数量的离散点，然后将式(3)进行中心差分离散，得到离散化后的导热微分方程如下：

17、

18、其中i、j分别为径向和周向离散点的下标，δr为二维横截面导热模型径向步长，δθ为二维横截面导热模型周向步长，ti,j为离散后的下标为i,j的网格点的温度，ri为下标为i的网格点的半径；

19、步骤1.3：确定高炉炉壁的边界条件；

20、包括炉壁内表面边界：

21、t1,j＝tin,j，r＝r0,0≤j≤n  (5)

22、耐火砖与冷却壁交界处边界：

23、

24、离散后：

25、

26、冷却壁外表面边界：

27、

28、离散后：

29、

30、式中，tin为内边界温度分布，即高炉炉壁内表面温度分布；tout为外边界温度分布，即高炉冷却壁与冷却管交接处的温度分布；tcw为冷却管中的冷却水的平均温度，tcw＝(tcw,in+tcw,out)/2，tcw,in为冷却水进口热电偶测温值，tcw,out为冷却水出口热电偶测温值；hx为冷却壁与冷却水的对流换热系数，v水为冷却水流速：

31、hx＝208+47.5v水  (10)

32、r0、r1、r2为高炉炉壁的物理参数，其中r0为二维横截面导热模型内环的内径，也就是高炉内壁的半径，r1为耐火砖与冷却壁交界处的半径，r2为二维横截面导热模型外环的半径，也就是冷却管位置处的半径；λ1为耐火砖的导热系数，λ2为冷却壁的导热系数。

33、步骤1.4：使用优化的pso算法求解炉壁内表面周向一维温度分布；

34、pso算法的基本公式为：

35、

36、式中，i为第i号粒子，i＝1,2,…,n；ω为惯性权重因子；vi(t)为第i号粒子迭代的第t代的速度；vi(t+1)为第i号粒子的第(t+1)代的速度；xi(t)为第i号粒子的第t代的位置；pi(t)为第t代时单个粒子的最佳位置；pg(t)为第t代时粒子群的全局最佳位置，c1、c2为加速度系数；r1、r2为[0,1]中均匀分布的随机数。

37、对pso方法中的惯性权重因子和加速度系数进行了重新定义，分别为

38、

39、式中ωini为ω的初值，ωfin为ω的最终值；m为整数，表示迭代次数，c1,ini和c1,fin分别表示c1的初始值和最终值；c2,ini和c2,fin分别表示c2的初始值和最终值，令ωini、c1,ini和c2,fin为较大值，ωfin、c1,fin和c2,ini为较小值。

40、粒子的差异性diversity(t)定义为：

41、

42、式中：lmax和lmin分别表示搜索空间的最大值和最小值，j表示每个粒子中的参数，i表示粒子编号。

43、根据差异性选择粒子下一时刻的搜索策略：

44、(1)diversity(t)＞ξ||r3＞a：

45、

46、(2)diversity(t)≤ξ&r3≤a：

47、

48、式中：r1、r2、r3、r4、r5表示[0,1]中均匀分布的随机数；ξ表示diversity(t)的限制值；a为控制粒子初始化概率的参数；

49、使用上述优化后的pso算法借助炉壁内部热电偶的真实值与计算值之差作为目标函数，对由步骤1.2得到的导热微分方程和步骤1.3得到的边界条件组成的高炉炉壁二维横截面导热模型进行求解，得到多层炉壁内表面周向一维温度分布；

50、步骤1.5：求解内表面一维纵向温度分布。

51、将高炉划分为炉喉和炉身上部，炉身下部，炉腰和炉腹3部分，然后利用三次曲线计算公式对纵向的每部分高炉的温度分布进行拟合：

52、tw(θj,h)＝β0,j+β1,jh+β2,jh2+β3,jh3，j＝1,2,...,n  (16)

53、由于每一部分的层数不一致，为了充分利用每一层的数据，采用式(17)进行拟合，得到高炉炉壁内表面不同角度的一维纵向温度分布。

54、

55、式中：h为高炉不同位置的高度，hi为第i层位置的高度，tw(θj,h)为高炉不同位置的温度，不同的角度、高度对应不同的内表面位置，β0,j、β1,j、β2,j、β3,j为三次曲线的待拟合系数，n为当前部分的层数。

56、至此，得到了炉壁内表面所有角度的一维纵向温度分布，组合在一起，即得到了炉壁内表面温度场分布。

57、步骤2：建立高炉炉身二维温度场分布模型；

58、步骤2.1：利用流体理论与n-s方程对高炉煤气进行描述，分析气相和固相的流动、传热和化学反应，通过建立描述高炉炉身压差、煤气速度、煤气温度和固体温度四个变量的偏微分方程建立高炉炉身二维温度场模型；将高炉炉身二维温度场模型根据耦合度划分为三个子模型：描述气体压力和速度的气体流动模型；描述气体温度的气体温度模型；描述固体温度的固体温度模型。

59、(1)所述气体流动模型由连续性方程与考虑了动量传递的动量守恒方程构成：

60、

61、

62、

63、式中：x、y分别为径向和轴向方向，u、v分别为径向和轴向速度，p为压力，t为时间，β为扩散项系数：μ为动力学黏性系数，ρ为气体密度，f1、f2为源项系数：

64、

65、式中ηg为气体粘度；εs为炉料的孔隙度；ψs为固体颗粒的形状因子；ds为固体颗粒直径；

66、(2)所述气体温度模型由能量守恒方程考虑能量传递和化学反应表示：

67、考虑能量传递之后的气体温度模型为：

68、

69、式中：tg为气体温度，ts为固体温度，s为源项对气体温度模型的影响，α为扩散项系数，k为傅里叶定律的导热系数，cv为比热，γ为传热系数：

70、

71、式中：ag-s为接触面积；kg为气体的导热系数；ds为平均粒度；re为雷诺数；prg为气体的对流传热普朗特数，代表间接还原反应对气体的能量影响，单位：w/m3；代表直接还原反应对气体的影响，单位：w/m3。

72、(3)所述固体温度模型由热传导方程考虑气固传热和化学反应表示：

73、

74、式中：代表间接还原反应对固体的能量影响，单位：w/m3；代表直接还原反应对固体的影响，单位：w/m3。

75、步骤2.2：确定高炉二维温度场模型的建模区域和边界条件，其中边界条件包括：

76、入口：风口为模型的入口边界；

77、出口：炉顶设置为压力出口，该边界的参数除了压力均从建模区域内部插值得到；

78、炉壁：壁面的流动设置为无滑移边界，同时，炉壁是无法穿透的，因此炉壁速度设置为0。传热设置为对流换热，即气体与炉壁进行对流换热；

79、对称轴：高炉为对称模型，对称轴处的垂直速度的梯度设为0，法向速度大小为0；

80、步骤3：求解高炉炉身二维温度场模型；

81、步骤3.1：对高炉炉身二维温度场分布模型的建模区域进行网格划分；

82、建模区域为鼓风口下部至炉喉十字测温之间，采用椭圆偏微分方程生成建模区域的网格，网格迭代是否收敛根据当前轮次与上一轮次的物理平面坐标的最大差值为依据，如下式所示：

83、

84、其中e为收敛误差，为当前轮次与上一轮次的物理平面坐标x轴的最大差值，为当前轮次与上一轮次的物理平面坐标y轴的最大差值。

85、通过求解泊松方程来生成bfc网格，设(x,y)为物理平面坐标系，(ξ,η)为计算平面坐标系，泊松方程进过变换，使用计算平面求解物理平面的方程为：

86、

87、式中：j、α、β、γ为系数：

88、

89、p(x,y)、q(x,y)为源项：

90、

91、其中边界上的φ、ψ的计算公式为

92、

93、求得边界上的φ、ψ后，采用线性拟合的办法求得内部点的φ、ψ值，再根据式(28)求得所有网格点的p、q，同时计算j、α、β、γ，得到物理平面和计算平面的关系。

94、步骤3.2：使用simple算法对高炉炉身二维温度场进行求解；

95、所述simple算法根据假设的压力分布借助动量方程求解出速度之后，使用压力修正方程进行迭代修正，使速度同时满足动量方程和连续性方程；

96、首先，修正前的速度由动量方程求得，以径向速度u进行说明，v与u同理。对动量方程离散有：

97、

98、式中为计算点修正前的速度项；为计算点相邻点修正前速度的耦合项；be为源项；为压力修正值梯度项。其中ae为计算点修正前的速度的系数，同理，anb和ae均为其后缀参数的系数，下标e/e、w/w、n/n和s/s表示的为主控制单元网格的左、右、下和上边界，p表示压力当前点；

99、修正后的速度和压力满足动量方程，有：

100、

101、其中ue′、unb′、pp′和pe′分别为计算点和相邻点的速度和压力修正值。

102、因此有：

103、

104、仅考虑压力修正值梯度项有：

105、aeue′＝ae(pp′-pe′)  (33)

106、因此有速度修正方程为

107、

108、式中：

109、

110、同时，修正后的压力和速度需要满足连续性方程，因此有：

111、(ue-uw)δy+(vn-vs)δx＝0  (36)

112、式中：δx为径向步长，δy为轴向步长。

113、将式(34)代入式(36)有压力修正方程为：

114、appp′＝aepe′+awpw′+anpn′+asps′+bp  (37)

115、式中：下标e、w、n和s表示计算点左、右、下和上方向的邻点，系数ap、ae、aw、an和as为：

116、

117、bp反映了靠近收敛的距离；

118、

119、步骤3.3：使用pod降阶算法加快高炉炉身二维温度场模型求解速度；

120、步骤3.3.1构建pod基

121、假设高炉炉身二维温度场温度模型时间总步数为t，空间在横向和纵向分别为n、m步，总节点个数为m；首先取l时间步数的simple算法求解过程的参数，包括径向速度u、轴向速度v、压力p、气体温度tg和固体温度ts分布，0＜l≤m,m＝mn，分别构成以下五个矩阵au＝(ui)m*l、av＝(vi)m*l、ap＝(pi)m*l、atg＝(tgi)m*l和ats＝(tsi)m*l，使用tg和ts代表气体温度和固体温度变量；

122、

123、

124、求取和的特征向量和特征值；

125、从特征向量构成的特征矩阵中取若干列向量构成pod基；

126、步骤3.3.2：定义并求解当前的pod特征值；

127、假设整个高炉炉身温度模型的解为：和

128、那么，它们收敛于

129、

130、定义pod特征值为：

131、

132、根据当前的前l时间段的数据以及对应的pod基求出前l时间段的pod特征值。

133、步骤3.3.3：求解下一时刻的pod特征值；根据高炉炉身二维温度场模型得：

134、

135、其中表示当前时间步和上述时间步的参数之差的集合，根据离散后的高炉炉身二维温度场模型得到。

136、根据式(40)和式(41)得到

137、

138、根据式(42)和式(43)有：

139、

140、其中

141、

142、根据式(44)求出下一时刻的pod特征值。

143、步骤3.3.4：求取下一时刻的参数，包括轴向速度、纵向速度、压力、气体温度和固体温度分布直至收敛；

144、根据式(41)得到式(46)，结合下一时刻的pod特征值以及pod基，得到下一时刻的参数分布

145、

146、然后继续求解直到收敛便得到高炉内部温度模型的目标值

147、步骤3.3.5：pod基更新；

148、有pod误差为

149、

150、得

151、

152、根据式(40)，有

153、

154、根据微分收敛界定有：

155、在保证离散方程稳定的前提下，根据式(48)和式(49)得

156、

157、其中和分别表示时间步为n和n-1，径向步长为i，轴向步长为j处的径向速度；和分别表示时间步为n和n-1，径向步长为i，轴向步长为j处的径向速度；δt为时间步长；δx为径向步长；δy为轴向步长，en+1和en分别为时间步为n+1和n的误差。

158、在保证方程稳定的条件性条件下，高炉二维温度分布模型的精度解与pod降阶算法的数值解之间存在如下误差估计：

159、

160、其中

161、

162、如果满足

163、

164、则δ≤1。

165、如果选择的mu、mv、mp、mtg和mts满足那么cn(δ)＝(1+δ)n-l便是确定pod基更新的决定性因素。因此有如下结论：

166、如果便需要使用新的快照矩阵来确定新的pod基来代替旧的pod基。

167、步骤4：根据需要进行实验的高炉的物理参数与理化参数，确定高炉炉壁内部热电偶情况，得到对应高炉的炉壁内表面温度分布和炉身二维温度分布。

168、采用上述技术方案所产生的有益效果在于：

169、(1)对高炉炉壁进行机理与数据分析，利用传热模型结合热电偶以及改进的粒子群算法计算高炉内壁温度，提高模型边界条件设置的精确度，进而提高模型的精确度。风口回旋区的燃烧反应复杂，采用红外摄像仪从风口小镜抓取图片测量风口回旋区温度，以此温度代替鼓风温度，使入口边界条件更加精确的同时，可以忽略风口回旋区内燃烧反应对温度的影响，简化计算，提高精确度。

170、(2)采用单相流仿真思路，使用单相流对高炉进行分析，并分析单相流直接的质量、动量和能量传递，以较低的计算成本更好的贴近现场的实际需要，提高实时性。

171、(3)在传统模型迭代求解的过程中，结合pod算法对模型进行降阶，加快模型的计算速度，减少计算成本，提高实时性。

172、总体来说，本发明从机理出发，结合图像、数据和数学方法，得到了一种在保证一定计算精度的前提下，提高计算速度的炉身二维温度场计算方法。