反应炉的温度分布控制状态空间法有限元模型

反应炉的温度分布控制状态空间法有限元模型

反应炉的操作就如高炉一样仍然依靠熟练操作人员对于炉内部现象和高温的复杂性的经验和直觉。在熟练操作人员减少和技术传承困难的情况下对于稳定操作有着极大的需求。本文旨在建立一个反应炉控制的数学描述。进一步，提出了反应炉的线性二次高斯控制系统，在这里使用炉壁附近测量的数据来估计炉内温度分布。炉温度分布的控制是基于估计的炉内温度分布边界条件的变化。并通过数值实验来检验提出的控制方法的可行性。

1．介绍

反应炉就如高炉一样在钢铁工业中发挥着重要作用。炉的性能伴随着炉设备的扩大和生产的增加已经有了明显的改进。但是，炉的操作仍然依靠熟练操作人员对于炉内部现象和高温的复杂性的经验和直觉。有许多主要因素，例如气体流量，化学反应，烧穿点和铁矿石的吻合和炉中填料的运动。

由于炉调整的困难，高炉的平稳操作就十分必要。最近，伴随着熟练操作人员减少和技术传承困难，出现了对于自动控制的极大需求。

在本研究中，我们研究反应炉控制系统的设计来估计和控制反应炉中的温度分布。正如前面提到的，反应炉有许多复杂的现象，如化学反应、炉内高压和高温。只有靠近炉壁的变量可以被测量和用来进行炉控制。因为一个反应炉体积大，并且炉内的环境不能直接改变。因而控制系统通过测量炉壁附近的数据来估计炉内温度分布，并基于估计的炉内温度分布边界条件的改变来建立反应炉温度分布的控制。

首先，建立反应炉的仿真，来对应于炉内气流和温度分布。数值模拟的炉只有在炉壁

附近有着测量仪器，并只能在边界操作控制输入。这项研究用反应炉仿真进行。为了进行控制系统设计，状态空间模型是从运用有限元方法（FEM）到反应炉仿真推导的。在此过程中，线性系统理论被运用在控制系统的设计中。进一步，采用线性二次高斯（LQC）控制来处理问题描述。接下来的部分将描述具体内容。

2．反应炉模型

高炉被用在将铁矿石变为生铁的过程中。图1是高炉的示意图，高炉的高度大约是40米，直径大约是20米。铁矿石和焦炭从高炉的上部交替进入。此外，有多个底部的高炉风口，通过这些风口吹约1200摄氏度的热空气。通过吹热空气，焦炭燃烧并产生一氧化碳。铁矿石在大约8小时中被分解为生铁，并且一氧化碳和生铁在高炉的底部堆积。另一方面，热空气则从高炉顶部出去。

就如我们前面提到的，高炉有许多复杂的现象，如化学反应、炉内高压和高温。然而，高炉的测量仪表只能安装在炉壁附近。

接下来，将会建立反应炉模型。

二维反应炉模型见图2。Uin1、Uin2、Vin1、Vin2是风口吹气的速度。节点号(i,j)如图2所示分配给节点，反应炉模型的特性如表1所示。

测量仪表的设立如图3，其中△是气体流量计，□是晶体温度计，这些仪表设置在反应炉模型的外面和固体层的顶部。

假设压力分布为Px,y,t，则气流分布为：

Vx,y,tUx,y,tiVx,y,tj （1）

Tgx,y,t炉内气相温度分布为

，炉内固相温度分布为Tsx,y,t。如表2所示，这些变量

通过节点号(i,j)和时间n来建立。

反应炉模型通过运用MAC法和有限差分法[1]来计算表2中的离散变量来建立接下来部分关于P、V、

Tg和Ts的控制方程

图1 高炉示意图 图2 二维反应炉模型 图3 测量仪表

表1 反应炉模型的特性

表2 节点号(i,j)在n时的变量

2.1 气流模型[2,3,4,5]

假设压力和气体流量通过连续性方程、Navier-Stokes方程和Erugan方程[5,6]来描

述，我们得到：

UV0 xy （2）

U tUUxVUyPx1Re2U2Ux2y2f221f2UVU VVVP12V2 tUxVyyReV22x2y2f1f2UVV 边界条件：

 炉底部和炉壁：

UV0,Pn0 n：法向量  出口：

Uni,21Uni,20,Vni,21Vni,20,

Pni,210i6,7,8  进口：

UnUn1,3in1,U13,3Uin2,

Vnn4,1Vin1,V10,1Vin2 3）

4）

（5）

（6）

（7）

（ （

Uin1其中P是压力分布，V是气流分布，Re是雷诺数，f1和f2是Erugan方程的系数，、Uin2、Vin1、Vin2是风口吹气的速度。

2.2 气相温度模型[2,3,4,5]

假设气相温度是通过气固两相间热量传递的能量平衡方程来描述，我们得到：

Tg

2Tg2Tg1UV2txyRePrxy2TgTgTgTs （8）

边界条件：

 炉底部：

Tg y0 （9）

 炉壁和出口：

Tg naTgTout n：法向量 （10）

 进口：

Tg1,3Tgin1,Tg13,3Tgin2,nnn

Tg4,1Tgin3,Tg10,1Tgin4n （11）

其中

TgToutTs和a是热传递系数，是气相温度分布，是固相温度分布，Pr是普朗特数，

是外部温度，

Tgin1、

Tgin2、

Tgin3和

Tgin4是风口吹气的温度。

2.3 固相温度模型[2,3,4,5]

假设固相温度是通过材料的反应热和气相间的热传递的热量传导平衡来描述[8]，我们得到：

2Ts2TsTsx2y2t TsTgf3U,V,Tg （12）

边界条件：

 炉底部：

Ts0 y （13）

 固相层的顶部：

Tg ybupTsTg （14）

 炉壁：

TsbTsToutn n：法向量 （15）

其中，，，

bup和b是热量传递系数，f3是材料的反应热，f3的方程为：

expE/RTgQ1f3U,V,TgQ2U2V2Dp1expTg1000/K1 （16）

其中

Dp是材料颗粒直径，Q1，Q2，E，R和K1是常数。

3．状态空间模型的求导

本章会基于方程（2）-（16）来建立状态空间模型。状态空间模型表示领域分布目标以及包括的状态方程和输出方程的特点。如果对状态空间模型进行求导，就可以使用线性控制理论[9,10,11]。

首先，通过方程（2）-（16）和需要的分布对线性化方程进行求导。然后，广义系统通过运用FEM微分来得到线性化方程。下一步，状态方程通过运用广义逆矩阵来求导得出广义系统。最后输出方程是通过图3中显示位置的测量值来建立。通过这个过程，状态空间模型就建立了。

3.1邻域分布目标的线性化[12]

方程（2）-（16）在附录A的情况下的一个稳态解见图4-7。

假设图4-7所示的稳态解所需的分布（Pr，Vr，

Tgr，Tsr）和风口吹气的稳态速率所需

~~~~TUUVV的分布（in1，in2，in1，in2）。定义摄动变量（P，U，V，g，Ts）和（Uin1，

Uin2，Vin1，Vin2）：

PPrP,UUrU,VVrV

TgTgrTg,TsTsrTs （17）

~~Uin1Uin1Uin1,Uin2Uin2Uin2~~VVV,VVin1in2in2Vin2 in1in1 （18）

将其代入方程（2）-（16），我们得到了线性化方程：

U xVy0 UtUrUUrUxUUrxyVVryPx12U2U22rVrUrVrRex22Uy2f1Uf2U22Uf22V2VrVrUrr VVrVVrVtxUUrxyVVryPy12V2VU2fr2V2rUrVrRex2y21Vf2U22Vf222UrVrUrVr 边界条件：

 炉底：

UV0,Pn0 n：法向量  出口：

（19）20）

21）

（22）

（ （

Uin,21Uin,20,Vi,n21Vi,n20,

Pi,n210i6,7,8 （23）

 进口：

Unn1,3Uin1,U13,3Uin2,

Vnn4,1Vin1,V10,1Vin2 TgTgrTgTgrTgtxUUrxyVVry1RePr2Tg2TgxT

2y2gTs 边界条件：

 炉底部：

Tg

y0  炉壁和出口：

Tg naTg n：法向量  进口：

（24）

25）

（26）

（27）

（

nnnTg1,3Tgin1,Tg13,3Tgin2,

Tg4,1Tgin3,Tg10,1Tgin4n （28）

2Ts2TsTs22TsTgt

边界条件：

 炉底部：

 固相层的顶部：

 炉壁：

xyf3f3UTUf3UgTVTgTVTTggTVgrUVrUgrUgTrVVUrrVgrUVrr Ts y0 Tg

ybupTsTg Ts nbTs 29）

（30） （31）

（32） （

图4 压力的稳态分布（t=5[-]） 图5 气流的稳态分布（t=6[-]）

图6 气相温度的稳态分布（t=90[-]） 图7 固相温度的稳态分布（t=80[-]）

3.2 有限元分析在线性方程组中的应用

本章，FEM将会应用到线性方程组中。摄动变量的定义位置和有限元分析的三角元素见图8-10，其中○是摄动变量的定义位置，其中，U、V和

Tg的摄动变量定义在相同

的位置。

图8-10中典型的三角元素e见图11，其中4,5,6是边线的中点。通过式（33），摄动变量定义在图11所示的边的节点和中点。另一方面，P的摄动变量和式（34）中所需的值定义在节点。

eUtTUe1Ue2UUe3Ue45Ue6eVtTVe1V2Ve3Ve4Ve5Ve6TeeeTgtTg1Teg2Teg3Teg4Tg5Tg6

eTstTTs1Tes2Tes3Tes4TeTes5s6 TeePtPe1P2P3eeeUrtTUr1Ur2Ur3VrtTVer1Ver2Ver3eTgrtTTegr1Tegr2Tgr3 eeeTsrtTTsr1Tsr2Tsr3 U、V、Tg、Ts、P、Ur、Vr、Tgr、Tsr通过e的逼近为：UeMUt,VeMVt

TegMTgt,TesMTst PeLPtUerLeUr,VrLVr

TegrLTgr,TesrLTsr 其中，M和L是：

33）

（34）

（35）

（36）

（

M2L11L12L21L22L31L3LL1L2L3L1A3A1A2,L,L23eeey1x11,A11A11xy222221x3y311x1y111,e11A31xy2222y1x1x1xx3x1x2x3y1yy3y1y2y34L1L24L2L34L3L1

（37）

从式（19）-（32），（35），（36），通过选择一组测试方程M和L及以下的要求可以得到逼近

UeVeeLxy

Tdxdy0 （38）

eeeeeeeUUUUPeeUTeerrMtxUUrxyVVryx12Ue2Ue2Rexy2f2UerVerUerVer22e2UerVereefUfU1222UerVer22Vdxdy0

（39）

eeeeVeVreVrPeVeVeeMtxUUrxyVVryxTe12Ve2Ve2Rexy2f2UerVerUerVer22eUer2VereVef1Vf222UerVer22Udxdy0

（40）

Me

TgTgrTgTgTgrUUVVrrtxxyy2Tg2Tg1TgTsdxdy022RePryxT （41）

Ts2Ts2TseMx2y2tTTsTg

f3f3f3UVTgdxdy0TgTgrTgTgrTgTgrUUVTgUUUrUUrrVVrVVrVVr （42）

其中式（38），式（39）和（40）第7部分，式（41）第6部分，式（42）第2部分用了部分积分法。式（42）的第4-6部分近似为：

Tf3UTgTgrUUrVVrfL3UfL3VeTgTgr1eUUr1eVVr1f3UUUre2VVre2TgTgre2f3eTT3eUUgUgrr3eVVr3f3eTT3eVUgUgrr3eVVr3Tf3VTgTgrUUrVVreUUr1eVVr1eTgTgr1f3VUUre2VVre2TgTgre2f3TgTgTgrUUrVVrfL3TgeTgTgr1eUUr1eVVr1f3TgUUre2VVre2TgTgre2f3eTgTgTgre3UUr3eVVr3

T （43）

根据式（38）-（43），考虑图8-10中的所有三角元素和边界条件，得到：

Ui0i:炉壁号Vi0i:炉壁号Pi0i2U26Uin1,U30Uin2,V32Vin1,V34Vin2,Tg260,Tg300,Tg320,Tg340

（44）

假设该状态向量x：（99×1）和控制向量u：（4×1）为：

xx1x1x11x2991Tx2P1x21x11Tg1Tg25x31U3x41V3uUin1U7V7Uin2U8V8P3P12111Tg27Tg28Tg29Tg31Tg33x31x41881Tg35311x21Ts1Ts25251U9U27V9V27Vin2T41U28V28V29161U29161Vin1

（45）

我们得到：

AxBu:E9999,A9999,B998Ex0xEx1,E1100x2A11A12B1A,BBAA22212

（46）

其中A220，由于E不是非奇异矩阵，式（46）不称为控制系统设计中广泛使用的状态方程。所以方程（46）称为广义系统[14]。因为在式（19）-（21）中P没有时间分化，所以建立了广义系统。

如果式（46）的广义系统脉冲可控且有限动态稳定，就可能设计动态系统[14]。但是式（46）不是脉冲可控。所以，在下面的章节，我们会将式（46）的广义系统通过广义逆矩阵 转化为状态方程[15]，从而能设计控制系统。

图8 P的摄动变量 图9 U的摄动变量

图10 Ts的摄动变量 图11 一个典型的三角元素

3.3 广义逆矩阵在广义系统中的应用[15]

从式（46），我们得到

A22x2A21x1B2u 因为A#220，广义逆矩阵A22为0，所以从式（46）和（47）我们知道：

47） （

x2A22A21x1B2u0 （48）

#然后代入式（46），我们得到状态方程：

x1E11A11x1E11B1u （49）

11本文假设E11A11为A，E11B1为B，考虑了图3所示的仪表的位置，我们得到状态空间模型：

AxBu:A8888,B888xyCx:C1388xx1131T11uUin1Uin2Vin1Ts10Vin2T41Ts11Ts15Ts21Ts25U3V3131（45）

TyTs1Ts6因为式（50）是稳定和可观测的，可以使用线性控制理论。

4. 反应炉控制系统

本文采用线性二次高斯控制（LQG）来处理下面描述的问题[10]。对于LQG控制的详细解释见附录B。

反应炉模型的LQG控制系统见图12。LQG控制系统通过式（50）来创建。其中K是

ˆ是反应炉模型估计的内部状态变量。下面是控制流优化调节增益，L是卡尔曼滤波增益，x程：

第一步：通过卡尔曼滤波和图3中仪表测量的数据来估计反应炉的内部状态变量。

第二步：基于优化调节和第一步估计的内部状态变量，修改风口吹气的速度来匹配反应炉模型所需的初始分配。

图12 LQG优化控制系统

Uin2Vin1Vin2Vx,y,0、g其中风口吹气速度为（Uin1、、、），初始分布为（Px,y,0、Tsx,y,0Tx,y,0、

），需要的分布为（Prx,y、Vrx,y、

Tgrx,y、Tsrx,y）和（Pr、Vr、

Tgr、Tsr）

是式（2）-（16）的一个稳态解。

5. 数值实验

在本研究中，反应炉模型的温度分布控制采用第2章中建立的反应炉模型。本章，将会提供LQG控制的数值结果。下面是问题的描述：

Uin2Vin1Vin2 确定了风口吹气速度（Uin1、、、），因而在图13-16初始的分布（Px,y,0、

Vx,y,0、Tgx,y,0、Tsx,y,0）能被控制到图4-7中所需的分布。

图15和16中初始的温度分布总的高于图6和图7中所需的。

假设附录B中的Q、W、R、V为：

56elements32elementsQWdiag10,,100.1,,0.1:8888RV7.0106diag1,1,1,1

因为本研究的目的是估计反应炉模型的内部温度，并和反应炉模型需要的内部温度匹配，因此定义了两个目标函数：

ˆtf

Tgiti1i26,30,32,3435ˆtTtTˆtTgitTgisisii125 （51）

ˆtTgi其中

和Tsit是反应炉模型在图9和10中节点的温度，和Tsit是在上述点通

ˆ过LQG控制得到的估计温度。方程（51）代表上述点的温度的估计误差。

ftTgi,jTgi,jTsi,jTsi,jnrnri1j1i1j113211313 （52）

其中

Tgi,j

r

和

Tsi,j

r

是在节点（i,j）所需的温度，方程（52）代表所有节点温度误差。

从图17和18可以看到风口吹风速度的改变。两个目标函数的变化见图19。从图17和18可以看出，风口吹气速度在接近15[-]的时候操作。当风口吹气，两个目标函数减少。优化的温度分布见图20和21。所以反应炉模型的内部温度变量就被估计出，且温度分布接近了通过LQG控制系统得到的所需的值。

图13 压力的初始分布

图14 气流的初始分布

图15 气相温度的初始分布 图16 物料温度的初始分布

图17 Uin1和Vin1

图18 Uin2和Vin2图19 fˆt和ft的变化

图20 气相温度的优化分布（t=50[-]） 图21 材料温度的优化分布（t=50[-]）

6. 总结

本文建立了反应炉的仿真，能够用来计算炉内气相温度、固相温度、气体流量和压力的分布。该仿真的测量值均来自炉壁附近和边界的控制输入。为了设计控制系统，状态空间模型通过有限元分析来推导。通过该过程，线性控制理论能够被用在控制系统的设计上。接下来采用LQG方法来处理描述的问题。我们通过LQG控制系统来对仿真的温度分布进

行估计和控制。

将来，我们需要对于估计的内部温度进行改进来调节温度分布在很短的时间内到需要的值。为了这个目标，我们将会采用其他的线性控制理论如H2控制和H∞控制。

参考文献

[1] T. Kawamura: Fluid analysis I, Asakura shoten(1996), 1-83.

[2] T. Shibuta: Master Thesis. Okayama Univ, (2003)

[3] K. Ishimaru: Master Thesis. Okayama Univ. (2005)

[4] T. Takeda: Graduation Thesis. Okayama Univ. (2006)

[5] K. Takatani, T. Inada and K. Takata: ISIJ International, 41-10 (2001), 1139-1145.

[6] I. Imai: Fluid Dynamics, Shokabo (1973), 259-284.

[7] S. Ergun: Chemical Engineering Progress, 48-2(1952), 89-94.

[8] M. Iri and Y. Iri: Partial Differential Equation, Asakura shoten (1983), 13-35.

[9] K. Z. Liu: Linear Robust Control, Coronasha(2002), 78-179.

[10] S. Hosoe: System and Control, Ohmsha (1997), 97-110

[11] A. Fujimori: Robust Control, Coronasha (2001), 42-139.

[12] O. M. Aamo and M. Krstic: Flow Control by Feedback, Springer (2003), 14-17.

[13] G. Yagawa: Beginning FEM of Flow and Heat Transfer, Baifukan (1983), 86-170.

[14] T. Katayama: Optimal Control of Linear System, Kindai Kagaku sha (1999), 125-180.

[15] H. Kimura: Linear Algebra, University of Tokyo Press (2003), 112-114.

附录A

在本研究中，采用有限差分方法来计算反应炉模型的气流和温度分布，图22是计算的流程图。气流的计算是采用MAC方法。这些计算会一直持续，直到满足收敛性条件。收敛性条件的定义为：

fpmaxPx,y,npPx,y,np10.0001fUmaxUx,y,nVUx,y,nV10.0001fVmaxVx,y,nVVx,y,nV10.0001np1,2fTgmaxTgx,y,10nTgTgx,y,10nTg10.1nV1,2nV1,2nTg1,21,2

fTsmaxTsx,y,10nTsTsx,y,10nTs10.1nTs （53）

这些条件见图23-25，表3和4来得到3.1章中图4-7的稳态解。反应炉模型在表3的条件下在dt（=0.01）的间隔时计算了表2中的离散值。

图22 数学模型计算的流程图 图23 所需分布的f1的条件

图24 所需分布的f2的条件 图25 所需分布的Dp的条件

表3 所需分布的条件1 表4 所需分布的f2的条件2

附录B

在LQG控制中，状态空间模型假设受到干扰：

AxBuw:Ann,Bmnx yCxv:Cln （54）

其中w是干扰，v是测量噪声。w和v是假设为零均值和协方差的不相关高斯随机过程。

Ewt0,Evt0wtTEwvt vTW0t0V （55）

其中W和V是常数阵，E{·}是期望运算符，δ是δ函数，LQG控制问题是为了找到反馈u来是成本功能最低。

T1TTJElimxQxuRudtTTt0 （56）

TTQQ0,RR0。该问题的唯一解为： 其中Q和R是恒定的加权矩阵满足

ˆ,KR1BTPuKxˆAxˆBuLyCxˆ,xT1ˆt00 （57） LSCV,x

其中P和S唯一解决代数黎卡提方程。

ATPPAQPBR1BTP0,TT1 ASSAWSCVCS0 （58）

如果满足以下两个条件，我们可以说明存在一个可接纳的LQG控制系统：

TT是稳定的，即： A,B,A,C（A.1）

rankSIABn,rankSIATCTn,Res0,sC （59）

（A.2）Q,A,W,AT是可观的，即： rankSIASIATQn,rankn,Res0,WsC 60） （