《化学二元反应蒸馏塔仿真mathematica实现程序.docx》由会员分享，可在线阅读，更多相关《化学二元反应蒸馏塔仿真mathematica实现程序.docx（11页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、(*已知参数*)F = 100;xF = 0.5;xD = 0.95;xW = 0.05;q = 1;han = 6944;af = 2; DD = 50;W = 50;erro = 10-8;time = 20;cur = 15;fun = NSolveyq*(q - 1) = q*xq - xF, yq*(af - 1)*xq + 1) = af*xq, xq, yq; Rmin = FromDigits(xD - Evaluateyq /. fun)/(Evaluateyq /. fun - Evaluatexq /. fun);(*控制部分*)(*可知F=D +W,所以w=0,newt
2、on法在此建立，即w-50=0,而w-time=LN-VN=L1+qF-Q/han-time=0=R*D+q*F-Q/han-50=0,即： R*D+q*F-Q/han-W=0 ,其中，R,Q的值影响整个反应体系，具体分离目标和它俩的取值有关，利用newton-rapphson迭代法xn+1=xn-f (xn)/f(xn)求取合适的R和Q 。因为R的取值范围一般为R=1.12.0Rmin,而Rmin=(xD-Yq)/(Yq-Xq) ,因此尽量满足R的值在此范围内*)pc = 1/3;ic = 1;dc = 1;R = 2;erropre1 = 5;erropre2 = 5;errold = 0
3、;flag = 1;caltime = 0;failtime = 0;Whileflag != 0, +caltime; L1 = R*DD; LN = L1 + F*q; Q = (LN - W)*han; VN = Q/han; V1 = VN + F*(1 - q); (*物料平衡微分方程36块板，一个冷凝器，一个再沸器，共38个微分方程*) function = NDSolve (*初始化方程，初始值可以任意给定*) x10 = 0, x20 = 0, x30 = 0, x40 = 0, x50 = 0, x60 = 0, x70 = 0, x80 = 0, x90 = 0, x100
4、 = 0, x110 = 0, x120 = 0, x130 = 0, x140 = 0, x150 = 0, x160 = 0, x170 = 0, x180 = 0, x190 = 0, x200 = 0, x210 = 0, x220 = 0, x230 = 0, x240 = 0, x250 = 0, x260 = 0, x270 = 0, x280 = 0, x290 = 0, x300 = 0, x310 = 0, x320 = 0, x330 = 0, x340 = 0, x350 = 0, x360 = 0, x370 = 0, x380 = 0, (*冷凝器方程*) x1t
5、= V1*(af*x2t)/(af - 1)*x2t + 1) - L1*x1t - DD*x1t, (*1-36块板的物料平衡方程*) (*精馏段方程*) x2t = V1*(af*x3t)/(af - 1)*x3t + 1) + L1*x1t - L1*x2t - V1*(af*x2t)/(af - 1)*x2t + 1), x3t = V1*(af*x4t)/(af - 1)*x4t + 1) + L1*x2t - L1*x3t - V1*(af*x3t)/(af - 1)*x3t + 1), x4t = V1*(af*x5t)/(af - 1)*x5t + 1) + L1*x3t -
6、L1*x4t - V1*(af*x4t)/(af - 1)*x4t + 1), x5t = V1*(af*x6t)/(af - 1)*x6t + 1) + L1*x4t - L1*x5t - V1*(af*x5t)/(af - 1)*x5t + 1), x6t = V1*(af*x7t)/(af - 1)*x7t + 1) + L1*x5t - L1*x6t - V1*(af*x6t)/(af - 1)*x6t + 1), x7t = V1*(af*x8t)/(af - 1)*x8t + 1) + L1*x6t - L1*x7t - V1*(af*x7t)/(af - 1)*x7t + 1),
7、 x8t = V1*(af*x9t)/(af - 1)*x9t + 1) + L1*x7t - L1*x8t - V1*(af*x8t)/(af - 1)*x8t + 1), x9t = V1*(af*x10t)/(af - 1)*x10t + 1) + L1*x8t - L1*x9t - V1*(af*x9t)/(af - 1)*x9t + 1), x10t = V1*(af*x11t)/(af - 1)*x11t + 1) + L1*x9t - L1*x10t - V1*(af*x10t)/(af - 1)*x10t + 1), x11t = V1*(af*x12t)/(af - 1)*x
8、12t + 1) + L1*x10t - L1*x11t - V1*(af*x11t)/(af - 1)*x11t + 1), x12t = V1*(af*x13t)/(af - 1)*x13t + 1) + L1*x11t - L1*x12t - V1*(af*x12t)/(af - 1)*x12t + 1), x13t = V1*(af*x14t)/(af - 1)*x14t + 1) + L1*x12t - L1*x13t - V1*(af*x13t)/(af - 1)*x13t + 1), x14t = V1*(af*x15t)/(af - 1)*x15t + 1) + L1*x13t
9、 - L1*x14t - V1*(af*x14t)/(af - 1)*x14t + 1), x15t = V1*(af*x16t)/(af - 1)*x16t + 1) + L1*x14t - L1*x15t - V1*(af*x15t)/(af - 1)*x15t + 1), x16t = V1*(af*x17t)/(af - 1)*x17t + 1) + L1*x15t - L1*x16t - V1*(af*x16t)/(af - 1)*x16t + 1), x17t = V1*(af*x18t)/(af - 1)*x18t + 1) + L1*x16t - L1*x17t - V1*(a
10、f*x17t)/(af - 1)*x17t + 1), x18t = V1*(af*x19t)/(af - 1)*x19t + 1) + L1*x17t - L1*x18t - V1*(af*x18t)/(af - 1)*x18t + 1), x19t = V1*(af*x20t)/(af - 1)*x20t + 1) + L1*x18t - L1*x19t - V1*(af*x19t)/(af - 1)*x19t + 1), x20t = V1*(af*x21t)/(af - 1)*x21t + 1) + L1*x19t - L1*x20t - V1*(af*x20t)/(af - 1)*x
11、20t + 1), x21t = V1*(af*x22t)/(af - 1)*x22t + 1) + L1*x20t - L1*x21t - V1*(af*x21t)/(af - 1)*x21t + 1), (*第21块板为进料板，即第22个方程*) x22t = VN*(af*x23t)/(af - 1)*x23t + 1) + L1*x21t - LN*x22t - VN*(af*x22t)/(af - 1)*x22t + 1) + F*xF, (*提馏段方程*) x23t = VN*(af*x24t)/(af - 1)*x24t + 1) + LN*x22t - LN*x23t - VN
12、*(af*x23t)/(af - 1)*x23t + 1), x24t = VN*(af*x25t)/(af - 1)*x25t + 1) + LN*x23t - LN*x24t - VN*(af*x24t)/(af - 1)*x24t + 1), x25t = VN*(af*x26t)/(af - 1)*x26t + 1) + LN*x24t - LN*x25t - VN*(af*x25t)/(af - 1)*x25t + 1), x26t = VN*(af*x27t)/(af - 1)*x27t + 1) + LN*x25t - LN*x26t - VN*(af*x26t)/(af - 1
13、)*x26t + 1), x27t = VN*(af*x28t)/(af - 1)*x28t + 1) + LN*x26t - LN*x27t - VN*(af*x27t)/(af - 1)*x27t + 1), x28t = VN*(af*x29t)/(af - 1)*x29t + 1) + LN*x27t - LN*x28t - VN*(af*x28t)/(af - 1)*x28t + 1), x29t = VN*(af*x30t)/(af - 1)*x30t + 1) + LN*x28t - LN*x29t - VN*(af*x29t)/(af - 1)*x29t + 1), x30t
14、= VN*(af*x31t)/(af - 1)*x31t + 1) + LN*x29t - LN*x30t - VN*(af*x30t)/(af - 1)*x30t + 1), x31t = VN*(af*x32t)/(af - 1)*x32t + 1) + LN*x30t - LN*x31t - VN*(af*x31t)/(af - 1)*x31t + 1), x32t = VN*(af*x33t)/(af - 1)*x33t + 1) + LN*x31t - LN*x32t - VN*(af*x32t)/(af - 1)*x32t + 1), x33t = VN*(af*x34t)/(af
15、 - 1)*x34t + 1) + LN*x32t - LN*x33t - VN*(af*x33t)/(af - 1)*x33t + 1), x34t = VN*(af*x35t)/(af - 1)*x35t + 1) + LN*x33t - LN*x34t - VN*(af*x34t)/(af - 1)*x34t + 1), x35t = VN*(af*x36t)/(af - 1)*x36t + 1) + LN*x34t - LN*x35t - VN*(af*x35t)/(af - 1)*x35t + 1), x36t = VN*(af*x37t)/(af - 1)*x37t + 1) +
16、LN*x35t - LN*x36t - VN*(af*x36t)/(af - 1)*x36t + 1), x37t = VN*(af*x38t)/(af - 1)*x38t + 1) + LN*x36t - LN*x37t - VN*(af*x37t)/(af - 1)*x37t + 1), (*再沸器物料平衡方程*) x38t = LN*x37t - W*x38t - VN*(af*x38t)/(af - 1)*x38t + 1) , (*待求未知量列表*) x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x
17、16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30, x31, x32, x33, x34, x35, x36, x37, x38, (*微分时间变量设置*) t, 0, 50; (*控制部分*) (* IfR=1.1*Rmin&R=2.0*Rmin,R=R-(R*DD+q*F-Q/han-W)/DD,Q=Q-(R*DD+q*F-Q/han-W)*(-1/han),R=(Q/han+W-F*q)/DD,*) (*算法在第二个if语句中，*) erro1 = FromDigitsxD - Evaluat
18、ex1time /. function; erro2 = FromDigitsxW - Evaluatex38time /. function; IfAbserro1 = erro & Abserro2 Abserropre1 & Abserro1 Abserropre1, +failtime, Iffailtime 0, Print不能实现控制 ., Break, failtime = 0, erropre1 = erro1, erropre2 = erro2 计算的次数 - caltime最小回流比 - Rmin计算的回流比 - R计算的回热量 - Q塔顶浓度误差 - erro1塔底浓度误
19、差 - erro2精度是否可达，存在扰动？ - failtime (* 清除已知变量的值 Clear* *) (*结果处理与显示*)fun2 = NSolve T1 = 3862/(11.6531 - Log9/(x1time /. function)*(E0.6932 - 1) + 1), T2 = 3862/(11.6531 - Log9/(x2time /. function)*(E0.6932 - 1) + 1), T3 = 3862/(11.6531 - Log9/(x3time /. function)*(E0.6932 - 1) + 1), T4 = 3862/(11.6531
20、- Log9/(x4time /. function)*(E0.6932 - 1) + 1), T5 = 3862/(11.6531 - Log9/(x5time /. function)*(E0.6932 - 1) + 1), T6 = 3862/(11.6531 - Log9/(x6time /. function)*(E0.6932 - 1) + 1), T7 = 3862/(11.6531 - Log9/(x7time /. function)*(E0.6932 - 1) + 1), T8 = 3862/(11.6531 - Log9/(x8time /. function)*(E0.
21、6932 - 1) + 1), T9 = 3862/(11.6531 - Log9/(x9time /. function)*(E0.6932 - 1) + 1), T10 = 3862/(11.6531 - Log9/(x10time /. function)*(E0.6932 - 1) + 1), T11 = 3862/(11.6531 - Log9/(x11time /. function)*(E0.6932 - 1) + 1), T12 = 3862/(11.6531 - Log9/(x12time /. function)*(E0.6932 - 1) + 1), T13 = 3862
22、/(11.6531 - Log9/(x13time /. function)*(E0.6932 - 1) + 1), T14 = 3862/(11.6531 - Log9/(x14time /. function)*(E0.6932 - 1) + 1), T15 = 3862/(11.6531 - Log9/(x15time /. function)*(E0.6932 - 1) + 1), T16 = 3862/(11.6531 - Log9/(x16time /. function)*(E0.6932 - 1) + 1), T17 = 3862/(11.6531 - Log9/(x17tim
23、e /. function)*(E0.6932 - 1) + 1), T18 = 3862/(11.6531 - Log9/(x18time /. function)*(E0.6932 - 1) + 1), T19 = 3862/(11.6531 - Log9/(x19time /. function)*(E0.6932 - 1) + 1), T20 = 3862/(11.6531 - Log9/(x20time /. function)*(E0.6932 - 1) + 1), T21 = 3862/(11.6531 - Log9/(x21time /. function)*(E0.6932
24、- 1) + 1), T22 = 3862/(11.6531 - Log9/(x22time /. function)*(E0.6932 - 1) + 1), T23 = 3862/(11.6531 - Log9/(x23time /. function)*(E0.6932 - 1) + 1), T24 = 3862/(11.6531 - Log9/(x24time /. function)*(E0.6932 - 1) + 1), T25 = 3862/(11.6531 - Log9/(x25time /. function)*(E0.6932 - 1) + 1), T26 = 3862/(1
25、1.6531 - Log9/(x26time /. function)*(E0.6932 - 1) + 1), T27 = 3862/(11.6531 - Log9/(x27time /. function)*(E0.6932 - 1) + 1), T28 = 3862/(11.6531 - Log9/(x28time /. function)*(E0.6932 - 1) + 1), T29 = 3862/(11.6531 - Log9/(x29time /. function)*(E0.6932 - 1) + 1), T30 = 3862/(11.6531 - Log9/(x30time /
26、. function)*(E0.6932 - 1) + 1), T31 = 3862/(11.6531 - Log9/(x31time /. function)*(E0.6932 - 1) + 1), T32 = 3862/(11.6531 - Log9/(x32time /. function)*(E0.6932 - 1) + 1), T33 = 3862/(11.6531 - Log9/(x33time /. function)*(E0.6932 - 1) + 1), T34 = 3862/(11.6531 - Log9/(x34time /. function)*(E0.6932 - 1
27、) + 1), T35 = 3862/(11.6531 - Log9/(x35time /. function)*(E0.6932 - 1) + 1), T36 = 3862/(11.6531 - Log9/(x36time /. function)*(E0.6932 - 1) + 1), T37 = 3862/(11.6531 - Log9/(x37time /. function)*(E0.6932 - 1) + 1), T38 = 3862/(11.6531 - Log9/(x38time /. function)*(E0.6932 - 1) + 1), T1, T2, T3, T4,
28、T5, T6, T7, T8, T9, T10, T11, T12, T13, T14, T15, T16, T17, T18, T19, T20, T21, T22, T23, T24, T25, T26, T27, T28, T29, T30, T31, T32, T33, T34, T35, T36, T37, T38;rat = x1time, x2time, x3time, x4time, x5time, x6time, x7time, x8time, x9time, x10time, x11time, x12time, x13time, x14time, x15time, x16t
29、ime, x17time, x18time, x19time, x20time, x21time, x22time, x23time, x24time, x25time, x26time, x27time, x28time, x29time, x30time, x31time, x32time, x33time, x34time, x35time, x36time, x37time, x38time /. function;tem = T1, T2, T3, T4, T5, T6, T7, T8, T9, T10, T11, T12, T13, T14, T15, T16, T17, T18,
30、 T19, T20, T21, T22, T23, T24, T25, T26, T27, T28, T29, T30, T31, T32, T33, T34, T35, T36, T37, T38 /. fun2;flow1 = L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, L1, LN, LN, LN, LN, LN, LN, LN, LN, LN, LN, LN, LN, LN, LN, LN, LN, LN, W;flow2 = V1, V1, V1, V1, V1, V1, V1
31、, V1, V1, V1, V1, V1, V1, V1, V1, V1, V1, V1, V1, V1, VN, VN, VN, VN, VN, VN, VN, VN, VN, VN, VN, VN, VN, VN, VN, VN, VN, VN;Plotx1t /. function, x38t /. function, t, 0, time, PlotStyle - Red, Blue, PlotRange - All, AxesLabel - 时间, 浓度x1time, x38time /. functionShowListPlottem, Joined - True, ListPlo
32、ttem, Filling - Axis, PlotStyle - Red, PointSizeLarge, AxesLabel - 塔板编号, 温度ShowListPlotaf*rat/(af - 1)*rat + 1), Joined - True, ListPlotaf*rat/(af - 1)*rat + 1), PlotStyle - Red, PointSizeLarge, ListPlot1 - af*rat/(af - 1)*rat + 1), Joined - True, ListPlot1 - af*rat/(af - 1)*rat + 1), PlotStyle - Bl
33、ack, PointSizeLarge, AxesLabel - 塔板编号, A、B汽相浓度ShowListPlotrat, Joined - True, ListPlotrat, PlotStyle - Red, PointSizeLarge, ListPlot1 - rat, Joined - True, ListPlot1 - rat, PlotStyle - Black, PointSizeLarge, AxesLabel - 塔板编号, A、B液相浓度ShowListPlotflow1, Joined - True, ListPlotflow1, PlotStyle - PointSizeLarge, Red, ListPlotflow2, Joined - True, ListPlotflow2, PlotStyle - PointSizeLarge, Black, AxesLabel - 塔板编号，气液相流量计算的次数 - 43最小回流比 - 1.7计算的回流比 - 1.3514计算的回热量 - .塔顶浓度误差 - -8.07045*10-9塔底浓度误差 - 7.94779*10-9精度是否可达，存在扰动？ - 00.95,0.05