第26卷 第2期 2006年6月
西安科技大学学报Vol.26 No.2Jun.2006
JOURNALOFXI.ANUNIVERSITYOFSCIENCEANDTECHNOLOGY
文章编号:1672-9315(2006)02-0286-04
用Excel快速求解一维非稳态对流扩散方程
郑 宁1,李建伟2,褚维盘1
(1.西安科技大学基础课部;2.西安科技大学化学与化工系,陕西西安 710054)
y
摘 要:对流扩散方程在物理及工程问题中具有广泛的应用。文中针对一维非稳态对流扩散方程,分析了5种常用差分方法的截断误差和稳定条件,给出了一种利用Excel软件的自动迭代功能代替编程进行求解的新方法,对不同离散格式的稳定性条件及截断误差进行了分析。最后以套管换热器的物流能量交换问题为例,分别用修正中心显式格式和Crank-Nicolson隐式格式介绍了Excel迭代求解过程与步骤,并对二种方法的求解结果进行了比较。关键词:对流扩散方程;Excel;迭代;数值解中图分类号:TK124 文献标识码:A
Solvingonedimensiondiffusion2convectionequationbyExcel
ZHENGNing1,LIJian2wei2,CHUWei2pan1
(1.Dept.ofBasicCourses,Xi.anUniversityofScienceandTechnology;
2.Dept.ofChemistryandChemicalEngineering,Xi.anUniversityofScienceandTechnology,Xi.an710054,China)
Abstract:Diffusion2convectionequationshavebeenappliedinmanyfieldsofphysicsorengineering.Inthisarticle,thecut2offerrorandstableconditionoffivediscretizationformsforonedimensiondif2fusion2convectionequationareanalyzed,andanewmethodsolvingtheequationwithexplicitorim2plicitmethodbyExcelisillustratedbyaheatexchangingexample.Theresultdataiswellconsistentwiththosebyprogrammingsolving,butthemethodusedhereismuchsimpleandcanbeeasilyhan2dled.
Keywords:diffusion2convectionequation;Excel;iteration;approachsolving
0 前 言
对流扩散方程描述了物质传输及热传递的综合过程,在水利工程、环境工程及化工、冶金、航空等研究领域里受到了充分重视,因此对流扩散方程的数值求解一直是人们关注的问题之一。目前求解对流扩散方程的方法虽然较多,但就其本质来说,最终都要编程实现[1-2]
。文中应用具有较强数据分析与运算能力的电子表格软件
Excel求解常系数的一维非稳态对流扩散方程,实现过程简捷方便,避免了繁琐的编程过程。
1 一维非稳态对流扩散方程的常用差分法
考虑如下形式的对流扩散方程
y收稿日期:2005-06-10
基金项目:陕西省自然科学基金项目(2004E223)
宁(),,,.
第2期 郑 宁等 用Excel快速求解一维非稳态对流扩散方程
287
5u5u52u
+c=(1)5t其中|x|0,a,c为常数,a>0。如果给定初始条件和边界条件,(1)式就成为一维对流扩散方程混合条件的求解问题。根据目前常用的5种有限差分方法对方程进行离散化处理,并利用泰勒级数和VonNeumann条件分别推出不同方法的截断误差及稳定性条件,结果见表1列(下文同)。
表1 一维对流扩散方程的不同离散格式
Tab.1 Variousdiscretizationmethodsforonedimensiondiffusion-convectionequation
名 称中心显式格式
离散格式
1nnnnun+-ujnuj+unj1-uj-1j+1-2uj+uj-1
+c=aS2hh21n
un+-ununa-unn
+c=2(unj+1-2uj+uj-1)Shh1n
un+-ununjjj+1-uj-+c
1un+-unujn-unjjj-+c
Sh
1
nn
S2unj+1-2uj+uj-c)h1
[3-7]
。其中上标表示时间序列,下标表示坐标序
截断误差O(S+h2)
a
稳定条件[h2
2
1迎风差分格式O(S+h)(a+
hcS)[h2
修正中心显式格式
=(a+
O(S+h2)O(S+h2)
2a
2[12+chh2aRS
[h2
12
指数型格式
1
=-
h(1-e)
Kh
nnn
(uj+1-2uj+uj-1)
Crank-Nicolson格式
1nnn+1n+1
un+-uncuj+1-uj-1uj+1-uj-1a2n+1jj
+(+=D(uj+unj)2h2
O(S2+h2)
无条件稳定
从表1中可以看出,上述5种差分格式中,前4种为显式差分格式,可由前一时间层结点值直接迭代求解,而最后的Crank-Nicolson格式是一种隐式格式,要用到后一时间层的结点值循环迭代求解。一般地,显式差分格式虽然编程较隐式差分格式简单,
但由上表可见其精度与稳定性不如隐式格式。
2 用Excel求解的过程
Excel具有通过对单元格的引用结合鼠标的拖放功能实现数据或公式的迭代计算,这部分功能在产生等差、等比数列和批量数据的转换处理过程得到了广泛的应用[8],而在求解偏微分方程方面尚未见报道。这里分别以一种显式格式与一种隐式格式的求解为例说明Excel的实现过程。
套管换热器是化工过程中实现物流能量交换的常见设备。图1表示了内管中的反应介质与套管中的蒸汽动态传热过程。由于内管管径很小,故在计算中均不考虑径向热传导,在内管中沿轴向取一微元(x,x+dx)进行能量计算,可得管内物流的热传导方程如下
5T2KK52T5T
=(Tw-T)=+-u(2)
5trQCpQCp5x25t
式中 K为套管壁处的对流换热系数;Q为反应物流的密度;K为反应物流的导热系数;Cp为物流比热;u为物流流速;r为内管半径;Tw为套管壁温度(近似为蒸汽温度);T0为流体入口温度。
将这些物性参数[9]95代入可得
5T5T52T
+3=0.001+2(Tw-T)5x
Tw=150,Tj0[x[1
(0)
图1 套管式换热器的传热过程
Fig.1 Heattransportationinthedoublepipeexchanger
=30,T0
x=1
(n)
=30
(3)
|5t
=0
下面分别以修正中心显式格式和Crank-Nicolson格式为例说明利用Excel软件求解显式格式或隐式格式的
2.1 构造迭代公式,引入边界条件
1
取S=0.01,h=0.1,将微分方程按以上差分格式离散化处理,经移项化简得到关于Tn+的迭代公式。j
1nn
修正中心显式格式 Tn+=0.888Tjn-0.104Tj+j1+0.196Tj-1+31Crank-Nicolson格式 Tn+=j
nnn+1n+1
(0.979Tjn-0.0745Tj+1+0.0755Tj-1-0.0745Tj+1+0.0755Tj-1+3)1.001
以步长为间断值在Excel表格的第1列和第1行分别输入网格节点的坐标,并将初边值条件T(x,0)=30和
n
n
T(0,t)=30引入到Excel表格的第2列和第2行中,在x=1处根据边界条件有T10=T9
表2 区域网格划分及初、边值录入格式
Tab.2 Discreitizationofsolvingregionandboundaryorinitialdatainput
00.010.02...tntn+I...
303030...3030...
Tnj-1
0.130
0.230
......
xj-130
xj30
xj+130
0.730
0.830
0.930
1=T9
Tnj
1Tn+j
Tnj-1
1Tn+j-1
1Tn+j-1
2.2 迭代求解
1nn
在C4中输入迭代公式,可以看出在修正中心差分格式中Tn+仅与前一时刻tn-1的Tnjj-1,Tj,Tj+1处的温度值有关,这些值均为已知,为显式格式,可直接进行迭代。在C4的迭代公式中可直接引用其所在的单元B3,C3,D3,并拖放至K4,然后选定C4到L4向下拉,即可得到不同时刻在x轴向上的温度分布(图2)。
1+1n+1
Crank-Nicolson差分格式中的Tn+除了与上述已知的温度值有关外,还与Tnjj-1或Tj+1(未知)有关,为五
点隐式格式,此时在迭代传递时Excel会出现/循环引用0警告,需要通过内钳的循环迭代计算实现。为避免循环引用造成的错误,打开菜单栏中的工具/选项中的重新计算标签页面,选中迭代计算,并输入最大迭代次数I或最大容许误差值E
。
图2 修正中心格式的求解结果图3 Crank2Nicolson隐式格式的求解式的求解结果
Fig.2 ResultsbymodifiedcentraldifferencemethodFig.3 ResultsbyCrank2Nicolsonimplicitmethod
按相同的方法在C4单元格中输入迭代公式,并对L列中各行的值按自然边界条件离散引用K列中的值。对于隐式格式,Excel须经过多次迭代才能得到数值解,开始时自动将Tnj+1的引用单元格置零而参与迭代,所有单max1n(i-1
|Tn+|
I,迭代终止。输入格式及结果如图3所示,通过比较可以发现,两种格式的迭代结果非常接近,最大误差
e,且随着时间的推移,误差呈减小趋势,与用各自编程求解的结果一致[9]97。
3 结果与讨论
提出了一种运用Excel迭代功能实现一维对流扩散方程的方法。它具有简单快捷、易于实现的特点,并且能方便地以图形格式输出结果。以套管换热问题为例,针对显、隐式格式分别介绍了它的求解步骤与实现过程。结果表明,两种方法得到的结果与实验数据及编程运算结果都非常接近。
参考文献:
[1] 李立康,於崇华,朱政华.微分方程数值解法[M].上海:复旦大学出版社,2003.[2] FerzigerJH,PericM.Computationalmethodsforfluiddynamics[M].Berlin:Springer,1996.[3] 陆金甫,顾丽珍,陈景良.偏微分方程差分方法[M].北京:高等教育出版社,1988.
[4] 王同科.一维非线性对流占优扩散方程的变网格特征差分方法[J].计算物理,2003,20(6):18-21.[5] 秦新强,李秋芳.对流方程的一种特征差分算法[J].西安理工大学学报,2001,17(4):32-35.
[6] 由同顺.非线性对流-扩散方程的高阶特征)))差分格式及其误差估计[J].数值计算与计算机应用,1994,15(4):321
-317.
[7] 阴继翔,李国君,李卫华.对流扩散方程不同格式的数值稳定性分析[J].太原理工大学学报,2004,35(2):121-125.[8] 张学辉,张桂月.迭代算法在Excel中的实现[J].计算机时代,2005,(1):13-17.[9] 方利国,陈 砺.计算机在化学化工中的应用[M].北京:化学工业出版社,2003:95-97.
(上接第195页)
[3] 汤友谊.对煤体结构形态及成因分类的改进和完善[J].焦作工学院学报(自然科学版),2004,23(3):161-164.[4] 龙王寅.潘一矿8煤层煤体结构特征及煤与瓦斯突出区域性预测[J].安徽地质,2003,13(1):70-73.
[5] 张国成.淮南矿区井田小构造对煤与瓦斯突出的控制作用[J].焦作工学院学报(自然科学版),2003,22(5):329-333.[6] 薛喜成.重庆市砚石台煤矿煤层流变特征及规律探讨[J].煤田地质与勘探,2003,31(5):20-23.[7] 王桂梁.论煤层流变[J].中国矿业大学学报,1988,17(3):16-25.
[8] 李 涛.煤岩流变分析[C].//煤炭科学研究总院西安分院.煤炭科学研究院地质勘探分院文集.西安:陕西人民出版
社,1987:81-85.
[9] 龙王寅,朱文伟,徐静,等.利用测井曲线判识煤体结构探讨[J].中国煤田地质,1999,11(3):64-69.
(上接第199页)参考文献:
[1] 林年丰,李静昌,钟佐棠.环境水文地质学[M].北京:地质出版社,1990.
[2] 阎碟瑞.利用脱硫微生物处理高浓度SO2-4制药废水的研究[C].//中国水文地质工程地质勘查院.环境地质研究.北
京:地震出版社,1993.
[3] 王兰生.山区城市地质环境演化中包气带的二次污染机制[J].水文地质与工程地质,1998,25(5):11.[4] 郭建伟.中原油田洒落石油对地面生态的污染与防治[J].地质技术经济管理,2004,(5):34.[5] 河南省地矿局.河南省商丘地区浅层地下水资源评价[M].北京:地质出版社,1993.[6] 朱庆陛,王世杰.供水水文地质手册[M].北京:地质出版社,1990.
第26卷 第2期 2006年6月
西安科技大学学报Vol.26 No.2Jun.2006
JOURNALOFXI.ANUNIVERSITYOFSCIENCEANDTECHNOLOGY
文章编号:1672-9315(2006)02-0286-04
用Excel快速求解一维非稳态对流扩散方程
郑 宁1,李建伟2,褚维盘1
(1.西安科技大学基础课部;2.西安科技大学化学与化工系,陕西西安 710054)
y
摘 要:对流扩散方程在物理及工程问题中具有广泛的应用。文中针对一维非稳态对流扩散方程,分析了5种常用差分方法的截断误差和稳定条件,给出了一种利用Excel软件的自动迭代功能代替编程进行求解的新方法,对不同离散格式的稳定性条件及截断误差进行了分析。最后以套管换热器的物流能量交换问题为例,分别用修正中心显式格式和Crank-Nicolson隐式格式介绍了Excel迭代求解过程与步骤,并对二种方法的求解结果进行了比较。关键词:对流扩散方程;Excel;迭代;数值解中图分类号:TK124 文献标识码:A
Solvingonedimensiondiffusion2convectionequationbyExcel
ZHENGNing1,LIJian2wei2,CHUWei2pan1
(1.Dept.ofBasicCourses,Xi.anUniversityofScienceandTechnology;
2.Dept.ofChemistryandChemicalEngineering,Xi.anUniversityofScienceandTechnology,Xi.an710054,China)
Abstract:Diffusion2convectionequationshavebeenappliedinmanyfieldsofphysicsorengineering.Inthisarticle,thecut2offerrorandstableconditionoffivediscretizationformsforonedimensiondif2fusion2convectionequationareanalyzed,andanewmethodsolvingtheequationwithexplicitorim2plicitmethodbyExcelisillustratedbyaheatexchangingexample.Theresultdataiswellconsistentwiththosebyprogrammingsolving,butthemethodusedhereismuchsimpleandcanbeeasilyhan2dled.
Keywords:diffusion2convectionequation;Excel;iteration;approachsolving
0 前 言
对流扩散方程描述了物质传输及热传递的综合过程,在水利工程、环境工程及化工、冶金、航空等研究领域里受到了充分重视,因此对流扩散方程的数值求解一直是人们关注的问题之一。目前求解对流扩散方程的方法虽然较多,但就其本质来说,最终都要编程实现[1-2]
。文中应用具有较强数据分析与运算能力的电子表格软件
Excel求解常系数的一维非稳态对流扩散方程,实现过程简捷方便,避免了繁琐的编程过程。
1 一维非稳态对流扩散方程的常用差分法
考虑如下形式的对流扩散方程
y收稿日期:2005-06-10
基金项目:陕西省自然科学基金项目(2004E223)
宁(),,,.
第2期 郑 宁等 用Excel快速求解一维非稳态对流扩散方程
287
5u5u52u
+c=(1)5t其中|x|0,a,c为常数,a>0。如果给定初始条件和边界条件,(1)式就成为一维对流扩散方程混合条件的求解问题。根据目前常用的5种有限差分方法对方程进行离散化处理,并利用泰勒级数和VonNeumann条件分别推出不同方法的截断误差及稳定性条件,结果见表1列(下文同)。
表1 一维对流扩散方程的不同离散格式
Tab.1 Variousdiscretizationmethodsforonedimensiondiffusion-convectionequation
名 称中心显式格式
离散格式
1nnnnun+-ujnuj+unj1-uj-1j+1-2uj+uj-1
+c=aS2hh21n
un+-ununa-unn
+c=2(unj+1-2uj+uj-1)Shh1n
un+-ununjjj+1-uj-+c
1un+-unujn-unjjj-+c
Sh
1
nn
S2unj+1-2uj+uj-c)h1
[3-7]
。其中上标表示时间序列,下标表示坐标序
截断误差O(S+h2)
a
稳定条件[h2
2
1迎风差分格式O(S+h)(a+
hcS)[h2
修正中心显式格式
=(a+
O(S+h2)O(S+h2)
2a
2[12+chh2aRS
[h2
12
指数型格式
1
=-
h(1-e)
Kh
nnn
(uj+1-2uj+uj-1)
Crank-Nicolson格式
1nnn+1n+1
un+-uncuj+1-uj-1uj+1-uj-1a2n+1jj
+(+=D(uj+unj)2h2
O(S2+h2)
无条件稳定
从表1中可以看出,上述5种差分格式中,前4种为显式差分格式,可由前一时间层结点值直接迭代求解,而最后的Crank-Nicolson格式是一种隐式格式,要用到后一时间层的结点值循环迭代求解。一般地,显式差分格式虽然编程较隐式差分格式简单,
但由上表可见其精度与稳定性不如隐式格式。
2 用Excel求解的过程
Excel具有通过对单元格的引用结合鼠标的拖放功能实现数据或公式的迭代计算,这部分功能在产生等差、等比数列和批量数据的转换处理过程得到了广泛的应用[8],而在求解偏微分方程方面尚未见报道。这里分别以一种显式格式与一种隐式格式的求解为例说明Excel的实现过程。
套管换热器是化工过程中实现物流能量交换的常见设备。图1表示了内管中的反应介质与套管中的蒸汽动态传热过程。由于内管管径很小,故在计算中均不考虑径向热传导,在内管中沿轴向取一微元(x,x+dx)进行能量计算,可得管内物流的热传导方程如下
5T2KK52T5T
=(Tw-T)=+-u(2)
5trQCpQCp5x25t
式中 K为套管壁处的对流换热系数;Q为反应物流的密度;K为反应物流的导热系数;Cp为物流比热;u为物流流速;r为内管半径;Tw为套管壁温度(近似为蒸汽温度);T0为流体入口温度。
将这些物性参数[9]95代入可得
5T5T52T
+3=0.001+2(Tw-T)5x
Tw=150,Tj0[x[1
(0)
图1 套管式换热器的传热过程
Fig.1 Heattransportationinthedoublepipeexchanger
=30,T0
x=1
(n)
=30
(3)
|5t
=0
下面分别以修正中心显式格式和Crank-Nicolson格式为例说明利用Excel软件求解显式格式或隐式格式的
2.1 构造迭代公式,引入边界条件
1
取S=0.01,h=0.1,将微分方程按以上差分格式离散化处理,经移项化简得到关于Tn+的迭代公式。j
1nn
修正中心显式格式 Tn+=0.888Tjn-0.104Tj+j1+0.196Tj-1+31Crank-Nicolson格式 Tn+=j
nnn+1n+1
(0.979Tjn-0.0745Tj+1+0.0755Tj-1-0.0745Tj+1+0.0755Tj-1+3)1.001
以步长为间断值在Excel表格的第1列和第1行分别输入网格节点的坐标,并将初边值条件T(x,0)=30和
n
n
T(0,t)=30引入到Excel表格的第2列和第2行中,在x=1处根据边界条件有T10=T9
表2 区域网格划分及初、边值录入格式
Tab.2 Discreitizationofsolvingregionandboundaryorinitialdatainput
00.010.02...tntn+I...
303030...3030...
Tnj-1
0.130
0.230
......
xj-130
xj30
xj+130
0.730
0.830
0.930
1=T9
Tnj
1Tn+j
Tnj-1
1Tn+j-1
1Tn+j-1
2.2 迭代求解
1nn
在C4中输入迭代公式,可以看出在修正中心差分格式中Tn+仅与前一时刻tn-1的Tnjj-1,Tj,Tj+1处的温度值有关,这些值均为已知,为显式格式,可直接进行迭代。在C4的迭代公式中可直接引用其所在的单元B3,C3,D3,并拖放至K4,然后选定C4到L4向下拉,即可得到不同时刻在x轴向上的温度分布(图2)。
1+1n+1
Crank-Nicolson差分格式中的Tn+除了与上述已知的温度值有关外,还与Tnjj-1或Tj+1(未知)有关,为五
点隐式格式,此时在迭代传递时Excel会出现/循环引用0警告,需要通过内钳的循环迭代计算实现。为避免循环引用造成的错误,打开菜单栏中的工具/选项中的重新计算标签页面,选中迭代计算,并输入最大迭代次数I或最大容许误差值E
。
图2 修正中心格式的求解结果图3 Crank2Nicolson隐式格式的求解式的求解结果
Fig.2 ResultsbymodifiedcentraldifferencemethodFig.3 ResultsbyCrank2Nicolsonimplicitmethod
按相同的方法在C4单元格中输入迭代公式,并对L列中各行的值按自然边界条件离散引用K列中的值。对于隐式格式,Excel须经过多次迭代才能得到数值解,开始时自动将Tnj+1的引用单元格置零而参与迭代,所有单max1n(i-1
|Tn+|
I,迭代终止。输入格式及结果如图3所示,通过比较可以发现,两种格式的迭代结果非常接近,最大误差
e,且随着时间的推移,误差呈减小趋势,与用各自编程求解的结果一致[9]97。
3 结果与讨论
提出了一种运用Excel迭代功能实现一维对流扩散方程的方法。它具有简单快捷、易于实现的特点,并且能方便地以图形格式输出结果。以套管换热问题为例,针对显、隐式格式分别介绍了它的求解步骤与实现过程。结果表明,两种方法得到的结果与实验数据及编程运算结果都非常接近。
参考文献:
[1] 李立康,於崇华,朱政华.微分方程数值解法[M].上海:复旦大学出版社,2003.[2] FerzigerJH,PericM.Computationalmethodsforfluiddynamics[M].Berlin:Springer,1996.[3] 陆金甫,顾丽珍,陈景良.偏微分方程差分方法[M].北京:高等教育出版社,1988.
[4] 王同科.一维非线性对流占优扩散方程的变网格特征差分方法[J].计算物理,2003,20(6):18-21.[5] 秦新强,李秋芳.对流方程的一种特征差分算法[J].西安理工大学学报,2001,17(4):32-35.
[6] 由同顺.非线性对流-扩散方程的高阶特征)))差分格式及其误差估计[J].数值计算与计算机应用,1994,15(4):321
-317.
[7] 阴继翔,李国君,李卫华.对流扩散方程不同格式的数值稳定性分析[J].太原理工大学学报,2004,35(2):121-125.[8] 张学辉,张桂月.迭代算法在Excel中的实现[J].计算机时代,2005,(1):13-17.[9] 方利国,陈 砺.计算机在化学化工中的应用[M].北京:化学工业出版社,2003:95-97.
(上接第195页)
[3] 汤友谊.对煤体结构形态及成因分类的改进和完善[J].焦作工学院学报(自然科学版),2004,23(3):161-164.[4] 龙王寅.潘一矿8煤层煤体结构特征及煤与瓦斯突出区域性预测[J].安徽地质,2003,13(1):70-73.
[5] 张国成.淮南矿区井田小构造对煤与瓦斯突出的控制作用[J].焦作工学院学报(自然科学版),2003,22(5):329-333.[6] 薛喜成.重庆市砚石台煤矿煤层流变特征及规律探讨[J].煤田地质与勘探,2003,31(5):20-23.[7] 王桂梁.论煤层流变[J].中国矿业大学学报,1988,17(3):16-25.
[8] 李 涛.煤岩流变分析[C].//煤炭科学研究总院西安分院.煤炭科学研究院地质勘探分院文集.西安:陕西人民出版
社,1987:81-85.
[9] 龙王寅,朱文伟,徐静,等.利用测井曲线判识煤体结构探讨[J].中国煤田地质,1999,11(3):64-69.
(上接第199页)参考文献:
[1] 林年丰,李静昌,钟佐棠.环境水文地质学[M].北京:地质出版社,1990.
[2] 阎碟瑞.利用脱硫微生物处理高浓度SO2-4制药废水的研究[C].//中国水文地质工程地质勘查院.环境地质研究.北
京:地震出版社,1993.
[3] 王兰生.山区城市地质环境演化中包气带的二次污染机制[J].水文地质与工程地质,1998,25(5):11.[4] 郭建伟.中原油田洒落石油对地面生态的污染与防治[J].地质技术经济管理,2004,(5):34.[5] 河南省地矿局.河南省商丘地区浅层地下水资源评价[M].北京:地质出版社,1993.[6] 朱庆陛,王世杰.供水水文地质手册[M].北京:地质出版社,1990.