Matlab在化工计算及其可视化中的应用

游乐天

发贴心情
[转帖]Matlab在化工计算及其可视化中的应用

李俊杰
(襄樊学院国有资产处,湖北  襄樊  441053)
    [摘  要] Matlab是可以完成各种化工计算和数据处理的可视化软件，具有求解效率高、计算稳定、使用方便等特点。以常遇到的微分方程、曲线拟合、线性方程组等工程实际为例，分析了Matlab在化工计算及其运算结果可视化中的应用。
    [关键词]Matlab；化工计算；可视化
    [中图分类号] O 643               [文献标识码] A                    [文章编号] 1003-5095(2007）05-0045-03
    化工实验素有数据多、处理麻烦、计算复杂等问题，具体计算涉及到曲线拟合、解线性/非线性方程、解常微分/偏微分方程、插值和求积分等。用通常的数值计算方法解决这些问题比较复杂、工作量大，而采用Matlab语言进行编程来处理，就会很简单，只需调用有关的功能函数即可，同时Matlab还提供了方便灵活的数据输入方法，具有强大的数值运算能力，可方便地实现矩阵的操作运算，且图形可视化方法能使计算结果图形显示。本文通过实际算例，探讨Matlab在化工数值运算中的应用,以解决化工计算中数学求解的繁琐与困难。
1  常微分方程数值解
    已知:简单蒸馏苯－甲苯混合液。进料速率10  kmol/h ,   苯摩尔分率x0 =0.32,塔釜加热维持釜液量      20 kmol。物料系统服从Raoult定律，相对挥发度α=2.48。试求: 开始蒸馏时间t∈[0,2] h以内, 釜液中苯浓度x随时间t的变化关系。
    解算: 通过物料平衡和气液相平衡关系,得出常微分方程[1]:
x′(t)=(3.2-20.1x)/(20 + 29.6x)
    此方程可运用Matlab进行数值求解，其求解步骤分为两步：(1)创建常微分方程的函数M文件  fun.m；(2)采用四阶/五阶Ronge-Kutta方法，调用ode45( )函数求解，并绘图。具体程序和命令如下，结果如图1。
    % fun.m
    function dx =fun(t,x);
    dx = (3.2-20.1*x)/(20 + 29.6*x);
    在命令行输入如下内容：
    ts=0: 0.1: 2; %设定积分时间上下限
    x0 =0.32;   %设定初始值
    [t,x]=ode45(′fun′, ts, x0 );
    plot(t,x,′ *′); %绘制图形
    axis([0,2, 0.18,0.34]);
    grid,xlabel(′t′),ylabel(′x′)。


图1  常微分方程解
2  曲线拟合及其可视化
    数据图形能反映离散数据的许多内在本质或某些趋势，Matlab提供有多项式最小二乘拟合，并在数据可视化方面有卓越的表现。
    已知: 常压下乙醇－水的平衡数据如表1，液相摩尔分率x,气相摩尔分率y，求 y-x气液平衡曲线。
表1  实验原始数据[2]
x  0.0  2.0 4.0  6.0  8.0  10.0 14.0 18.0 25.0 30.0 35.0 40.0
y 0.0 1 7.5 27.3 34.0 39.2 43.0 48.2 51.3 55.1 57.5 59.5  61.4
x 45.0 50.0 55.0 60.0 65.0 70.0 75.0 80.0 85.0 90.0 95.0 100.0
y 63.5 65.7 67.8 69.8 72.5 75.5 78.5 82.0 85.5 89.8 94.2 100.0
    解算: 采用Matlab拟合曲线的准则是最小二乘法。调用多项式曲线拟合函数polyfit()对以上数据进行线性拟合、三次多项式拟合，并生成图形，将计算结果可视化。其程序为：
    x=[0.0  2.0  4.0  6.0  8.0  10.0 14.0 18.0 25.0 30.0 35.0 40.0 45.0 50.0 55.0 60.0 65.0 70.0 75.0 80.0 85.0 90.0 95.0 100];%构造向量矩阵x。
    y=[0.0  17.5 27.3 34.0 39.2 43.0 48.2 51.3 55.1 57.5 59.5 61.4 63.5 65.7 67.8 69.8 72.5 75.5 78.5 82.0 85.5 89.8 94.2 100]; %构造向量矩阵y。
    c1=polyfit(x,y,1)  %线性拟合系数矩阵
    a1=c1(1); a0=c1(2);  %取出线性拟合系数
    y1=a1*x+a0;  %产生一阶方程式
    s1=sum((y-y1).^2) %计算线性拟合方差
    c4=polyfit(x, y,4)  %三次多项式拟合系数
    a4=c4(1);a3=c4(2);a2=c4(3);a1=c4(4);
    a0=c4(5);  %取出三次多项式拟合系数
    y2=a4*x.^4+a3*x.^3+a2*x.^2+a1*x+a0;
    %构造三次方程式：
    s2=sum((y-y2).^2)  %三次多项式拟合方差
    plot(x,y,′ *′,x,y1,x,y2);  %绘图
    grid,xlabel(′x′ ),ylabel(′y′);
    gtext(′N=1′); gtext(′N=4′);
    运行程序得到结果如图2，线性拟合方程y1= 0.705 8x+28.716 5，方差为s1=1523.8；三次多项式拟合方程y2= 0.005 1x3-0.248 2x2+5.803 8 x+ 4.876 3，方差为s2=27.885 9，显然三次多项式拟合的结果要比线性拟合的好。


3  线性方程组的计算
    设某化合物的蒸气压与温度的关系为[3]：lgP=A+BlgT+C/T，试根据表中3组温度与蒸气压数据，求该化合物在T=15 ℃的蒸气压值。
表2  某化合物的蒸气压与温度数据
         序号   温度/℃ 蒸汽压/Pa  绝对温度/T
           1   -30       3.55×10-2    243.2
           2    0       1.90×10-1    273.2
           3   40         1.02    313.2
    解算: 由题意得下列线性方程组
    A +Blg(243.2) +C/243.2 =lg(3.55×10-2)
    A +Blg(273.2) +C/273.2 =lg(1.90×10-1)
    A +Blg(313.2) +C/313.2 =lg(1.02)
    此方程组的求解可直接在Matlab命令行中，通过矩阵表示法输入如下内容，即得参数A,B,C的值，于是T=15 ℃时该化合物的蒸气压可代入公式可求出。
    C=[1, log(243.2)/log(10), 1/243.2;
      1, log(273.2)/log(10), 1/273.2;
      1, log(313.2)/log(10), 1/313.2];
    D=[log(3.55e-2)/log(10);
      log(1.90e-1)/log(10);
      log(1.02)/log(10) ]。
    在Matlab命令行中输入：C \ D
    运行后得结果:
    ans =
    1.0 e+003
    0.0151  (A 的值)
    - 0.0035　 (B 的值)
    - 2.0017　 (C 的值)
    将A、B、C和T=15 ℃的值代入公式lgP=A+BlgT+C/T，并做等效变形，可求得:
    P=exp((0.0151-0.0035×log(288.2)-2.0017/288.2)×log(10))
    即P=0.9990。
4  结束语
    Matlab作为新一代的科学和工程计算语言,有着传统的计算机语言所无法比拟的优点，不需了解求解的算法，对于化工计算过程，具有开放性、可视性、易操作、易推广的特点。本文只是对Matlab在微分方程求解、曲线拟合、解线性方程组中的应用做了一点介绍，随着使用者水平的提高和对软件应用领域的开发, 它将发挥更大的作用，成为化工计算中教师和学生的得力助手。