第六章第六章 化工过程动态模拟与分析化工过程动态模拟与分析主要内容主要内容一、化工过程的动态模型一、化工过程的动态模型二、二、CSTR的动态模型的动态模型三、状态空间分析法三、状态空间分析法四、常微分方程（组）求解四、常微分方程（组）求解6.1 化工过程的动态模型化工过程的动态模型1、动态特性是化工过程系统最基本的特性之一、动态特性是化工过程系统最基本的特性之一间歇过程、开（停）工过程、系统的老化 2、化工过程动态模型分类、化工过程动态模型分类模型类型集中参数模型分布参数模型数学形式代数常微分方程组代数偏微分方程组应用实例理想搅拌罐反应器动态模型管式反应器动态模型多级集中参数模型代数常微分方程组混合模型串联CSTR 动态模型上述方程组的混合形式 多个单元过程组合而成的系统3、确定性动态模型的数学处理、确定性动态模型的数学处理（1）正问题模型方程组的求解；）正问题模型方程组的求解；（2）逆问题模型参数的估计；）逆问题模型参数的估计；（3）过程系统的定性分析。）过程系统的定性分析。定常态的多重性和稳定性、参数敏感性定常态的多重性和稳定性、参数敏感性 6.2 常微分方程的数值求解常微分方程的数值求解一、一阶常微分方程初值问题一、一阶常微分方程初值问题?dy? f (x,y)?dx?y(a) ? y0x? a,b计算出解函数 y(x) 在一系列节点 a = x0 x1 xn= b 处的近似值yi? y(xi)(i ? 1, . ,n)节点间距hi?xi?1?xi(i ?0,.,n?1)为步长，通常采用等距节点，即取 hi= h (常数)。I 欧拉法欧拉法?欧拉公式： 向前差商近似导数y(x1)? y(x0)y? (x0) ?h记为y(x1) ? y(x0)?hy? (x0) ? y0?h f (x0, y0)x0x1y1y1y2yi?1? yi?h f (xi, yi)(i ?0,.,n?1)y0x0x1x2II 隐式欧拉法隐式欧拉法向后差商近似导数y(x1)? y(x0)y? (x1) ?hy(x1)? y0?h f (x1, y(x1)迭代求解迭代求解x0x1yi?1? yi?h f (xi?1, yi?1)(i ? 0,.,n?1)一般先用显式计算一个初值，再迭代求解。III 改进欧拉法改进欧拉法yi?1? yi? h f (xi, yi)hyi?1? yi? f (xi, yi)? f (xi?1, yi?1)(i ? 0,.,n?1)2IV 龙格库塔法龙格库塔法将每两点之间分为若干点，以各点的处的斜率将每两点之间分为若干点，以各点的处的斜率平均值作为中值，并根据误差分析，可以给出下列平均值作为中值，并根据误差分析，可以给出下列常微分初值问题的数值计算公式：常微分初值问题的数值计算公式：yi?1?K1K2K3K4?yi? (K1?2K2?2K3?K4)h6f(xi, yi)hhf(xi?2, yi?2K1)f(xi? , yi? K2)f(xi?h, yi?hK3)h2h2例 6.1求解如下常微分初值问题：求解如下常微分初值问题：dT? 0.146?0.095?cos?t 12?0.0114TdtT?t ? 0? 93.3改进欧拉法程序：改进欧拉法程序：t = 0:10:100; h = 10; T0 = 93.3;Ts = zeros(1,11); Ts(1) = T0;for i = 1:10T_ = Ts(i) + h*f(t(i), Ts(i) ;Ts(i+1) = Ts(i) + 0.5*h*( f( t(i), Ts(i) ) + f( t(i+1), T_ ); end龙格龙格-库塔法程序：库塔法程序：t = 0:10:100; h = 10; T0 = 93.3;Ts = zeros(1,11); Ts(1) = T0;for i = 1:10xi = t(i); yi = Ts(i);K1 = f(xi,yi);K2 = f(xi+h/2, yi+h*K1/2);K3 = f(xi+h/2, yi+h*K2/2);K4 = f(xi+h, yi+h*K3);Ts(i+1) = yi + h*(K1+2*K2+2*K3+K4)/6;end二、一阶常微分方程组初值问题二、一阶常微分方程组初值问题一阶常微分方程初值问题的改进欧拉法：一阶常微分方程初值问题的改进欧拉法：y? y ? h f (x , y )i?1iiihyi?1? yi? f (xi, yi)? f (xi?1, yi?1)(i ? 0,.,n?1)2?y ? f?x,y,z?, y?x0? y0一阶常微分方程组：一阶常微分方程组：?z ? g?x,y,z?, z?x0? z0?yi?1? yi?hf?xi,yi,zi?zi?1? zi?hg?xi,yi,zi?h?y? y ?f x ,y ,z ? f x,y,iiiii?1i?1zi?1?i?12?h?zi?1? zi?g?xi,yi,zi? g?xi?1,yi?1,zi?1?2?一阶常微分方程组初值问题的改进欧拉法一阶常微分方程组初值问题的改进欧拉法例 6.2某一温度为某一温度为T=360 K的间歇式反应器中进行反应的间歇式反应器中进行反应A B C，设，设xA、xB分别为分别为A、B的浓度，反应速率方的浓度，反应速率方程如下，试求时间段程如下，试求时间段 0，1反应器中反应器中A、B浓度变化情浓度变化情况。况。?dxA2? ?4000?exp ?2500T ?x ,x 0 ?1AA?dt?dx2?B? 4000? exp?2500T?xA?620000?exp?5000T?xB,xB?0? 0? dtxA0 = 1; xB0 = 0; % 初始值t = 0:0.1:1; % 时间划分h = 0.1; % 步长A = zeros(1,11); A(1) = xA0;B = zeros(1,11); B(1) = xB0;for i = 1:10dA,dB = fs(A(i),B(i);A_ = A(i) + h*dA;B_ = B(i) + h*dB;dA_,dB_ = fs(A_, B_);A(i+1) = A(i) + 0.5*h*(dA + dA_);B(i+1) = B(i) + 0.5*h*(dB + dB_);endfunction dA,dB = fs(xA,xB)T = 300; % 反应温度dA = -4e3*exp( -2500/T )*xA2;dB = 4e3*exp( -2500/T )*xA2 -620000*exp( -5000/T )*xB;三、高阶常微分方程组初值问题三、高阶常微分方程组初值问题思路 转化为一阶常微分方程组求解转化为一阶常微分方程组求解?y ? f x,y,y?y?x0? y0,y?x0? y0?z? y?y ? z, y?x0? y0?z? f?x,y,z?, z?x0? y0可用改进欧拉法或龙格库塔法求解。可用改进欧拉法或龙格库塔法求解。6.3 CSTR 的动态特性的动态特性例249.9608 249.9666 249.9715 249.9757 249.97936.4求解常微分初值问题：如图所示储水罐，底面积求解常微分初值问题：如图所示储水罐，底面积为为A2.5，初始高度，初始高度H0=1，进水量，进水量Fi=10，k=0.04，求，求液面高度液面高度H随时间变化关系。随时间变化关系。FidHFik100.04根据质量守恒：根据质量守恒：dt?A?AH ?2.5?2.5?HHA参数参数Fi，k，AFo=k不同，最终解不同不同，最终解不同HdHFik?dt?A?A?H ? 0? H ? 250yi?1? yi?h f (xi?1, yi?1)(i ? 0,.,n?1)CSTR模型模型cA, f, FcA, F根据物料、能量守恒：根据物料、能量守恒：R? k0e?E RT?CA, fdcA?E RT?V? F CA,f?CA?V ?k0e? CA,fdtVdTQV ?CP?dt? F ?CP?Tf?T?UA?T ?TC?VR?H?ABTc为冷却剂温度，R为反应速度U、A为反应物与冷却剂之间总传热系数和总换热面积理想混合液相均相反映给定系统参数系统参数：cA, f, V, F, U, S, Tc, E, Cp 和r给定初始条件：t=0时，cA=cA0, T = T0, 则可求解：TTs显然，初始条件和设备、物性参显然，初始条件和设备、物性参数决定求解轨迹。数决定求解轨迹。0t若在数值求解时，不限制最终时间，则何时程序终止？若在数值求解时，不限制最终时间，则何时程序终止？当程序自动终止时的状态即为定常态当程序自动终止时的状态即为定常态 (Steady State)当系统参数已知当系统参数已知，cA和T为未知量，则定常态的求解：dcA?E RT?V? F?CA, f?CA?V ?k0e?CA, f? 0dtdTV ? CP? F ? CP?Tf?T?UA?T ?TC?VR?H? 0dt显然，只要系统参数确定，则上述非线性方程组显然，只要系统参数确定，则上述非线性方程组的解（的解（系统的定常态系统的定常态）只与选定的初始解（）只与选定的初始解（ 初始条件初始条件）有关。有关。多重定常态：非线性方程组有多个解多重定常态：非线性方程组有多个解定常态的稳定性分析：这些定常态皆稳定吗？定常态的稳定性分析：这些定常态皆稳定吗？该图何用？该图何用？该该相相图图如如何何做做？(cA0, T0)相轨迹相轨迹多多维维问问题题应应如如何何分分析析？上述上述CSTR中可能出现的相图（状态空间分析）中可能出现的相图（状态空间分析）定常态定常态A, B, C对应同一系统参数、不同初始条件对应同一系统参数、不同初始条件