面向本质安全化的化工过程设计：多稳态及其稳定性分析

第3章化工过程多稳态点的稳定性分析

3.1引言化工过程的稳定性在安全生产中发挥着重要作用，许多重大事故都是由于系统运行不稳定，人为操作失误，引发事故或者使得事故进一步恶化，进而造成财产损失和人员伤亡。有文献统计［1］表明，每10000起难于维持稳定的事故中，大约500起会造成财产的损失，有接近100起造成轻微的人身伤害，其中可能有1起事故会造成人员的严重伤害或死亡，如图31所示。由此可见，化工系统的稳定性在化工安全生产中占有重要位置，为减少重大事故的发生概率，一个有效的方法是提高系统的稳定性。

图31化工过程中各类事故的比例示意图

稳定性的概念出现在力学研究中，用于描述一个刚体运动的平衡状态。如果一个平衡状态是稳定的，那么这个刚体在受到干扰从原来位置微微移动后，最终仍能回到它原来的位置。反之，如果这个平衡状态不稳定，那么当这个刚体受到干扰的时候，它会趋于一个新位置，远离最初的平衡态，如图32所示。

图32力学中的稳定性示意图

在实际应用中，常常使用微分方程来描述系统的变化规律。在建立微分方程的过程中，我们只能考虑影响该过程的主要因素，忽略一些次要因素，而这些因素可以认为是干扰因素。干扰因素是不可避免的，可以瞬时起作用，也可以持续起作用。从数学上来看，瞬时干扰引起初值的变化，而持续干扰则会引起微分方程本身的变化。在某些系统中，系统初始条件或者微分方程的微小变化都会引起稳态解的巨大变化，因此，对于新设计的系统进行稳定性分析可以及时发现问题，避免设计方案实施后系统运行的不稳定。运动系统稳定性的概念是力学中平衡稳定性的扩展。李雅普诺夫定义下的运动稳定性理论主要研究微小干扰性因素对于系统运动的影响。微小的干扰因素普遍存在，不可避免，而且不确定。对于一些系统，微小的干扰因素的影响并不显著，因此，受干扰的运动与不受干扰的运动差别很小，这类运动系统称为稳定的； 对于另外一些运动，无论干扰多么小，随着时间的推移，受干扰的运动与不受干扰的运动总是相差巨大，这类运动系统称为不稳定的。由于干扰不可避免，所以运动稳定性的问题有着重要的理论和实际意义，在自然科学与工程技术领域受到了普遍关注。众所周知，化工过程是强非线性［2，3］过程，而非线性的系统通常存在多个稳态解［4~9］，这些稳态解的稳定性一般并不相同。在实际生产中，化工系统的操作条件受到人为操作、不确定因素的影响而不断变化，人们通常关心系统在某个操作条件不断变化时表现出来的特性，例如对于一个反应器，逐渐调整进料的流量，观察系统达到稳态时的特征。在现有的研究中，通过计算不同参数下系统的稳态解，然后判断每个稳态解的稳定特性，进而确定系统稳态解是否稳定。这种逐点判断稳定性的方法耗时巨大，本章将介绍通过奇异点划分区域快速判断稳定性区域的方法。本章首先介绍稳定性的概念和常用的判别方法，在此基础上提出通过奇异点划分区域快速判断稳定性区域的方法，之后使用1,3丙二醇厌氧发酵体系和苯乙烯聚合反应说明这种判断方法的有效性。3.2稳定性的概念化工过程是非线性很强的过程，在反应过程中伴随着物质的消耗和生成，同时有大量的能量释放和消耗。描述化工过程的系统中常常包含物料平衡、动量平衡以及能量平衡，这些过程都具有非线性的特性。为了更好地描述系统的特性，通常情况下使用动态方程来描述系统中变量随时间的变化关系，所谓动态方程即系统中的变量随时间变化的微分方程组。通常情况下，系统中具有可以改变的操作变量，随着操作条件的变化，系统的稳态解会发生变化，稳态解的稳定性也会发生变化。这样问题就抽象成为： 在含有参数的动态常微分方程组中，随着参数的不断变化，如何求解出系统中的稳态解并且判断其稳定性。首先给出稳定性严格的数学定义。对于微分方程F： DRn→Rn，如式(31)所示：

dxdt=F(x)x(0)=x0 (31

李雅普诺夫稳定性的原始定义［10］如下： 对于给定的初值x0，令x(x0,t)为微分方程的解。若对任何ε>0，存在δ>0，使得当初值x1满足|x1-x0|

|x(x1,t)-x(x0,t)|

则称x(x0,t)是李雅普诺夫稳定的。如果

limt→∞|x(x1,t)-x(x0,t)|=0(33

则称x(x0,t)是渐近李雅普诺夫稳定的。如果x(x0,t)是渐近李雅普诺夫稳定的，而且x0可以取任何值，则称x(x0,t)是全局渐近稳定的。李雅普诺夫稳定性主要研究系统稳态点在受到扰动时的运动特性，主要涉及稳定、渐近稳定、大范围渐近稳定和不稳定。而本章研究在参数条件变化时化工过程稳态点的稳定性，这些参数可以是操作参数，也可以是设计参数。研究前者可以分析系统在运行过程中对于扰动的耐受程度，后者可以在设计阶段就提高系统本身对于扰动的耐受程度。研究中将稳定性分为两类： 稳定，不稳定。带参数的化工过程系统的动态方程可以描述为

dxdt=F(x,λ) (34

系统的稳态解就是方程F(x,λ)=0的解，也就是说当λ=λ，x=x时，如果方程F(x,λ)=0，那么x就是系统(34)在λ=λ时的稳态解。稳定性与系统的稳态解有关，对于系统方程

dxdt=F(x,λ

x(0)=x Δx (35

已知x是系统的稳态解，Δx是小的扰动。此时，如果常微分方程的解x(t)满足limt→∞x(t)=x的条件，那么x(t)是稳定的解，相对应的x是一个稳定的稳态解，否则x就不是一个稳定的稳态解。应该指出，求解非线性方程组是一项比较困难的工作，通常不大可能求出解析解，因此常使用数值的方法进行求解。而数值方法在求解过程中存在两个问题： 一是迭代求解的过程有时不收敛，原因在于选取的初值不合适； 二是不能求出系统的所有解，这是由于非线性系统通常具有多个解，而在一般情况下迭代求解只能求出一个解。针对这个问题，上一章引入的同伦延拓法［6，11］可用来解决这个问题，该方法可以迅速求解动态非线性方程组的稳态解并且判定其稳定性。3.3稳定性的判断方法李雅普诺夫对稳定性问题提出了两种方法，及时种是级数展开法，第二种是通过构造李雅普诺夫函数来判断稳定性。虽然第二种方法由于不用求解稳态解得到了更广泛的应用，但是对于非线性系统而言，不存在构造李雅普诺夫函数的通用方法，而且对同一个问题可能构造出许多不同的李雅普诺夫函数，也可能很难构造出李雅普诺夫函数。因此，本书在判断单个稳态点的稳定性时采用李雅普诺夫及时方法。3.3.1李雅普诺夫判断方法由于李雅普诺夫及时方法在稳态点将系统线性展开，所以首先介绍线性微分系统的稳定性判断方法。对于线性微分方程组

dxidt=ai1x1 ai2x2 … ainxn(i=1,2,…,n) (36

矩阵形式如下：

X =AX(37

其中X=［x1x2…xn］T为n维向量，A为(aij)n×n的矩阵，那么方程组的解为

X=b eΛt (38

其中b=［b1b2…bn］T为n维向量，Λ是矩阵A的特征值。由于可能出现特征值相同的情况，不同特征值的重数记为n1,n2，…，nm，这里指代数重数。下面介绍几何重数，选取非奇异矩阵P，使得X=PY。对A做相似变换，问题转化为Y =P-1APY=JY。由于J与A相似，因此，具有相同的特征值和重数。不妨设J为约当(Jordan)标准型，即

J=J10…00J2…000…Jm (39

其中非对角线上子矩阵的元素为0，对角线上共有m个非零子矩阵，每个子矩阵也是ni×ni(i=1,2,…,m)对角型分块矩阵，即

Ji=Ji10…00Ji2…000…Jiα1(αi≤ni,i=1,2,…,m)(310

其中子矩阵

Jik=λi1000λi10001000λi(k=1,2,…,αi)(311

称为对应特征值的约当块，且满足ni1 ni2 … niαi=ni。特征值λi共有αi个约当块，其中1≤αi≤ni，αi称为几何重数。特别地，如果αi=ni，那么特征值λi的代数重数等于几何重数。对于线性系统X =AX： (1) 如果A的所有特征值具有负实部，即负实根或者负实部的复根，则系统是稳定的。(2) 如果特征值中有一个根有正实部，即正实根或者正实部的复根，则系统不稳定。(3) 如果没有带正实部的根，但是有实部为零的单根，即零根或一对纯虚根，则系统的解是稳定的，但不是渐近稳定。(4) 如果没有带正实部的根，但是有多重零根或多重虚根，此时每个重根的代数重数与几何重数相等，则系统为稳定的； 如果至少有一个重根的几何重数小于代数重数，那么系统是不稳定的。李雅普诺夫及时方法用于非线性的微分系统，通过坐标变换，将稳态点变换为零点，设微分方程为

dxidt=fi(x1,x2,…,xn)(i=1,2,…,n)(312

其中，fi(x1,x2,…,xn)为xi(i=1,2,…,n)的非线性函数，将fi在稳态点处线性展开，记为

fi(x1,x2,…,xn)=∑ni=1aijxj Xi(x1,x2,…,xn)(313

其中∑ni=1aijxj是一次近似项，Xi(x1,x2,…,xn)是高阶项。李雅普诺夫及时方法的定理如下： (1) 如果一次近似项的所有特征值都具有负实部，那么原非线性系统的稳态点是渐近稳定度，与高阶项无关。(2) 如果一次近似项的特征值至少有一个具有正实部，那么原非线性系统是不稳定的，与高阶项无关。(3) 如果一次近似项的有实部为零的特征值，而其余的特征值实部为负，那么原系统在稳态点的稳定性取决于高阶项，稳态点有可能稳定也有可能不稳定。3.3.2用奇异点判断系统的稳定性相比于传统的李雅普诺夫方法计算每一个点的稳定性，本章使用一种基于奇异点分析的方法来判断化工系统稳态点的稳定性，这种方法只判断部分点的稳定性就可以迅速判断出系统的稳定特性。如前所述，化工过程可以用动态方程(34)来表示，其中x是系统状态变量，x∈Rn，λ是可变化的参数。假设Fx为方程F(x,λ)的雅可比矩阵，x0为方程F(x,λ)=0在λ0处的系统的稳态解，即方程F(x,λ0)=0的解。那么当参数变化时，稳态解求解及稳定性判断的算法描述如下。(1) 计算系统的稳态解。求解不同λ值下的方程F(x,λ)=0的解x。(2) 计算系统雅可比矩阵的奇异值。对于每一个λ以及在该λ下的稳态解x，将它们代入Fx得到一个矩阵，判断矩阵Fx是否奇异，即矩阵是否满秩。如果矩阵奇异，那么记录下这个奇异点的数值λs和xs。在奇异点附近系统的稳定性可能发生变化，根据这个特性，通过判断奇异点两侧的稳态解的稳定性就可以确定被奇异点划分的区域的稳定性。(3) 将奇异值作为分界点，判断其两侧解的稳定性。计算λs左右两侧的稳定性，确定不同区域的稳定性。(4) 根据检验点的测试结果，划定稳定区域不同的稳定特性。重复以上过程计算出所有的奇异点并且判断出奇异点两侧的稳定性。详细的算法框图如图33所示。

图33用奇异点判断稳定性算法框图

对于复杂的问题，使用同伦延拓法求解系统存在的稳态解。对于复杂的系统尝试使用该算法来解决问题。3.4案例一发酵反应过程3.4.1发酵反应过程的数学模型

微生物发酵生产1，3丙二醇模型［12，13］如下，这是一个简化的模型，实际生产过程需要考虑更多的因素。这里仅仅使用该模型说明系统的稳定性及其判断方法。

dXdt=X(μ-D

dCsdt=D(Cs0-Cs)-Xqs

dCpdt=Xqp-DCp(314

其中，

μ=μmaxCsCs Ks1-CsCs1-CpCp

qs=ms μYms ΔqmsCsCs Ks

qp=mp Ympμ ΔqmpCsCs Kp(315

主要变量的含义： X为菌浓度，Cs为底物浓度，Cp为产物浓度，D为稀释速率。为方便计算起见，无因次化结果如下：

dxdτ=x(u-d

dydτ=d(y0-y)-x1β1 γ1u yy α1

dzdτ=x2β2 γ2u yy α2-zd

u=yy α3(1-y)(1-z)(316

其中主要符号的含义： x为无因次菌浓度，y为无因次底物浓度，z为无因次产物浓度，d为无因次稀释速度，y0为无因次进料底物浓度。3.4.2稳态点的稳定性判断首先计算稀释速度d=0.2时，无因次产物浓度z随无因次进料的底物浓度y0的变化关系，如图34所示，这里求出了该系统所有稳态解随参数变化的情况，在稳态解曲线上存在两个奇异点，它们是极限点(limit point)，用LP表示。

图34d=0.2时产物浓度随进料底物浓度的变化关系图

以下判断系统的稳定性。分析可知，该发酵反应体系中，稳态解的雅可比矩阵为3阶矩阵，有3个特征值。考虑到系统稳态点的稳定性与特征值的实部相关，我们绘制特征值实部随进料底物浓度的变化关系图，如图35、图36所示，主要关注特征值实部在底物进料浓度逐渐改变过程时的正负变化情况。

图35及时个特征值实部的变化曲线

图36第二个和第三个特征值实部的变化曲线

由图35和图36(a)可见，及时个和第二个特征值的实部在所有变化区域内均小于零，根据李雅谱诺夫及时方法判断可知，这两个特征值不对系统的稳定性产生决定性的影响。由图36(b)可见，第三个特征值的实部在经过两个LP点时符号发生变化，当第三个特征值的实部小于零时，系统的三个特征值的实部都小于零，此时，对应系统的稳态解是稳定的。随着操作参数变化，稳态点对应的雅可比矩阵的特征值也逐渐变化，当第三个特征值的实部大于零时，系统的三个特征值中有一个大于零，对应的稳态点为不稳定。由上述内容可知，系统的稳态点随操作参数的变化而变化，当参数经过奇异点后系统的稳态点的稳定性发生变化。奇异点将稳态解曲线划分为三段，这三段具有不同的稳定性。在这个例子中，三个特征值的奇异点是相关的，即特征值曲线都有两个奇异点，而且这两个奇异点对应的横坐标(操作变量)的数值