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

1、反应炉的温度分布控制状态空间法有限元模型反应炉的操作就如高炉一样仍然依靠熟练操作人员对于炉内部现象和高温的复杂性的经验和直觉。在熟练操作人员减少和技术传承困难的情况下对于稳定操作有着极大的需求。本文旨在建立一个反应炉控制的数学描述。进一步，提出了反应炉的线性二次高斯控制系统，在这里使用炉壁附近测量的数据来估计炉内温度分布。炉温度分布的控制是基于估计的炉内温度分布边界条件的变化。并通过数值实验来检验提出的控制方法的可行性。1介绍反应炉就如高炉一样在钢铁工业中发挥着重要作用。炉的性能伴随着炉设备的扩大和生产的增加已经有了明显的改进。但是，炉的操作仍然依靠熟练操作人员对于炉内部现象和高温的复杂性的经

2、验和直觉。有许多主要因素，例如气体流量，化学反应，烧穿点和铁矿石的吻合和炉中填料的运动。由于炉调整的困难，高炉的平稳操作就十分必要。最近，伴随着熟练操作人员减少和技术传承困难，出现了对于自动控制的极大需求。在本研究中，我们研究反应炉控制系统的设计来估计和控制反应炉中的温度分布。正如前面提到的，反应炉有许多复杂的现象，如化学反应、炉内高压和高温。只有靠近炉壁的变量可以被测量和用来进行炉控制。因为一个反应炉体积大，并且炉内的环境不能直接改变。因而控制系统通过测量炉壁附近的数据来估计炉内温度分布，并基于估计的炉内温度分布边界条件的改变来建立反应炉温度分布的控制。首先，建立反应炉的仿真，来对应于炉内气

3、流和温度分布。数值模拟的炉只有在炉壁附近有着测量仪器，并只能在边界操作控制输入。这项研究用反应炉仿真进行。为了进行控制系统设计，状态空间模型是从运用有限元方法（FEM ）到反应炉仿真推导的。在此过程中，线性系统理论被运用在控制系统的设计中。进一步，采用线性二次高斯（LQC ）控制来处理问题描述。接下来的部分将描述具体内容。2反应炉模型高炉被用在将铁矿石变为生铁的过程中。图1是高炉的示意图，高炉的高度大约是40米，直径大约是20米。铁矿石和焦炭从高炉的上部交替进入。此外，有多个底部的高炉风口，通过这些风口吹约1200摄氏度的热空气。通过吹热空气，焦炭燃烧并产生一氧化碳。铁矿石在大约8小时中被分解

4、为生铁，并且一氧化碳和生铁在高炉的底部堆积。另一方面，热空气则从高炉顶部出去。就如我们前面提到的，高炉有许多复杂的现象，如化学反应、炉内高压和高温。然而，高炉的测量仪表只能安装在炉壁附近。接下来，将会建立反应炉模型。二维反应炉模型见图2。U in 1、U in 2、V in 1、V in 2是风口吹气的速度。节点号(i , j 如图2所示分配给节点，反应炉模型的特性如表1所示。测量仪表的设立如图3，其中是气体流量计，是晶体温度计，这些仪表设置在反应炉模型的外面和固体层的顶部。假设压力分布为P (x , y , t ，则气流分布为：V (x , y , t =U (x , y , t i +V

5、(x , y , t j （1） 炉内气相温度分布为T g (x , y , t ，炉内固相温度分布为T s (x , y , t 。如表2所示，这些变量通过节点号(i , j 和时间n 来建立。反应炉模型通过运用MAC 法和有限差分法1来计算表2中的离散变量来建立接下来部分关于P 、V 、T g 和T s 的控制方程 图1 高炉示意图 图2 二维反应炉模型 图3 测量仪表表1 反应炉模型的特性 表2 节点号(i , j 在n 时的变量 2.1 气流模型2,3,4,5假设压力和气体流量通过连续性方程、Navier-Stokes 方程和Erugan 方程5,6来描述，我们得到：U V+=0 （2

6、） x yU U U P 12U 2U 22 =- U +V -+-f +f +V U （3）1222 x x Re x t y y V V V P 12V 2V 22 =- U +V -+-f +f +V V （4）1222 t y y Re x y x边界条件： 炉底部和炉壁：U =V =0, 出口： 进口：nU 1n , 3=U in 1, U 13, 3=U in 2,(P=0 n ：法向量 （5） nU i n , 21=U i n , 20, V i , n 21=V i , n 20, Pn i , 21=0(i =6, 7, 8（6）Vn4, 1=V in 1, Vn 10,

7、1=V in 2（7）其中P 是压力分布，V 是气流分布，Re 是雷诺数，f 1和f 2是Erugan 方程的系数，U in 2、U in 1、V in 1、V in 2是风口吹气的速度。2.2 气相温度模型2,3,4,5假设气相温度是通过气固两相间热量传递的能量平衡方程来描述，我们得到：2T g 2T g 1 =-U -V +2 t x y Re Pr x y 2T gT gT g边界条件： 炉底部： 炉壁和出口： 进口：-(T g -T s （8） T g y=0 （9）T g n=a (T g -T out n ：法向量 （10）T g 1, 3=T g in 1, T g 13, 3=

8、T g in 2, T g 4, 1=T g in 3, T g 10, 1=T g in 4nnn n（11）其中T g 是气相温度分布，T s 是固相温度分布，Pr 是普朗特数，和a 是热传递系数，T out 是外部温度，T gin 1、T gin 2、T gin 3和T gin 4是风口吹气的温度。2.3 固相温度模型2,3,4,5假设固相温度是通过材料的反应热和气相间的热传递的热量传导平衡来描述8，我们得到：2T s 2T s T s= +22 t y x边界条件： 炉底部： 固相层的顶部： -(T s -T g +f 3(U , V , T g （12）T s=0 （13）y炉壁：T

9、 g y=-b up (T s -T g （14）T s=b (T s -T out n ：法向量 （15） n其中，b up 和b 是热量传递系数，f 3是材料的反应热，f 3的方程为：exp -E /(R T g Q 1f 3(U , V , T g = Q 22+V 2 （16）D p 1+exp T g -1000/K 1其中D p 是材料颗粒直径，Q 1，Q 2，E ，R 和K 1是常数。3状态空间模型的求导本章会基于方程（2）-（16）来建立状态空间模型。状态空间模型表示领域分布目标以及包括的状态方程和输出方程的特点。如果对状态空间模型进行求导，就可以使用线性控制理论9,10,11

10、。首先，通过方程（2）-（16）和需要的分布对线性化方程进行求导。然后，广义系统通过运用FEM 微分来得到线性化方程。下一步，状态方程通过运用广义逆矩阵来求导得出广义系统。最后输出方程是通过图3中显示位置的测量值来建立。通过这个过程，状态空间模型就建立了。3.1邻域分布目标的线性化12方程（2）-（16）在附录A 的情况下的一个稳态解见图4-7。假设图4-7所示的稳态解所需的分布（P r ，V r ，T gr ，T sr ）和风口吹气的稳态速率所需的分布（U in 1，U in 2，V in 1，V in 2）。定义摄动变量（P ，U ，V ，T g ，T s ）和（U in 1，U in 2

11、，V in 1，V in 2）：P =P r +P , U =U r +U , V =V r +V T g =T g r +T g , T s =T s r +T s（17）U in 1=U in 1+U in 1, U in 2=U in 2+U in 2（18） V in 1=V in 1+V in 1, V in 2=V in 2+V in 2将其代入方程（2）-（16），我们得到了线性化方程：U V+=0 （19）x yU r U U U r U P=- U +U +V +V -r r t x y y x x2222 （20）+1Re U U 2U r +V r U r V r x 2

12、+y 2-f 1U -f 222U -f 222V r +V r r +V rV V r V V r V t =- xU +U r x +y V +V r y -Py 12+ Re V x 2+2Vy 2-f 1V -f U 2r +2V 2r U r V r 222V -f 222U r +V r r +V r 边界条件： 炉底：U =V =0, Pn=0 n ：法向量 出口：U n n i , 21=U n i , 20, V n i , 21=V i , 20, Pn i , 21=0(i =6, 7, 8 进口：n U n1, 3=U in 1, U 13, 3=U in 2,Vnn

13、4, 1=V in 1, V10, 1=V in 2T g =- T g rT g T g r T g tx U +U r x +y V +V ry 22+1 T g T g Re Pr 2+-(T x y 2g -T s 边界条件： 炉底部：T g y=0 炉壁和出口：21）22） 23）24）25）26）（ （ （ （ （ 进口：T g n=a T g n ：法向量 （27）T g 1, 3=T g in 1, T g 13, 3=T g in 2, n n（28）T nng 4, 1=T g in 3, T g 10, 1=T g in 4T s2T 2s T s t = 2+xy 2-

14、(T s -T g +f3f 3UT T g =T grV +f 3T T gU g =T U +U T V =U grV V rV =U V rg U g =T V =V U grrrr r边界条件： 炉底部：T sy=0 固相层的顶部：T g y=-b up (T s -T g 炉壁：T sn=b T s 图4 压力的稳态分布（t=5-） 图5 气流的稳态分布（t=6-）29）30）31）32） （ （ （ （ 图6 气相温度的稳态分布（t=90-） 图7 固相温度的稳态分布（t=80-）3.2 有限元分析在线性方程组中的应用本章，FEM 将会应用到线性方程组中。摄动变量的定义位置和有限元

15、分析的三角元素见图8-10，其中是摄动变量的定义位置，其中，U 、V 和T g 的摄动变量定义在相同的位置。图8-10中典型的三角元素e 见图11，其中4,5,6是边线的中点。通过式（33），摄动变量定义在图11所示的边的节点和中点。另一方面，P 的摄动变量和式（34）中所需的值定义在节点。U (t TV (t TT (t T T (t Tg s=V=U 1ee 1ee U 2e U 3e U 4e U 5e U 6=TV 2e T g 2T s 2e eV 3e T g 3T s 3e eV 4e T g 4T s 4e eV 5e T g 5T s 5e eV 6eee=T g 1e s

16、1T g 6T s 6（33）P (t T =P 1e P 2e P 3e U (t T =U r 1e U r e 2U r 3e e e eV (t T =V r 1V r 2V r 3 （34） T (t T =T gr 1e T gr e 2T gr 3e T (t T =T sr 1e T sr e 2T sr 3e r r gr srU 、V 、T g 、T s 、P 、U r 、V r 、T gr 、T sr 通过e 的逼近为：U e =M U (t , V e =M V (t T g =M T g (t , T s =M T s (t ee（35）P e =L P (t eU

17、r =L U r , V r =L V reerT g r =L T g , T s r =L T s re（36）4L 2L 34L 3L 1其中，M 和L 是：M =(2L 1-1L 1(2L 2-1L 2(2L 3-1L 3L =L 1L 2L 3L 1=4L 1L 2A 1, e 11A 1=12111A 3=121L 2=x x 2x 3x 1x 2xA 3A 2, L =3e e y 11y 2, A =122y 31y 111ey 2, =12y 1（37）x 1x x 3x 1x 2x 3y 1y y 3y 1y 2y 3从式（19）-（32），（35），（36），通过选择一组

18、测试方程M 和L 及以下的要求可以得到逼近U e V e e L x +y dxdy =0 （38）Tee e e eU e U r e U r P e U e U e eM t +x U +U r x +y V +V r y +xT12U e 2U e-+2Re x y 2+f 2U e r V e r e r +V e r22e2U e r +V e r eU e +f 1U +f 222e r +V e r22（39）V dxdy =0ee e e eV e V r e V r P e V e V e eM t +x U +U r x +y V +V r y +xT12V e 2V e-

19、+2Re x y 2+f 2U e r V e r e r +V e rTU e r +2V e r eV e +f 1V +f 222U e r +V e r22（40）22U dxdy =0eT g T g r T g T g T g r+U +U r +V +V re M t x x y y T g T g 1 -+(T g -T s dxdy =022Re Pr y x 22（41）T s 2T s 2T s- e M x 2+y 2 tT+(T s -T g （42） f f f-3T =T U -3T =T V -3T g dxdy =0g gr g gr T =T U U =U

20、r V U =U r T g g gr U =U r V =V r V =V rV =V r 其中式（38），式（39）和（40）第7部分，式（41）第6部分，式（42）第2部分用了部分积分法。式（42）的第4-6部分近似为：f 3UT g =T gr =U r V =V rf =L 3U f =L 3V e=U r 1eV =V r 1eT g =T gr 1f 3UU =U r e 2V =V r e 2T g =T gr e2f 3eT g =T gr 3e U U =U r 3eV =V r 3f 3eT =T 3e V U g =U gr r 3eV =V r 3TTf 3VT g

21、=T gr U =U r V =V reT g =T gr 1eU =U r 1eV =V r 1f 3VU =U r e 2V =V r e 2T g =T gr e2（43）f 3T gT g =T gr U =U r V =V rf =L 3T g eT g =T gr 1eU =U r 1eV =V r 1f 3T gU =U r e 2V =V r e 2T g =T gr e2f 3e T g =T gr 3T g e U =U r 3e V =V r 3T根据式（38）-（43），考虑图8-10中的所有三角元素和边界条件，得到：U i =0(i :炉壁号V i =0(i :炉壁号

22、P i =0(i =2U 26=U in 1, U 30=U in 2, （44） V 32=V in 1, V 34=V in 2, T g 26=0, T g 30=0, T g 32=0, T g 34=0假设该状态向量x ：（99×1）和控制向量u ：（4×1）为：x =x 1x 1=x 11x 2(991Tx 2=P 1x 21x 11=T g 1 T g 25x 31=U 3x 41=V 3u =U in 1我们得到：P 3 P 12(111T g 27T g 28T g 29T g 31T g 33x 31x 41(881T g 35(311x 21=T s

23、1 T s 25(251U 7V 7U in 2U 8V 8V in 1U 9 U 27V 9 V 27V in 2T (41U 28V 28（45）V 29(161U 29(161=Ax +Bu :E (9999, A (9999, B (998E xx =x 1E 110, E =00x 2A 12A B 1A =11, B =B A 21A 222（46）其中A 22=0，由于E 不是非奇异矩阵，式（46）不称为控制系统设计中广泛使用的状态方程。所以方程（46）称为广义系统14。因为在式（19）-（21）中P 没有时间分化，所以建立了广义系统。如果式（46）的广义系统脉冲可控且有限动态稳

24、定，就可能设计动态系统14。但是式（46）不是脉冲可控。所以，在下面的章节，我们会将式（46）的广义系统通过广义逆矩阵 转化为状态方程15，从而能设计控制系统。 图8 P 的摄动变量 图9 U 的摄动变量 图10 T s 的摄动变量 图11 一个典型的三角元素3.3 广义逆矩阵在广义系统中的应用15从式（46），我们得到A 22x 2=-A 21x 1-B 2u （47） 因为A 22=0，广义逆矩阵A 22为0，所以从式（46）和（47）我们知道：x 2=A 22(-A 21x 1-B 2u =0 （48）#然后代入式（46），我们得到状态方程：x 1=E 11A 11x 1+E 11B 1

25、u （49）-1-1A ，E 11B 1为B ，考虑了图3所示的仪表的位置，我们得到状态空间本文假设E 11A 11为模型：-1-1=Ax +Bu :A (8888, B (888xy =Cx :C (1388x =x 1(131T（45）u =U in 1U in 2V in 1T s 10V in 2T (41T s 11T s 15T s 21T s 25U 3V 3(131Ty =T s 1 T s 6因为式（50）是稳定和可观测的，可以使用线性控制理论。4. 反应炉控制系统本文采用线性二次高斯控制（LQG ）来处理下面描述的问题10。对于LQG 控制的详细解释见附录B 。反应炉模型的

26、LQG 控制系统见图12。LQG 控制系统通过式（50）来创建。其中K 是是反应炉模型估计的内部状态变量。下面是控制流优化调节增益，L 是卡尔曼滤波增益，x程：第一步：通过卡尔曼滤波和图3中仪表测量的数据来估计反应炉的内部状态变量。 第二步：基于优化调节和第一步估计的内部状态变量，修改风口吹气的速度来匹配反应炉模型所需的初始分配。 图12 LQG 优化控制系统其中风口吹气速度为（U in 1、U in 2、V in 1、V in 2），初始分布为（P (x , y , 0、V (x , y , 0、，需要的分布为（P r (x , y 、V r (x , y 、T g (x , y 、T s

27、r (x , y ）和T g (x , y , 0、T s (x , y , 0）r （P r 、V r 、T g 、T s r ）是式（2）-（16）的一个稳态解。r5. 数值实验在本研究中，反应炉模型的温度分布控制采用第2章中建立的反应炉模型。本章，将会提供LQG 控制的数值结果。下面是问题的描述：确定了风口吹气速度（U in 1、U in 2、V in 1、V in 2），因而在图13-16初始的分布（P (x , y , 0、V (x , y , 0、T g (x , y , 0、T s (x , y , 0）能被控制到图4-7中所需的分布。图15和16中初始的温度分布总的高于图6和图

28、7中所需的。 假设附录B 中的Q 、W 、R 、V 为：56elements 32elements (Q =W =diag 10, , 100. 1, , 0. 1:8888 R =V =7. 0106diag (1, 1, 1, 1因为本研究的目的是估计反应炉模型的内部温度，并和反应炉模型需要的内部温度匹配，因此定义了两个目标函数：(t = fi =1i 26, 30, 32, 3435(t +(t -T (t （51） T gi (t -T gi si sii =125(t 是在上述点通过(t 和T 其中T gi (t 和T si (t 是反应炉模型在图9和10中节点的温度，T si gi

29、LQG 控制得到的估计温度。方程（51）代表上述点的温度的估计误差。 f (t =rrTi =1j =11321n g i , j -T g i , j +s i , j -T s i , j （52）rnri =1j =11313其中T g i , j 和T s i , j 是在节点（i,j ）所需的温度，方程（52）代表所有节点温度误差。从图17和18可以看到风口吹风速度的改变。两个目标函数的变化见图19。从图17和18可以看出，风口吹气速度在接近15-的时候操作。当风口吹气，两个目标函数减少。优化的温度分布见图20和21。所以反应炉模型的内部温度变量就被估计出，且温度分布接近了通过LQG

30、 控制系统得到的所需的值。 图13 压力的初始分布 图15 气相温度的初始分布 图17 U in 1和V in 1 图14 气流的初始分布 图16 物料温度的初始分布 图18 U in 2和V in 2 图 19 f (t 和 f (t 的变化 图 20 气相温度的优化分布（t=50-） 气相温度的优化分布（ ） 图 21 材料温度的优化分布（t=50-） 材料温度的优化分布（ ） 6. 总结 本文建立了反应炉的仿真，能够用来计算炉内气相温度、固相温度、气体流量和压力 的分布。该仿真的测量值均来自炉壁附近和边界的控制输入。为了设计控制系统，状态空间 模型通过有限元分析来推导。通过该过程，线性控

31、制理论能够被用在控制系统的设计上。接 下来采用 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

32、 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 Differentia

33、l 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. Yaga

34、wa: 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 方法。这些计算会一直持续，直到满足收敛性条件。收敛 性条件的定义为： f p = max P (x, y, n p P(x, y, n p 1 < 0.0001 f V = max V (x, y, nV V ( x, y, nV 1 < 0.0001 f U = max U