MATLAB在热物理学中的应用

MATLAB在热物理学中的应用

学号：2010110114 姓名：贾海龙专业：物理学1011 摘要：本文阐述了基于MATLAB的数值计算、可视化图形处理、开放式以及可扩充体系结构的特点,并介绍了高性能语言 MATLAB 在大学物理热物理学中的一些应用，包括在麦克斯韦速率分布、理想气体定容比热回归分析和化工热力学中的应用。

关键词：MATLAB;麦克斯韦速率分布;理想气体;热力学

Application of MATLAB in thermal physics Abstract:Based on MATLAB’s features of numerical calculation, visualization of graphics processing,opening and scalable architecture,introduced the applications of language of MATLAB with high-performance in thermal physics of university physics,include in the Maxwell speed distribution、specific heat at constant volume regression analysis of ideal gas,and the chemical industry thermodynamics.

Key W ords: MATLAB; Maxwell speed distribution;ideal gas; thermodynamic

引言

大学物理学是工科学生的一门必修课, 由于大学物理数学处理比较复杂,恰当地使用可视化以展现数学公式的物理图像, 使其变得直观、形象。MATLAB 是一套高性能的数值计算和可视化软件,下面我们从一些典型的实例出发, 介绍 MATLAB 在热物理学方面的具体应用。

1 MA TLAB-PDEtool介绍

MATLAB-PDEtool提供了一个功能强大使用灵活的二维有限元偏微分方程求解环境，其图形用户界面更是使用十分方便、直观一般来说，MATLAB-PDEtoo l包括3个步骤：

（1）定义一个PDE问题，它包括确定二维求解区域、边界条件和PDE系数。MATLAB-PDEtool能够求解的PDE型式有：椭圆型、抛物线型、双曲线型、特征值型。当使用GUI时，可以在画图模式下确定求解区域；在边界模式下选择方程形式和设置方程系数。

（2）数值求解，它包括剖分、离散方程和得到一个数值解。在GUI中，在剖

分模式下形成满意的网格；在求解模式下通过选择数值计算方法求解。

（3） 图形化显示结果。通常用于表现有限元计算结果的图形有：变形网格图、云图、等值线图、矢量图、网格图、表面图、流线图等。

2 MA TLAB 在麦克斯韦速率分布中的应用

气体动力学理论中麦克斯韦速率分布律是大学物理讲授与学习中的一个难点和重点。这是因为公式比较复杂抽象，数学推导证明比较繁琐。如果借助 Matlab 就可以比较方便地解决这些问题[1]。首先, 推导三种速率和归一化条件。已知分布函数表达式为：

KT

mV

v e

v kT

m f

222

3

)

()

2(

4-=ππ （1）

最大概然速率分布可由下式求出： 0)(=dv

df v （2）

平均速率的定义是：

dv vf v o

v ?

∞

=)(

（3）

方均根速率为： 2

1

)(0

2

)

(2

dv f v

v rms

v

v

?

∞

==

（4）

归一化条件是：

dv vf v ?

∞

)( （5）

后三项求解比较复杂,其中用到Gamma 函数,传统方法是查数学用表得到结果。如果应用 Matlab 的符号计算功能, 只需要简单几行语句就可以解决这些问题。用到相关函数有：符号变量创建函数 syms ，求微分函数diff ， 求积分函数 int ，符号化简函数 simple ，字符串转化函数 eval 。

具体 Matlab 语句如下：

Syms mktvanp

g=exp(- (m*v^2)/(2*k*t))*v^2；

a=m/(2*k*t);f=int(exp(- a*v^2)*v^n,v ,0,inf)； b=f*4*pi*(a/pi)^(3/2)；

pretty(solve(diff(g,v)))； %求解最概然速率 n=2;eval(simple(eval(b)))，%证明归一化条件 n=3;pretty(simple(eval(b)))，%求解平均速率

n=4;pretty(eval(simple(eval(b^0.5))))， %求解方均根速率 程序中第五行语句给出最大概然速率： m

kT v p

2=

（6）

其物理意义是“若把整个速率范围分成许多相等的小区间，则 vp 所在的区间的分子数占分子总数的百分比最大”[2]。所以可知在( 0，+∞) 速率区间的分子数占分子总数的百分比(对应着曲线与X 轴所围的面积) 恒等于1，即满足归一化条件，第六行语句证明了这个结论。

第七行语句给出平均速率：

m

kT v 8=

（7）

第八行语句给出方均根速率：

m

kT v

rms

3=

（8）

其次,运用 Matlab 强大的画图功能，可以画出同一种气体分子( 氮气) 在不同温度下的分布曲线( 图 1)和不同气体分子( 二氧化碳、氧气、甲烷) 在同一温度下的曲线( 图 2)。程序中用到画图函数 ezplot 。

图1、图2中实线为麦克斯韦速率分布曲线，虚线的 X 轴坐标为最大概然速率

p v

的值。从图中曲线变化情况可以验证p v 是气体分子质量 m 的减函数， 是温度 T

的增函数。由此可以说明在满足归一化条件下，温度升高时曲线变得平坦些，并向高速区域扩展。即温度越高，速率大的分子越多，这就是通常所说的温度越高，分子运动越剧烈的真正含义。

3基于MATLAB 的理想气体定容比热容的回归分析

比热容是气体重要的热力学性质之一，在工程热力学和化工热力学的计算中，常常需要用到各种理想气体的定容热容，并通常使用某种数值计算公式进行计算[3]。为了在工程应用中能更准确地进行热力学计算，本文用MATLAB 对理想气体的定容比热容计算公式进行回归分析。 3.1 已有计算公式存在的问题

分析各种相关的文献，理想气体的定容比热容的计算公式存在公式不统一，且计算结果误差偏大的问题。以氧气为例，在273K~1 800K 的温度范围内，文献[4]中给出的计算公式为：

)/(10

312.110

062.510

2.15167.173

9

2

6

3

K kmol kJ T T

T C V ??+?+?+=---

（9）

文献[5]中给出的计算公式为：

)/(108644.3109874.12571.17263K kmol kJ T T C V ??-?+-- （10）

文献[6]中给出的计算公式为：

)

/(10

10

02748.010

181.010

4341.08056.03

3

9

2

6

3

K kg kJ T T

T C V ???+?-?+=----（11）

由此看出,氧气的定容比热容的计算公式之间存在差异，其计算结果和误差也各不相同（见表1），这种情况在其它气体中也普遍存在。

表1 o 2的计算定容比热容Cv 值与查表Cv 值kJ/(kmol ·K)

由表1可看出:按计算公式计算出的定容比热容值与查表值相比，其数值普遍相差较大。文献[5]提供的计算公式的最大误差为12. 26%，平均误差9.92%，且计算值均大于查表值，表明该计算公式不能很好地从数理统计方面对原始数据进行合理的描述；文献[4]和[6]的最大误差竟分别高达129.9%和45.44%，平均误差分别为63.81%和44.61%，疑是公式引用有误(例如公式中“-”误为“+”)，但因为文献[4]和[6]不仅没有勘误表，也未注明所引用公式的出处，还不能下此结论。考虑到上述文献所具有较大的影响力，为避免以讹传讹并提高回归计算公式的数理统计效果，有必要对理想气体的定容比热容重新进行数据回归分析，以得到更精确的数值及计算公式。

3.2 数据回归分析

理想气体的定容比热容为温度的单值函数，这种函数关系既可能是线性的，也可能是非线性的[7]。本文利用MA TLAB尝试建立

C与温度T的线性关系。以氧气为例，

V

用MATLAB6.5编程进行回归得：

b= 0. 6267 0. 0001

bint= 0. 61710. 6363

0. 00010. 0001

stats=0.9855746. 14330. 000

s=0.0057

图3

o

2

散点分布图 图4

o

2

残差分布图

由以上结果可见，相关系数：

9855.02

=R

（12）

则R=0. 9927>0. 9，统计量：

4

10

,1433.746-<=P F （13）

则此线性回归模型成立,回归公式为：V C =0. 6267+0. 0001TkJ/(kg ·K)，剩余标准差为0. 0057。回归公式的计算值与查表值的对比见表2。可以看出，计算值与查表值的最大误差为2. 91%，平均误差(绝对值)1. 93%，明显优于上述文献提供的计算公式。

同理可得其它常见理想气体的定容比热容的回归公式，见表3。

表2

o

2

的回归公式计算的定容比热容V C 值与查表V C 值kJ/(kmol ·K)

由上述回归分析结果看，最小的相关系数为0. 9781[9]，远大于0. 9，剩余标准差

也接近于0，从数理统计的角度看，采用线性模型是完全合理的。同时，上述回归公式的计算值与相应的查表值吻合得很好，最大误差及平均误差都不大于5%。

定容比热容是气体的重要热力学性质之一，其合理的数值计算公式的提出，不仅为计算气体的定容比热容提供了方便，也为热力学的电算提供了方便。借助数值计算软件MATLAB 的回归分析功能，得到了几种常见理想气体定容比热容与温度的线性回归公式，数理统计分析的结果证明了这些回归公式的合理性，回归公式的计算误差很小，最大误差低于5%[10]，能很好地满足科研和工程技术中热力学计算的要求。

表3 几种理想气体定容比热容回归公式

4 Matlab 在处理化工热力学数据中的应用

4.1 问题的描述

常压下,DMC (碳酸二甲酯)和CH3OH (甲醇)二元物系的Wilson 活度模型是一个带二元参数1A 和2A 的超越方程。求得该二元参数的值,是计算二元物系相平衡的关键。我们将对所求问题经过适当的变形后直接给出最终的方程形式：

),,,()(j i j i i

c i e A A x x G y y = )

exp(

)()

exp(

i

j j j j j i i j

i i j j i s x A x A x x A x p x A x A x x p +++=

（14）

（其中i= 1或2,且i= 1时j= 2；i= 2时j= 1）

上式中i e y 、i c y 为汽相组分摩尔百分率的实验值和计算值，且1c y +2c y =1，2

1e e y y +=1。

i x ，j x 为液相组分摩尔百分率的实验值，且i x +j x =1；i s p 为不同组分的饱和蒸汽压，

为了避免叙述过多的热力学专业知识，在不影响分析问题的前提下,我们将把根据物性常数计算得到的i s p 值直接给出；实验压力恒为p=101.325kPa ，具体数据见表4。由于上式是非线性的，用一般消去法求解是无效的，所以本文将采用非线性最小二乘法来寻求最佳参数。

表4 待处理数据

注：表4中x 1

e , y 1

e 和p 是文献[5]中的实验数据，p 1

s ，p 2

s 是由物性常数计算得到的。

4.2 利用Matlab 处理实验数据的原理及方法

非线性最小二乘法的数学模型为:

)(min 1

2

x f m

i i ∑

==

（15）

其中),...,,(321n x x x x x =是n 维向量，,))(),...(),(()(21T m x f x f x f x f =则上式又可表示为：

2

2||)(||min x f （16）

求解非线性最小二乘问题的算法主要有Gauss-Newton 法和Levenberg-Marquadt 法。Gauss-Newton 法通过在每一次迭代步骤求解下列线性最小二乘问题来获得搜索方向k d

22

||)()(||min k k k x F d x J - （17） 搜索方向k d 可以用于一维搜索,以保证每次迭代都使f(x)减小。Levenberg-Marquadt 法(又称阻尼牛顿法)则是用下式求搜索方向k d 。

)()())(((k k k k T

k x F x J d I x J x J -=+λ （18）

其中k λ为阻尼因子。它可以控制k d 的大小和方向。当k λ=0时,即为Gauss-Newton 法；当k λ→∞时，即为最速下降法。

在Matlab 的优化工具箱中有求解非线性最小二乘问题的函数Lsqnonlin():[x,resnorm,residual,exitflag,output,lanbda,jacobian]=Lsqnonlin(‘fun’,x0,LB ，UB,options,p1,p2,…)。在这个函数中，输入变量中只有‘fun ’和0x 是不可缺省的，fun

是给出向量函数的m-文件的文件名，0x 是x 的初值。输出宗量中只有x 是不可缺省的，

它是问题的解。针对不同的问题，可以选择不同的算法来计算。默认时，该函数将采用Gauss-Newton 法进行一维搜索，在残差较小时，这种算法有较好的效果和较快的速度。也可通过将输入变量https://www.360docs.net/doc/8a4699538.html,rgescale 设置为‘off’来使用Levenberg-Marquardt 算法进行一维搜索[11]。下面我们就来运用这一函数来寻求最佳配偶参数1A 和2A ，并将得到的参数代入方程(1)来检验计算值和实验值的误差。 (1)建立以汽相组成误差的平方和作为目标函数:

])()[(2

222114

1

1

j c e j c j e y y y y

F -+-=

∑= （19）

其中1e y 和2e y 为实验值，1

c y 和2c y 为待计算值。

(2)建立jscs.m 文件:

%计算配偶参数

Function z=jscs(a)

x1=[0.000 00 0.009 8 0.049 5 0.081 1 0.102 00.149 7 0.249 4 0.321 3 0.450 1 0.690 50.911 2 0.965 6 0.981 8 0.999 9];

x2=1-x1;

y1e=[0.000 00 0.014 80 0.060 30 0.095 10 0.110 20 0. 145 40 0. 198 80 0.

227 10 0.262 50 0. 324 00 0. 555 10 0. 795 30 0.899 90 0.999 90];

y2e=1-y1e;

p=[101.325 101.325 101.325 101.325 101.325 101.325 101.325 101.325 101.325 101.325 101.325 101.325 101.325 101.325];

ps1=[43.264 9 42.641 7 41.720 9 41.417 6 41.266 7 41. 266 7 41. 569 0 42. 026 0 43.421 9 49. 747 0 71. 294 4 85. 881 8 92.128 2 101.240 1];

ps2=[101.746 2 100.153 1 97.802 5 97.029 1 96.644 4 96.644 4 97.415 2 98.580 9 102.147 7 118.409 4 174.799 0 213.639 8 230.401 8 254.975];

for i=1:14

y1c(i)=exp(a(1)*x2(i)/(x1(i)+a(1)*x2(i)))*

ps1(i)*x1(i)/(p(i)*(x1(i)+a(1)*x2(i))*exp((a(2)*x2(i)/(x2(i)+a(2)*x1(i)))));

y2c(i)=exp(a(2)*x1(i)/(x2(i)+a(2)*x1(i)))*

ps2(i)*x2(i)/(p(i)*(x2(i)+a(2)*x1(i))*exp((a(1)*x1(i)/(x1(i)+a(1)*x2(i)))));

z(i)=(y1e(i)-y1c(i))^2+(y2e(i)-y2c(i))^2; end

%计算平均偏差

y1c=(ps1.*x1.*exp(a(1).*x2./(x1+a(1).*x2)))./(p.*(x1+a(1).*x2).*exp(a (2).*x2./(x2+a(2).*x1)));

n=size(ps);

meandev=sum(1.0-y1c./y1e)/n(2) (3)在Matlab 命令窗口输入:

0a =[1 1]；

%0a 为参数1A ,2A 的初值 x=lsqnonlin(‘jscs’, 0a ) 运行后,具体结果见表5。

表5 计算结果

注：平均偏差=[∑-14

1

0.1(计算值÷实验值]]÷14

4.3 小结

汽相组分的计算值与实验值的误差能够满足实验的要求，达到了计算的目的，说明matlab 的计算程序是有效可行的。本文得到的计算值与文献的计算值相比误差更小，并且编程简单，体现了matlab 进行科学计算的优越性。

5 总结

由以上论述可以看出MATLAB 在热物理学方面有很大的用处，方便快捷清晰形象的解决了热物理学中一些复杂的问题。在实际应用中，更多的复杂问题都可以用MATLAB 方便快捷的解决，应该根据具体情况，灵活运用MATLAB 工具。