放射性气体扩散的预测
摘要
摘要中要把文章中模型的方法、思想、技巧、结论体现出来。
日本福岛核电遭遇日本近海发生9.0级地震并引发了大海啸的破坏，发生了核泄漏，核事故发生后，气载放射性污染物经大气扩散，在短期内对环境产生大范围的影响，日本核污染扩散问题引起了国际社会的广泛关注。进行核事故应急响应决策时，气载放射性污染物大气扩散过程的模拟非常重要。本文结合高斯烟羽模型、线性拟合，以及微分方程模型，运用MATLAB软件，分析泄漏源强度、风速、大气稳定度参数等因素对放射性气体扩散的影响，预测了放射性气体浓度在不同时间，不同地区的浓度变化，并且本文模型中数据可以根据不同的实际情况而加以改变，因而是本文的应用范围大大增加，可以适用于具有较强的应用型。对于问题一，讨论在无风的情况下，放射性气体以s m/s的匀速在大气中向四周扩散。本问中由于不考虑风力的影响，且扩散出来的气体匀速向四周散开，这样经过任意时刻t，扩散的气体围成一个半径为st的球，且距球心位置不同的地方浓度值不同。采用列数列的表现方法，设定相同时间段t，把条件进行整理，并经过简单计算得出每段时间所预测得到的扩散距离r和浓度C。利用MATLAB 软件对数据进行线性拟合，采用微分方程模型得到核电站周边放射性气体在不同地区，不同时间段的浓度变化，得出随着离泄漏源距离的延伸，最后放射性物质的浓度越来越小，趋近于零，即当L趋向无穷时，C(x,y,z,t)趋向于零；当时间趋向于无穷时，C(x,y,z,t)也趋于无穷。
对于问题二，要探究风速对放射性物质浓度分布的影响。风速的处理是问题的中心，采用大气污染的经典高斯扩散模型,实现了高斯烟团气体扩散模型的动态预测,分析计算了气体扩散过程中的各关键参数。
对于问题三，本文在问题二的基础上，结合考虑风速和放射性物质扩散速度在空间中的矢量运算，将在上风和在下风不同情况下与传播速度s之间的比较的分析，利用高斯烟羽模型对核电站周边地区的浓度进行预测，然后，利用MATLAB 软件，将相关数据代入程序，我们得到核电站周边地区的浓度分布的等高曲线。对于问题四，本文参阅整理大量气象、地理、新闻资料，选择我国东海岸典型地域---山东半岛作为研究对象，综合考虑对应海域平均风速及风向、地理距离、海水对放射性物质扩散的部分反射系数等因素，集合核电站周边的浓度等高线，可
以直观地看出日本福岛核电站对我国东海岸以及美国斯海岸的影响。
关键词：放射性气体扩散浓度变化高斯修正模型预测
一、问题重述
日本福岛核电遭遇日本近海发生9.0级地震并引发了大海啸的破坏，发生了核泄漏，释放出大量具有放射性的物质，日本核污染扩散问题引起了国际社会的广泛关注。因此，正确的测出大气中放射性物质的浓度在环境监测以及安全评估中具有重要意义。
设核电站发生泄漏，放射性气体匀速以浓度为0p排出，速度为m kg/s，在无风的情况下，匀速在大气中向四周扩散, 速度为s m/s。
问题一，建立一个描述核电站周边不同距离地区、不同时段放射性物质浓度的预测模型。
问题二，当风速为k m/s时，预测核电站周边放射性物质浓度的变化情况。
问题三，当风速为k m/s时，分别预测出上风和下风L公里处，放射性物质浓度的预测模型。
问题四，将建立的模型应用于福岛核电站的泄漏，计算出福岛核电站的泄漏对日本本土、我国东海岸、及美国西海岸的影响。
二、问题分析
对于问题一，讨论在无风的情况下，放射性气体以s m/s的匀速在大气中向四周扩散。核电站源源不断泄漏引起气体扩散船舶可以看做在无穷空间有连续点源导致的扩散过程，能够由二阶抛物型偏微分方程描述放射性气体扩散过程中浓度变化规律。本问中由于不考虑风力的影响，且扩散出来的气体匀速向四周散开，这
样经过任意时刻t，扩散的气体围成一个半径为st的球，且距球心位置不同的地方浓度值不同。利用MATLAB软件对数据进行线性拟合，采用微分方程模型得到核电站周边放射性气体在不同地区，不同时间段的浓度变化。
问题二，针对在有风速k m/s时，选取高斯连续点源扩散模型进行分析，这也是在所有有风速时，放射性物质扩散模型中最常见也最方便的一种泄漏方式。
不考虑垂直速度，并且假设空间放射性气体云的浓度服从高斯分布的情况下，运用高斯模型可以较合理的计算出核电站周边地区的放射性气体浓度。
对于问题三，本文在问题二的基础上，结合考虑风速和放射性物质扩散速度在空间中的矢量运算，将在上风和在下风不同情况下与传播速度s之间的比较的分析，得出必须满足风速大于自然扩散速度，放射性物质才能到达上风口。
对于问题四，本文参阅整理大量气象、地理、新闻资料，选择我国东海岸典型地域---山东半岛作为研究对象，综合考虑对应海域平均风速及风向、地理距离、海水对放射性物质扩散的部分反射系数等因素进行预测。在以上分析中，我们可知，通过收集福岛核电站和我国东岸以及美国西岸之间的距离，以及和两地之间的风速，风速数据，可以大致判断出福岛核电站的泄漏事故对我国东岸以及美国西岸的影响，然后进行仿真模拟得出比较准确的结果。
三、基本假设
考虑到放射性气体扩散的复杂性，为简单起见，在讨论扩散模型时都作了如下假设:
（1）瞬时泄漏假定瞬时完成，连续泄漏假定泄漏速率恒定；
（2）不考虑地面的地形、地物对气体扩散的影响，即气云在平整、无障碍物的地面上空扩散；
（3）为水平风向，且在气体扩散的过程中保持不变
（4）气体的扩散看作空中某已连续点向四周强度的释放气体，放射性气体在无（5）穷空间的扩散过程中不发生性质变化
（6）气体的传播服从扩散规律，级单位时间通过单位法相面积的流量与它的浓度梯度成正比
（7）定常态，及所有的变量不随时间变化
（8）假定释放的气体的密度与空气像差不多（不考虑重力与浮力的作用），且气体扩散过程没有发生化学反应 （9）扩散气体的性质与空气相同 （10）本文引用数据、资料均真实可靠；
四、 符号说明
4.1. 符号说明 符号
说明
s
无风时，放射性气体的传播速度 k 风速，单位/m s
H
泄漏点O 距有效地面的高度 t
任意扩散时刻
(,,,)C x y z t
空间任意一点的放射性物质浓度 0ρ
放射性气体排出时的浓度 m
放射性气体排出时的速度
(,,)i i x y z δ= 空间任意一点的放射性物质的扩散系数
Ω
空间域 V 空间域其体积 S 一规则的球面面积 1Q 在(,)t t t +?内通过Ω的流量
2Q Ω内放射性物质的增量
0Q
从泄漏源泄漏的放射性物质的总量 h ? 附加高度
s T 核泄漏出口处的温度 0T
环境温度
α
设地面反射系数 Q
源强，单位为/kg s
,,x y z σσσ
分别为用浓度标准差表示的,,x y z 轴上的扩散参数 Vs
沉降速度，单位为/m s d W
地面干沉积率 ?
冲洗系数
0.5T
放射性核素的半衰期
五、 模型的建立与求解
5.1问题一模型的建立与求解 5.1.1 问题一模型的建立
核电站源源不断泄漏引起气体扩散船舶可以看做在无穷空间有连续点源导致的扩散过程，能够由二阶抛物型偏微分方程描述放射性气体扩散过程中浓度变化规律。本问中由于不考虑风力的影响，且扩散出来的气体匀速向四周散开，这样经过任意时刻t ，扩散的气体围成一个半径为st 的球，且距球心位置不同的地方浓度值不同。
讨论在无风的情况下，放射性气体以s m/s 的匀速在大气中向四周扩散。核电站源源不断泄漏引起气体扩散船舶可以看做在无穷空间有连续点源导致的扩散过程，能够由二阶抛物型偏微分方程描述放射性气体扩散过程中浓度变化规律。本问中由于不考虑风力的影响，且扩散出来的气体匀速向四周散开，这样经过任意时刻t ，扩散的气体围成一个半径为st 的球，且距球心位置不同的地方浓度值不同。利用MATLAB 软件对数据进行线性拟合，采用微分方程模型得到核电站周边放射性气体在不同地区，不同时间段的浓度变化。
将气体从泄露时刻记作t=0，泄漏点选为坐标原点，时刻t 无穷空间中任一点（x,y,z ）的气体浓度记为C(x,y,z,t).```
由质量守恒定律Q 1=？Q 2 ···
这结果表明，对于任意时刻t 烟雾浓度C 的等值是球面2222R z y x =++,并且随着球面半径R 的增加C 的只是连续减少的；当R →∞或t →∞时，C(x,y,z,t) →0
5.1.2 问题一模型的求解
···我们通过将
5.2问题二模型的建立与求解 5,.2.1 问题二模型的建立 (1)泄漏源为近地面 （2）泄漏源为高架点
设泄漏源有效高度为H ，取其在地面投射为坐标原点，x 轴指向风向，考虑地面放射作用，可以得到烟羽模型的浓度分布计算公式为：
（3）扩散参数的确定
应用高斯模型的关键是确定扩散参数，··· ···
在进行和计算时，首先从表1-2中确定大气稳定度，我国在标准GB/T13201-1991中规定，当确定稳定度级别后，实际的扩散参数选择见中国有关环境评价标准采用的系数值-P-G 扩散曲线幂函数数据。
5,.2.2 问题二模型的求解
？？为上风向浓度分布，？？为下风向浓度分布。由两图我们可以直观地看出核电站周边地区的放射性气体的浓度变化规律。该程序中的源强、风速、大气稳定度参数、地面粗糙度参数和计算精确度等都可根据实际情况分析需要设置。
···
5.2.2.1模型参数的确定
在事故过程模型中所需参数的选取及确定十分重要，通常情况下气象参数的选取是利用该地区多年气象资料，采取工业安全与环保统计的方法进行有关参数的确定，而其他扩散参数是以实际测定为基础的。
a:大气稳定度的计算
根据国家标准(GB/ T 13201 -- 1991) )制定地方大气污染物排放标准的技术方法的规定，划分大气稳定度的级别，共分为6级A ~ F，A为极不稳定；F为极稳定。首先，根据释放源所在地的经度和纬度以及泄漏的日期和时间计算当时的太阳高度角h0；然后，由太阳高度角h0和云量查出太阳辐射等级；最后，再根据地面风速确定当时的大气稳定度，计算细节可参考文献。
表 1 大气稳定度的级别参考表
地面风速
1 () m s-
?白天太阳辐射阴天的白天或夜
间
有云的夜晚
强中弱薄云遮天或低云≧
0.5
云量≦
0.4
<3 A A -B B D - -
2-3 B B C D E F
3-5 B-
C
B-C C D D E
5-6 C-
D
C-D D D D D >6 C D D D D D b：扩散参数的计算
有风时的扩散参数，(,)
y z
σσ的确定采用Briggs给出一套扩散参数，如表1和
表2所示。
表2 Briggs 扩散参数（开阔平原田野）
大气稳定度
y σ
z σ
B 21)0001.01(16.0x x + x 12.0
C 21)0001.01(11.0x x +
21
)0002.01(08.0x x +
D 2
1
)0001.01(08.0x x +
21
)0015.01(06.0x x +
E 21
)0001.01(06.0x x + )0003.01(03.0x x + F
21
)0001.01(04.0x x +
)0003.01(016.0x x +
表3 Briggs 扩散参数（城市）
大气稳定度
x σ
y σ
C 2
1
)0004.01(22.0-+x x x 20.0
D 2
1
)0004.01(16.0-+x x
2
1
)0003.01(14.0-+x x E-F
2
1
)
004.01(11.0-+x x
2
1
)
0015.01(08.0-+x x
5.3.问题三模型的建立与求解 5.3.1 问题三模型的建立
本文在问题二的基础上，结合考虑风速和放射性物质扩散速度在空间中的矢量运算，将在上风和在下风不同情况下与传播速度s 之间的比较的分析，得出必须满足风速大于自然扩散速度，放射性物质才能到达上风口。
在考虑风速的条件下，由于我们在前面的假设中认为风向与x 轴相同，故我们要求计算上风和下风处的放射性物质浓度，只需在···，这样就可以求出当风速一定时，上风和下风1公里处，放射性物质浓度的估计模型。上风I 处浓度计
算式如下：···
下风1公里处浓度计算式如下：
···
根据扩散系数的定义可知，？=，？=，代入上两式可以求出上风和下风I处的浓度
5.3.2问题三模型的求解
我们通过对问题二的分析可以从图中的处在居上风处和下风初L公里时可以应付的得出相应点的气体浓度，方便直观。比如当L=100km时，通过观察上下风浓度等高线可以很直观，和方便地得出：在下风向的浓度是千分之四；在上风向时千分之五。而在实际中，则要根据不同实际情况，确定好各个参数，更准确地得出相应点的气体浓度。我们通过查询大亚湾核电站周围放射性物质的浓度数据用matlab进行仿真模拟，得出实际的图形与模拟的图形如下（？？），（？？），通过观察图形我们可以知道，两图拟合较好，说明以上所建模型的准确性较高。
···
5.4问题四的求解
从空气动力来说，我国东部的降水主要源于西太平洋的暖湿气团，在其影响下，广东所受影响要高于美国西海岸；从洋流来看，西太平洋盛行暖性洋流，沿我国东海岸经日本沿岸顺时针流向太平洋东部，美国西岸受影响较大。
在结合模型二、模型三的条件下，在参阅整理大量气象、地理、新闻资料，选择我国东海岸——山东半岛作为研究对象，综合考虑对应海域平均风速及风向、地理距离、海水对放射性物质扩散的部分反射系数等因素，并通过C++编程模拟计算，预测出放射性核物质将经过6.5天到达我国东海岸——山东半岛，且131I浓度值为：0.1003
?与，与实际情况比较吻合。
mBq m-
由于本文中，我们只需通过检测福岛的核泄漏是否对我国东岸以及美国西岸有影响，故我们只需在前面建立的模型的基础上，通过网络查询有关的泄漏在我国东
岸以及美国西岸的风向，风速以及之间的距离，带入建立的模型中，并用matlab 进行仿真模拟即可。
六、模型的分析
1.1.假设的合理性分析
1.2.灵敏度（稳定性、可靠性）分析（给出模型的适用范围）
1.3.原理误差分析（给出模型的误差范围）
1.4.……
七、模型的检验
用模型计算数据与实际数据进行比对，计算误差大小，结合模型的分析说明误差产生的原因，以及误差是否在模型估计和实际许可的范围之内。
八、模型的推广
模型在实际问题中的应用：可以用于解决哪一类的问题，用于不同类别问题时应对模型做出怎样的变化。主要阐述本文模型的广泛适用性。
九、模型的评价与优化
1.1.模型的优缺点分析
模型的优点
模型的缺点：（高斯烟团气体扩散模型气体覆盖范围在二维平面上的投影为椭圆形,在城市建筑物密度较高的地区,这显然是不太符合实际的,所以高斯烟团模型适合在建筑物较低矮或建筑物较少的平原或郊区。而化学工业的厂址一般都选择在郊区或人口密度较低的地区,说明高斯烟团气体扩散模型对这些地方有实值。通过试验可以发现,地理信息系）
1.2.模型的优化
1.2.1.模型的优化方案
1.2.2.优化模型的建立、求解与分析
参考文献：
[1]XXX，XXXXXXXXXXXXXXXXXXX，XXXXXXX，XXXXX；
[2]XXX，XXXXXXXXXXXXXXXXXXXXXXXXXXXXX，XXXXXXXXXXXX，XXXXX；
[3]XXX，XXXXXXXXXXXXXXXXXXX，XXXXXXX，XXXXX；
[4]XXX，XXXXXXXXXXXXXXXXXXXXXXXXXXXXX，XXXXXXXXXXXX，XXXXX。[编号] 作者，书名：起止页码，出版地：出版社，出版年。
[编号] 作者，论文名，杂志名，卷期号：起止页码，出版年。
[编号] 作者，资源标题，网址，访问时间（年月日）。
附件
1.附件一
2.附件二
徐州工程学院 课程设计报告 课程名称偏微分方程数值解 课题名称对流扩散方程 的迎风格式的推导和求解专业信息与计算科学 班级10信计3 姓名学号 指导教师杨扬 2013年 5 月23 日
一、实验目的： 进一步巩固理论学习的结果，学习双曲型对流扩散方程的迎风格式的构造 方法，以及稳定的条件。从而进一步了解差分求解偏微分方程的一些基本概念，掌握数值求解偏微分方程的基本过程。在此基础上考虑如何使用Matlab 的软件进行上机实现，并针对具体的题目给出相应的数值计算结果。 二、实验题目： ?? ? ??-=-==<<<<+=+);2/1exp(),1();exp(),0();2/exp()0,(10,10,11t t u t t u x x u t x f u b u a u xx x t 其中a1=1,b1=2, ) 2/exp(),(t x t x f --=。 用迎风格式求解双曲型对流扩散方程,观差分解对真解的敛散性()2/exp(t x u -= 三、实验原理： 1、用迎风格式求解双曲型对流扩散方程，迎风格式为： ) 01(21 1 )01(2112 1 1112 1 11 1<++-=-+->++-=-+--+++-+-+a f h u u u b h u u a u u a f h u u u b h u u a u u n j n j n j n j n j n j n j n j n j n j n j n j n j n j n j n j τ τ 若令,/*1,/*12h b h a r τμτ== 则迎风格式可整理为： > <<++-+-+=><>++++--=-+++-+2)01()()21(1)01()()21(111111a f u u r u r u a f u u r u r u n j n j n j n j n j n j n j n j n j n j τμμμτμμμ2、稳定条件： ） （） （01),*11*2/(01),*11*2/(2 2<-≤>+≤a h a b h a h a b h ττ（*） 四、数值实验的过程、相关程序及结果： 本次的实验题目所给出的边界条件是第一边界条件，直接利用所给的边界条件，我们可以给出界点处以及第0层的函数值，根据a1的正负性，使用相应的<1>或者<2>式，求出其他层的函数值。误差转化成图的形式，并输出最大值。 针对三种不同的输入对应输出结果 ：
质子交换膜燃料电池气体扩散层憎水性衰减机理研究 于书淳，李晓锦*，李进，邵志刚，衣宝廉 （中国科学院大连化学物理研究所，大连，辽宁，116023,Email:xjli@dicp.ac）作为质子交换膜燃料电池的重要组件之一，典型的双层气体扩散层由基底层和微孔层构成。其中，基底层通常由憎水处理过的碳纸构成，微孔层通常由碳粉和憎水剂构成。具有良好化学稳定性的聚四氟乙烯（PTFE）是气体扩散层中最常用的憎水剂。在燃料电池中，气体扩散层必须具有合适的憎水性能以实现良好的导气和排水功能[1]。然而，电池在长时间运行后，尤其是在较为苛刻的工作环境下（频繁的启动/停车、动态工况等），气体扩散层的憎水性会逐步变得下降[2]。陈等人采用恒电位氧化法对气体扩散层进行耐久性研究，发现氧化后微孔层表面的接触角显著下降[3]。Lee等人在研究气体扩散层耐久性的实验中也观察到同样的现象[4]。但是，文献中并没有对憎水性下降的原因进行深入的研究。 鉴于气体扩散层憎水性的下降会引起电极水淹并最终降低电池性能，有必要对憎水性下降的原因进行深入的研究，但是目前有关这方面的报道很少[5]。因此，我们的主要工作是在模拟的电池环境下考察气体扩散层憎水性下降的原因。实验分为恒电位氧化及酸浸泡两部分。恒电位氧化实验是以N2饱和的0.5M H2SO4为电解液，采用相对于饱和甘汞电极为1.25V的恒电位对气体扩散层进行氧化处理，酸浸泡实验是将气体扩散层浸泡在70℃、air饱和的1M H2SO4溶液中1200h.通过扫描电镜、红外光谱、X射线光电子显微镜、热重仪等手段对氧化前后的扩散层的特性进行分析。恒电位氧化实验结果发现，无论是对基底层还是整平层，氧化之后表面的形貌发生了改变；此外，碳的氧化不但导致亲水性氧化物的生成而且导致了碳材料及PTFE的流失。这也正是恒电位氧化条件下气体扩散层憎水性下降的原因。在酸浸泡实验条件下观察到同样的现象。 图1XPS全谱（a）氧化前整平层（b）氧化后整平层（c）氧化前基底层（d）氧化后基底层
1 扩散动力学方程——菲克定律 1.1 菲克第一定律 1.1.1宏观表达式 1858年，菲克（Fick ）参照了傅里叶（Fourier ）于1822年建立的导热方程，建立定量公式。 在t ?时间内，沿x 方向通过x 处截面所迁移的物质的量m ?与x 处的浓度梯度成正比： t A x C m ???∝ ? 即 )(x C D Adt dm ??-= 根据上式引入扩散通量概念，则有： x C D J ??-= （7-1） 图7-1 扩散过程中溶质原子的 分布
式（7-1）即菲克第一定律。 式中J 称为扩散通量，常用单位是mol /（)2s cm ?； x C ??浓度梯度； D 扩散系数，它表示单位浓度梯度下的通量，单位为2cm /s 或s m /2； 负号表示扩散方向与浓度梯度方向相反见图7-2。 1.1.2微观表达式 微观模型： 设任选的参考平面1、平面2上扩 散原子面密度分别为n 1和n 2，若n 1＝n 2，则无净扩散流。 假定原子在平衡位置的振动周期为τ，则一个原子单位时间内离开相对平衡位置跃迁次数的平均值，即跃迁频率Γ为 τ 1 = Γ (7-2) 由于每个坐标轴有正、负两个方向，所以向给定坐标轴正向跃迁的几率是Γ6 1。 设由平面l 向平面2的跳动原子通量为J 12，由平面2向平面1的跳动原 图7-2 溶质原子流动 的方向与浓度降低的方 向相一致 图7-3 一维扩散的微观 模型
子通量为J 21 Γ=11261n J (7-3) Γ=221 6 1 n J (7-4) 注意到正、反两个方向，则通过平面1沿x 方向的扩散通量为 ()21211216 1n n J J J -Γ=-= (7-5) 而浓度可表示为 δδn n C =??= 11 (7-6) 式(7-6)中的1表示取代单位面积计算，δ表示沿扩散方向的跳动距离（见图7-3），则由式(7-5)、式(7-6)得 ()dx dC D dx dC C C C C J -=Γ-=-Γ-=-Γ=2122116 1)(6 16 1δδδ (7-7) 式(7-7)即菲克第一定律的微观表达式，其中 26 1 δΓ=D (7-8) 式（7-8）反映了扩散系数与晶体结构微观参量之间的关系，是扩散系数的微观表达式。 三维情况下，对于各向同性材料（D 相同），则 C D x C k x C j x C i D J J J J z y x ??-=??+??+??-=++=)( （7-9） 式中：x k x j x i ?? +??+??=?为梯度算符。 对于各向异性材料，扩散系数D 为二阶张量，这时，
扩散系数计算 WTD standardization office【WTD 5AB- WTDK 08- WTD 2C】
扩散系数 费克定律中的扩散系数Ｄ代表单位浓度梯度下的扩散通量，它表达某个组分在介质中扩散的快慢，是物质的一种传递性质。 一、 气体中的扩散系数 气体中的扩散系数与系统、温度和压力有关，其量级为52 10/m s -。通常对于二元气体Ａ、B 的相互扩散，Ａ在B 中的扩散系数和B 在A 中的扩散系数相等，因此可略去下标而用同一符号Ｄ表示，即AB BA D D D ==。 表７－１给出了某些二元气体在常压下（5 1.01310Pa ?）的扩散系数。 对于二元气体扩散系数的估算，通常用较简单的由富勒（Fuller ）等提出的公式： 1/31/32 [()()]A B D P v v = +∑∑ （７－１９） 式中，D －Ａ、B 二元气体的扩散系数，2 /m s ； P －气体的总压，Pa ； T －气体的温度，Ｋ； A M 、 B M －组分Ａ、B 的摩尔质量，/kg kmol ； A v ∑、B v ∑－组分Ａ、B 分子扩散体积，3 /cm mol 。 一般有机化合物可按分子式由表７－２查相应的原子扩散体积加和得到，某些简单物质则在表７－２种直接列出。 5 1.01310Pa ?
式７－１９的相对误差一般小于１０％。 二、 液体中的扩散系数 由于液体中的分子要比气体中的分子密集得多，因此也体的扩散系数要比气体的 小得多，其量级为92 10/m s -。表７－３给出了某些溶质在液体溶剂中的扩散系数。 式估算： 15 0.6()7.410 T B AB A M T D V -φ=?μ 2/m s （７－２１） 式中，AB D －溶质Ａ在溶剂Ｂ中的扩散系数（也称无限稀释扩散系数），2 /m s ； T －溶液的温度，Ｋ； μ－溶剂Ｂ的粘度，.Pa s ；
一、扩散方程稳态扩散与非稳态扩散 1．稳态扩散下的菲克第一定律（一定时间内，浓度不随时间变化dc/dt=0） 单位时间内通过垂直于扩散方向的单位截面积的扩散物质流量（扩散通量）与该面积处的浓度梯度成正比 即J=－D（dc/dx） 其中D：扩散系数，cm2/s，J：扩散通量，g/cm2·s ，式中负号表明扩散通量的方向与浓度梯度方向相反。 可见，只要存在浓度梯度，就会引起原子的扩散。 x轴上两单位面积1和2，间距dx，面上原子浓度为C1、C2 则平面1到平面2上原子数n1=C1dx ，平面2到平面1上原子数n2=C2dx 若原子平均跳动频率f, dt时间内跳离平面1的原子数为n1f·dt 跳离平面2的原子数为n2fdt，但沿一个方向只有1/2的几率，则单位时间内两者的差值即扩散原子净流量。 令，则上式 2．扩散系数的测定：
其中一种方法可通过碳在γ-Fe中的扩散来测定纯Fe的空心园筒，心部通渗碳气氛，外部为脱碳气氛，在一定温度 下经过一定时间后，碳原子从内壁渗入，外壁渗出达到平衡，则为稳态扩散单位时单位面积中碳流量： A：圆筒总面积，r及L：园筒半径及长度，q：通过圆筒的碳量 则： 即： 则： q可通过炉内脱碳气体的增碳求得，再通过剥层法测出不同r处的碳含量，作出C-lnr曲线可求得D。 第一定律可用来处理扩散中浓度不因时间变化的问 3．菲克第二定律：解决溶质浓度随时间变化的情况，即dc/dt≠0 两个相距dx垂直x轴的平面组成的微体积，J1、J2为进入、流出两平面间的扩散通量，扩散中浓度
变化为，则单元体积中溶质积累速率为 （Fick第一定律） （Fick第一定律） （即第二个面的扩散通量为第一个面注入的溶质与在这一段距离内溶质浓度变化引起的扩散通量之和） 若D不随浓度变化，则 故： 4．Fick第二定律的解：很复杂，只给出两个较简单但常见问题的解 a. 无限大物体中的扩散 设：1）两根无限长A、B合?金棒，各截面浓度均匀，浓度C2>C1 2）两合金棒对焊，扩散方向为x方向 3）合金棒无限长，棒的两端浓度不受扩散影响 4）扩散系数D是与浓度无关的常数 根据上述条件可写出初始条件及边界条件 初始条件：t=0时, x>0则C=C1，x<0, C=C2 边界条件：t≥0时, x=∞,C=C1, x=－∞, C=C2
一个基于碳纸和碳布的气体扩散层对质子交换膜燃料电池 性能的影响 Sehkyu Park ?, Branko N. Popov电化学工程中心，化学工程系，南卡罗来纳大学，哥伦比亚，SC29208，USA 关键字：质子交换膜燃料电池；气体扩散层；复写纸；碳布；微孔层。 摘要：一个市售的基于碳纸和碳布的气体扩散层就如一个可以通过各种物理和电化学测量方法的大孔基板；压汞法，表面形态分析法，接触角测量法，水渗透测量法，偏振技术，和交流阻抗谱。和基于碳布的ELAT-LT-1400W相比，基于碳纸的SGL 10BB的双孔径分布和高水流动阻力是因为大孔基板不太透水，多疏水性和紧密的微孔层。当空气作为氧化剂时，用SGL10BB制作的膜-电极-组件表现出一种优越的燃料电池性能。交流阻抗响应表明一个具有大量的微孔和疏水性的微孔层更容易允许氧朝催化剂层扩散因为可以有效的除去水在催化剂层的气体流路。 1.序言 在质子交换膜燃料电池中，气体扩散层被嵌入催化剂层和气体流路之间。气体扩散层的主要功能：1，气体扩散；2，一个电流收集器；3，一个物理支持，从而确定催化剂的利用率和整体性能。它也允许水蒸气到膜和液态水从催化剂层出来。一个气体扩散层防湿透过可防止水倒流和提高反应物到催化活性位点。 一个气体扩散层包括一个大孔基板和一个有炭黑的微孔层。编制碳布或非编制碳纸由于其较高的透气性和电子导电性被广泛的运做大孔基板。一个微孔层可以减少催化剂层和大孔基板之间的欧姆电阻，在催化沉寂中提供非渗透性支持和管理液态水流动。 一个单层气体扩散层（如：碳纸和碳布）在燃料电池性能上的效果已经被几个研究人员研究，他们表明碳布可导致更高的性能主要由于较高的孔隙率和较低的水饱和度。此外，丰富的工作已经在进行研究这种微孔层的性能如何能像1，碳粉类型；2，碳载量（或微孔层厚度）和3，疏水剂的浓度在质子交换膜燃料电池中控制水的管理。然而，在大量文献中，大孔基板在气体扩散层关于孔隙特征的反应和产物运输中的作用还没有解决。我们在此项工作中的目标就是表征用碳纸或碳布制备的市售气体扩散层的物理属性和研究气体扩散层属性如何影响水的管理和氧气在质子交换膜燃料电池中流动的途径。 2.实验 2.1.气体扩散层的物理特性 多孔结构的气体扩散层被用一个压汞仪分析。为了进行分析，一小片的气体扩散层被称重并被加载上一个覆盖着金属箔的玻璃毛细管制成的样品杯，然后，在真空中从气体扩散层出气。之后，自动灌满水银。孔径分布曲线（PSD）从水银进入开始测定即水银的体积与贯穿孔所施加的压力。在所有的毛孔都是圆柱形的假设下，孔直径dp用一个众所周知的毛细管公式从P的值开始算：dp=4γcosθ/p (1)其中，γ和θ分别表示水银的表面张力和与样本中水银的接触角度
对流扩散方程的定解问题是指物质输运与分子扩散的物理过程和黏性流体流动的数学模型，它可以用来描述河流污染、大气污染、核污染中污染物质的分布，流体的流动和流体中热传导等众多物理现象。关于对流扩散方程的求解很也备受关注，因此寻找一种稳定实用的数值方法有着重要的理论与实际意义。 求解对流扩散方程的数值方法有多种，尤其是对流占优扩散方程，这些方法有迎风有限元法，有限体积法，特征有限体积法，特征有限差分法和特征有限元法，广义差分法，流线扩散法，以及这些方法与传统方法相结合的方法如迎风广义差分法，迎风有限体积法有限体积——有限元法等这些方法数值求解效果较好，及有效的避免了数值震荡，有减少了数值扩散，但是一般计算量偏大 近年，许多研究者进行了更加深入的研究，文献提出了对流扩散方程的特征混合元法，再次基础上，陈掌引入了特征间断混合元方法，还有一些学者将特征线和有限体积法相结合，提出了特征有限体积元方法（非线性和半线性），于此同时迎风有限元也得到较大的发展，胡建伟等人研究了对流扩散问题的Galerkin 部分迎风有限元方法和非线性对流扩散问题的迎风有限元，之后又有人对求解发展型对流扩散问题的迎风有限元法进行了理论分析 有限差分法和有限元是求解偏微分方程的常用数值方法，一般情况下考虑对流占优的扩散方程，当对流项其主导作用时，其解函数具有大梯度的过渡层和边界层，导致数值计算困难，采用一般的有限元或有限体积方法虽然具有形式上的高精度，不能解决数值震荡的问题，虽然我们不能简单的将对流占优扩散方程看做对流方程，但由于次方程中含有一阶不对称的导数，对流扩散方程仍会表现出“对流效应”，从而采用迎风格式逼近，尽量反应次迎风特点，此格式简单，克服了锋线前沿的数值震荡，计算结果稳定，之前的迎风格式只能达到一阶精度，我们采用高精度的广义迎风格式，此格式是守恒的，精度高，稳定性好，具有单调性，并且是特征线法的近似，有效的避免了锋线前沿的数值震荡。 有限体积是求解偏微分方程的新的离散技术，日益受到重视。有限体积与有限差分、有限元法最大的区别及优点在于有限体积将求解区域内的计算转化到控制体积边界上进行计算，而后二者均是直接（或间接）在域内计算，故有限体积有着明显的物理涵义，在很大程度上减少计算工作量又能满足计算精度要求，加快收敛速度。由于此方法讲散度的积分化为子域边界积分后子啊离散，数值解满足离散守恒，而且可以采用非结构网格，所以在计算物理特别是计算流体力学领域上有限体积有广阔的前景。 间断Galerkin（DG）方法是在1973年，Reed和Hill在求解种子迁移问题时，针对一阶双曲问题的物理特点提出的。之后C.Johnson，G.R.Richter等人对双曲问题的DG方法做了进一步的研究，并且得到了该机的误差分析结果，由于这种方法具有沿流线从“上游”到“下游”逐层逐单元计算的显示求解的特点，并且可以进行并行计算，所以被广泛应用于各类方程的求解。最近Douglas等人在{25}中处理二阶椭圆问题时，得到DG方法的有限元空间不需要满足任何连续性条件，因此空间构造简单，具有较好的局部性和并行性。DG发展的一个重要方面是对对流占优扩散方程的应用。G.R.Richter等在1992年提出利用DG方法求解定长对流扩散问题 近年DG方法有了新的发展，其中YeXiu提出间断体积元方法备受人们关注，2004年，她将有限体积法与DG相结合，提出了椭圆问题的间断有限体积法，此方法解除了逼近函数在跨越边界上连续的限制，之后更多的研究者应用到Stokes问题，抛物问题，双曲问题，并得到了较好的结果，该方法不但继承了有限体积的高精度计算简单及保持物理间局部守恒等优点，而且有限元空间无需满足任何连续性要求，空间构造简单，有较好的局部和并行性。 当对流扩散方程中的对流项占主导地位时，方程具有双曲方程的特点，这是由于对流扩散方程中的非对称的对流项所引起的迎风效应使对流扩散方程的数值求解更困难，用传统的中心差分法和标准的有限元求解会差生数值的震荡，从而使数值模拟失真，为了克服这一困难，早在20世纪50年代，就有人提出了迎风思想，由于使用迎风技巧可以有效的消除数值解不稳定性，因此吸引了众多学者的关注，从1977年，Tabata等人就针对对流扩散方程提出了三角形网格上的迎风格式{42,38}，并进行了深入的研究，梁栋基于广义差分法，提出并分析了一类建立在三角网格上的广义迎风差分格式，袁益让2001年就多层渗流方程组合系统提出并分析了迎风分数步长差分方法，以上均是讨论的线性对流扩散问题，胡建伟等通过引入质量集中算子，构造并分析了一类基于三角网格的质量集中型的部分有限元方法处理线性和非线性对流扩散问
国内氢燃料电池气体扩散层GDL研究进展实现国家燃料电池技术规模化的生产，就需要解决燃料电池里面核心材料国产化。气体扩散层GDL就是燃料电池电堆里面的一个核心材料，目前还是主要依靠进口来解决气体扩散层这样一个材料问题。气体扩散层总体来说，也是我们燃料电池技术发展里面一个卡脖子技术。国内南方科技大学氢能与燃料电池研究团队在过去两年多的发展里面，优化了燃料电池气体扩散层的结构，开发了气体扩散层的生产工艺，也实现了气体扩散层小批量生产。 首先从燃料电池单电池的结构开始，燃料电池电堆是由很多单电池串联组成的。对燃料电池的单电池有一个七层的结构，中间是质子交换膜，质子交换膜的两侧是催化层，催化层的两侧是有两片气体扩散层，再加上双极板就构成了燃料电池单电池这样一个结构。一个燃料电池单电池需要两片气体扩散层，虽然目前市场上气体扩散层的价格比质子交换膜要便宜一些，但是它的用量比较大，需要两片，所以它对燃料电池的价格有非常大的影响。美国能源部DOE按照现有的燃料电池技术所进行测算，气体扩散层在燃料电池电堆里面价格的分布，在小批量的时候，基本上价格占到21%，比质子交换膜还要贵一些。随着量的增加，气体扩散层的价格占比会减少一些，但即使到批量化的生产，它也占有一定的比例。当规模达到50万辆的时候，它还有6%的价格比例，所以气体扩散层对燃料电池的价格
有非常大的影响。气体扩散层影响到燃料电池里面的传质、传热以及欧姆电阻，因为它对水管理有非常大的作用，它间接的影响到燃料电池里面的动力学。因为它影响燃料电池的性能，所以间接影响到燃料电池的成本。尤其是我们现在燃料电池技术发展是朝着高电流操作这样一个方向来发展，高电流操作的时候，传质就非常重要。 如何有一个好的气体扩散层的材料？首先要分析一下它在燃料电池里面所起的作用。首先第一个作用是传质的作用，它要把气体均匀的扩散到催化剂上进行反应，同时要把这个产物这个水能够带出来，这是它的第一个作用。同时，燃料电池里面反应产生的热量也要通过气体扩散层带出去。另外所产生的电流也要通过气体扩散层导到外面的电路，传质、传热、导电是气体扩散层在燃料电池里面主要的功能。同时它需要有一定的机械强度来支撑我们膜电极，膜电极比如是CCM机械强度不够，需要有一个支撑，这个是气体扩散层四个功能。因为有这样一个功能的要求，所以我们对气体扩散层的材料选择也是有一定的要求。首先它必须是多孔的介质，需要良好的导热性，有良好的导电性，同时有一定的机械强度。另外因为它在燃料电池里面工作，所以它需要有良好的化学稳定性。正因为有这样一个要求，所以气体扩散层的选择并不是很多。早期的时候有些人用金属网或者金属泡沫来做气体扩散层，因为它的化学稳定性不高，所以后来基本上气体扩散层材料的选择只有两
7.1 扩散定律(1) 7.1.1 菲克第一定律（Fick’s First Law） 扩散过程可以分类为稳态和非稳态。 在稳态扩散中，单位时间内通过垂直于给定方向的单位面 积的净原子数（称为通量）不随时间变化，即任一点的浓度 不随时间变化。在非稳态扩散中，通量随时间而变化。研究 扩散时首先遇到的是扩 散速率问题。 菲克(A. Fick)在1855 年提出了菲克第一定律， 将扩散通量和浓度梯度 联系起来。菲克第一定律 指出，在稳态扩散（即） 的条件下，单位时间内通过垂直于扩散方向的单位面积 的扩散物质量（通称扩散通量）与该截面处的浓度梯度 成正比。为简便起见，仅考虑单向扩散问题。设扩散沿x 轴方向进行（图7-1），菲克第一定律的表达式为 （7-1） 式中：J为扩散通量(atoms/(m2·s)或kg/(m2·s))；D为扩散 系数(m2/s)；为浓度梯度(atoms/(m3·m)或kg/(m3·m)) （图7-2为浓度梯度示意图）；“-”号表示扩散方向为浓度梯度的反方向，即扩散由高浓度向低浓度区进行。此方程又称为扩散第一方程。 当扩散在稳态条件下应用（7-1）式相当方便。 7.1.2 菲克第二定律(Fick’s Second Law) 实际上，大多数重要的扩散是非稳态的，在扩散过程中扩散物质的浓度随时间而变化，即dc/dx≠0。为了研究这种情况，根据扩散物质的质量平衡，在菲克第一定律的基础上推导出了菲克第二定律，用以分析非稳态扩散。在一维情况下，菲克第二定律的表达式为 （7-2） 式中：为扩散物质的体积浓度（atoms/m3或kg/m3）；为扩散时间（s）；为扩散距离（m）。（7-2）式给出c=f(t,x)函数关系。式（7-2）又称为扩散第二方程。由扩散过程的初始条件和边界条件可求出（7-2）式的通解。利用通解可解决包括非稳态扩散的具体扩散问题。 7.1.3 扩散方程的求解
燃料电池气体扩散层设计与选型 作为膜电极的重要组成部件，气体扩散层的设计与选型需根据电堆水管理特性、极板尺寸、单体目标厚度等因素因地制宜。 气体扩散层(GDL)是一类疏水多孔介质材料，位置介于流场板和催化层，担当水气输运、热量传递、电子传导的载体，并在装配和运行过程中提供结构支撑。GDL通常由大孔基底层(Macroporous substrate:MPS)和微孔层(Microporous layer:MPL)组成。其中，基底层通常由碳纤维各向异性堆叠组成，直接与流场板接触；微孔层由碳基粉末和憎水剂混合而成，直接与催化层接触。 气体扩散层关键特性
气体扩散层通常由多孔、非编织性和大孔结构的碳基材组成，基材经PTFE 疏水处理后，并涂覆有单层或多层微孔层(MPL)。一般，质子交换膜燃料电池用气体扩散层材料应具有反应气扩散、产物水扩散传输、导电、导热和机械支撑等关键特性。 反应气扩散气体扩散层的首要任务是传送反应气氢气和氧气，确保足够的反应物质快速和均匀扩散至催化层。因此，气体扩散层的孔径在一定范围内应足够大，且孔隙需具备足够的疏水特性以避免燃料电池的产物水阻塞孔道。产物水扩散与 传输 一方面，气体扩散层需有效将液态水自催化层移至流场板(或极板)，以避免液态水阻塞反应物扩散通道引起传质极化增加。另一方面，排水特性需进行最佳设计。排水能力过强，将导致质子膜过度干燥产生“脱水”现象，质子传导率下降。 导电 气体扩散层材料导电能力高有助于降低电子传导过程中的欧姆损失。但调整气体扩散层的其他物理特性会影响到材料的导电特性，如：增加气体扩散层的孔隙率及PTFE含量时，通常导电率将下降。一般，碳基材料的导电特性可依据碳材料的热处理温度进行改善。导热
《计算材料学》实验讲义 实验二：分子动力学模拟-水分子扩散系数 一、前言 分子动力学模拟的基本思想是将物质看成是原子和分子组成的粒子系统（many-body systems ），设置初始位能模型，通过分析粒子的受力状况，计算粒子的牛顿运动方程，得到粒子的空间运动轨迹，可以求得复杂体系的热力学参数以及结构和动力学性质。分子动力学模拟的理论是统计力学中的各态历经假说(Ergodic Hypothesis)，即保守力学系统从任意初态开始运动，只要时间足够长，它将经过相空间能量曲面上的一切微观运动状态，系统力学量的系综平均等效力学量的时间平均，因此可以通过计算系综的经典运动方程来得到力学量的性质。比如，由N 个粒子组成的系综的势能计算函数为： int U U U VDW += (1-1) VDW U 表示粒子内和粒子之间的Van der Waals 相互作用；int U 表示粒子的内部势能（键角弯曲能，键伸缩能、键扭转能等）；根据经典力学方程，系统中第i 个粒子的受力大小为： U k z j y i x U F i i i i i ??? ? ????+??+??-=-?= (1-2) 那么第i 个粒子的加速度可以通过牛顿第二定律得到： ()()i i i m t F t a = (1-3) 由于体系有初始位能，每个粒子有初始位置和速度，那么加速度对时间进行积分，速度对时间积分就可以获得各个任意时刻粒子的速度和位置： i i i a v dt d r dt d ==22 (1-4) t a v v i i i +=0 (1-5) 2002 1t a t v r r i i i i ++= (1-6) i r 和v 分别是系统中粒子t 时刻的位置和速度，0i r 和0i v 分别是系统中粒子初始时刻的位置和速度。依据各态历经假说，可获得任意物理量Q 的系综平均，因此得到体系的相关性质：
对流扩散方程有限差分方法 求解对流扩散方程的差分格式有很多种，在本节中将介绍以下3种有限差分格式：中心差分格式、Samarskii 格式、Crank-Nicolson 型隐式差分格式。 3.1 中心差分格式 时间导数用向前差商、空间导数用中心差商来逼近，那么就得到了（1）式的中心差分格式]6[ 2 1 11 1122h u u u v h u u a u u n j n j n j n j n j n j n j -+-+++-=-+-τ （3） 若令 h a τ λ=，2h v τ μ=，则（3）式可改写为 )2()(2 111111 n j n j n j n j n j n j n j u u u u u u u -+-+++-+--=μλ （4） 从上式我们看到，在新的时间层1+n 上只包含了一个未知量1 +n j u ，它可以由时间层n 上的值n j u 1-，n j u ，n j u 1+直接计算出来。因此，中心差分格式是求解对 流扩散方程的显示格式。 假定),(t x u 是定解问题的充分光滑的解，将1 +n j u ，n j u 1+，n j u 1-分别在),(n j t x 处 进行Taylor 展开： )(),(),(211ττO t u t x u t x u u n j n j n j n j +??? ?????+==++ )(2),(),(3 22211 h O x u h x u h t x u t x u u n j n j n j n j n j +????????+????????+==++ )(2),(),(3 22211 h O x u h x u h t x u t x u u n j n j n j n j n j +????????+????????-==-- 代入(4)式，有 2 111 1122),(h u u u v h u u a u u t x T n j n j n j n j n j n j n j n j -+-+++---+-= τ )()()(2222 h O v x u v h O a x u a O t u n j n j n j ?-????????-?+????????++????????=τ )()()(222h O v a O x u v x u a t u n j n j n j ?-++????????-??? ?????+????????=τ
第十一讲扩散定律 1.扩散第一定律 考点再现：08、09年出现了证明扩散第一定律的题目，而10年出现了用误差函数解决扩散第二定律的问题，所以按照往年的经验，扩散第一定律和扩散第二定律必考其一，所以这部分比较重要，分值会在8至10分，题型还是简答。 考试要求：理解，记忆，并且要求会推导出扩散第一定律。 知识点 在气体或者液体中，物质的传输方式为（扩散）和（对流）。★★★★ 在固体中，物质的传输方式为（扩散）。★★ 菲克第一定律，条件——稳态扩散，即材料内部各处的浓度不随时间而变（dc/dt=0）★★★★★ 单位时间内通过垂直于扩散方向单位截面的物质流量（称为扩散通量J）与该出的浓度梯度成正比。 J为扩散通量；D为扩散系数；dc/dx为浓度梯度。 由扩散第一定律可以得到一下几点结论：★★★ （1）只要有浓度梯度，就会有扩散。 （2）扩散通量的大小与浓度梯度成正比 （3）扩散方向与浓度梯度正方向相反，即扩散的宏观流动总是从溶质浓度高的向浓度低的方向进行。 2.扩散第二定律
考点再现：10年已经出现了扩散第二定律的内容，相对来说11年考的可能性就要小一些，但是不能完全忽视，毕竟08，09还都考了扩散第一定律了，考试方式很固定，即误差函数法解扩散第二定律。 考试要求会用误差函数法解题，会计算，知道每个字母所代表的意义，对于一些题目，能够从中抽象出问题。 知识点 扩散第二定律的表达★★★ 条件，非稳态扩散，即材料内部溶质浓度随时间改变。dc/dt≠0 因为这个公式相对比较复杂，所以对于这个公式的推导并不作为考试的要求，这一部分我们只需要把公式记住就可以了。 利用误差函数分布作为扩散第二定律的解★★★★★ 现对书中的例题进行讲解
1扩散动力学方程一一菲克定律 1.1菲克第一定律 1.1.1宏观表达式 1858年，菲克（Fick ）参照了傅里叶（Fourier ）于1822年建立 的导热方程，建立定量公式 在t 时间内，沿x 方向通过x 处截面所迁移的物质的量 m 与x 处 分布 的浓度梯度成正比： C m A t x 即如 D （_C ） Adt x 根据上式引入扩散通量概念，则 1 . 1 1 1 1 有: (7-1) 图7-1扩散过程中溶质原子的 ( C-C) 繞扩ft 石 原始状畚 盘蚌#态
式(7-1)即菲克第一定律。 式中J 称为扩散通量，常用单位是mol / ( cm 2 s); 散原子面密度分别为n 1和n 2,若n 1=巳，则无净扩散流。 假定原子在平衡位置的振动周期为 T 则一个原子单位时间内离 开相对平衡位置跃迁次数的平均值，即跃迁频率 为 丄 (7-2) 由于每个坐标轴有正、负两个方 向，所以向给定坐标轴正向跃迁的几率 是 1 。 6 设由平面I 向平面2的跳动原子通 量为J 12，由平面2向平面1的跳动原 模型 -C 浓度梯度； x D 扩散系数，它表示单位浓度梯度下的 通量，单位为 cm 2/s 或 m 2 / s ； 负号表示扩散方向与浓度梯度方向相 反见图7-2。 1.1.2微观表达式 微观模型： 设任选的参考平面1、平面2上扩 图7-2溶质原子流动 的方向与浓度降低的方 向相一致 图7-3 一维扩散的微观
子通量为J 21 (见图7-3)，贝卩由式(7-5)、式(7-6)得 式(7-8)反映了扩散系数与晶体结构微观参量之间的关系，是扩散 系数的微观表达式。 三维情况下，对于各向同性材料(D 相同)，则 C C C /-7 C\ J J x J y J z D(i j k ) D C (7-9) XXX 式中： i j k 为梯度算符。 x x x 对于各向异性材料，扩散系数 D 为二阶张量，这时， J 12 1 6n i 1 6n 2 (7-3) 注意到正、反两个方向，则通过平面 J i J 12 J 21 而浓度可表示为 1沿x 方向的扩散通量为 (7-5) (7-6) 式(7-6)中的1表示取代单位面积计算, 表示沿扩散方向的跳动距离 J 1 C 1 C 2 1 6(6 C 1) 2 dC dx 式(7-7)即菲克第一定律的微观表达式，其中 dC D (7-7) dx
（1 菲克第一定律 扩散流密度与扩散组元浓度梯度间关系 称为菲克第一定律。扩散流密度与在扩散介质中的浓度梯度成正比, 比例常数称为扩散系数。 菲克第二定律 稳态扩散特征是0 dc dt =。在物质的浓度随时间变化的体系中，即0 dc dt ≠，体系中发生的是非稳态扩散。在一维体系中，单位体积单位时间浓度随的变化等于在该方向 上通量，这既是菲克第二定律，其数学表达式为,A A x c t J x ????= )A A A ( x c D t c x ??????= 若D A 为常数, 即可以忽略D A 随浓度及距离的变化, 在x-y-z 三维空间中, 则菲克第二定律的表示式为 （2）掌握 D 为常数时费克第二定律的几个特解 扩散偶问题 如图4-1-2 初始条件 ｔ=0，x >0，c =0 ; 边界条件 t >0, x =0, c = c 02 ; x =∞, c = 0 解方程????c t D c x =2 2 ，得 ) d π2 1(2202 0ξξ?--=Dt x e c c 不同扩散时间后,扩散偶中扩散组元的浓度分布 ξ ξd π 2202 ?-Dt x e 为积分函数 。 (式中Dt x 2= ξ)称为误差函数, 记作Dt x 2erf 。 于是 )2e r f 1(2 ),(0Dt x c t x c -=
注：误差函数有如下主要性质 erf(x )= λλ d π 22 -?e x erf(-x )= - erf(x ) erf(0)=0, erf(∝) =1 1-erf(x )= erfc(x ) erfc(∝)=0, erfc(0)=1 式中 erfc(x )称为余误差函数。 若初始条件变为t =0, x ＞0，c =c 1则解为 )2erf 1(2),(101Dt x c c c t x c --+= 几何面源问题 数学模型1 初始条件： t =0, x =0, c =c 0；x ≠0, c =0 Vc 0 =Q 式中V ? 极薄扩散源的体积; Q ? x =0处扩散组元的总量。如图4-1-2所示。 边界条件：t >0, x →∞, c =0; x →-∞, c = 几何面源、全无限长一 维扩散 (a) 边界条件; (b) 浓度分布曲线（扩散时间 t =1, 14 , 164 , 横坐标距离x 为任意长度位置） 由初始及边界条件得到的菲克第二定律的解为 Dt x e Dt Q c 42 2- =π 数学模型2 初始条件： t = 0，x = 0，c =c 0，Ｑ=Vc 0; x >0，c = 0 边界条件： t > 0，x =∞，c = 0 所得的菲克第二定律的解为 c Q D t e x D t = -π2 4 数学模型3 t = 0，x ≥0，c = c b 0 < t ≤ t e ，x =0，c =c s ; x =∞，c =c b 菲克第二定律的解为 c c c c x D t --=-b s b er f 12( ) 或
7、2、2扩散系数 费克定律中的扩散系数Ｄ代表单位浓度梯度下的扩散通量,它表达某个组分在介质中扩散的快慢,就是物质的一种传递性质。 一、气体中的扩散系数 气体中的扩散系数与系统、温度与压力有关,其量级为5 2 10/m s -。通常对于二元气体Ａ、B 的相互扩散,Ａ在B 中的扩散系数与B 在A 中的扩散系数相等,因此可略去下标而用同一符号Ｄ表示,即AB BA D D D ==。 表７－１给出了某些二元气体在常压下(5 1.01310Pa ?)的扩散系数。 对于二元气体扩散系数的估算,通常用较简单的由富勒(Fuller)等提出的公式 : 1/31/32 [()()]A B D P v v = +∑∑ (７－１９) 式中,D －Ａ、B 二元气体的扩散系数,2 /m s ; P －气体的总压,Pa ; T －气体的温度,Ｋ; A M 、 B M －组分Ａ、B 的摩尔质量,/kg kmol ; A v ∑、B v ∑－组分Ａ、B 分子扩散体积,3 /cm mol 。 一般有机化合物可按分子式由表７－２查相应的原子扩散体积加与得到,某些简单物质则在表７－２种直接列出。 5 1.01310Pa ?
式７－１９的相对误差一般小于１０％。 二、液体中的扩散系数 由于液体中的分子要比气体中的分子密集得多,因此也体的扩散系数要比气体的小得多, 其量级为92 10/m s -。表７－３给出了某些溶质在液体溶剂中的扩散系数。 对于很稀的非电解质溶液(溶质Ａ＋溶剂Ｂ),其扩散系数常用Wilke-Chang 公式估算: 15 0.6()7.410 T B AB A M T D V -φ=?μ 2/m s (７－２１) 式中,AB D －溶质Ａ在溶剂Ｂ中的扩散系数(也称无限稀释扩散系数),2 /m s ;
DOI：10.3969/j.issn.1000-4874.2010.02.001
水动力学研究与进展A辑2009年第2期 128 1 引言 对流-扩散方程是描述粘性流体运动的非线性Burgers方程的线性化模型,它可以刻画许多自然现象,如：水体和大气中污染物的输移、扩散和降解，海水盐度和温度的扩散，流体流动与传热和电化学反应等。研究对流-扩散模型具有重要的理论价值和实际意义,它已经广泛应用于环境科学、能源开发、流体力学和电子科学等领域。总的来说，目前关于对流-扩散方程的研究大致可以分为两个方面。一方面是在给定初边值条件下，通过不同的数值计算方法求解对流-扩散方程，以模拟研究对象(例如：温度、盐度和污染物等)在时间和空间上的发展演化，这类问题可以统称为正问题。迄今为止已经有很多成熟方法求解对流-扩散方程,如有限差分方法(FDM)[1,2,3]、有限体积方法(FVM)[4,5,6]和有限元方法(FEM)[7,8,9]等。 另外一方面是关于对流-扩散方程反问题的研究，即通过所研究对象的观测资料来估计和识别方程中的参数、源项、边界和初始条件等。从某种意义上讲，反问题的求解是对流-扩散模型研究中一个更重要的问题,因为它的正确与否直接影响模型的可靠性。由于偏微分方程反问题固有的非线性和不适定性[10], 对流-扩散方程反问题的求解会存在巨大困难，通常的方法常常导致求解失败。近年来, 国内外学者关于对流-扩散反问题开展了广泛研究。Andreas Kirsch对一维扩散方程逆过程反问题进行了稳定性分析，并给出了误差估计公式[11]。Yildiz[12]、刘继军等[13-16]对相关问题采用Tikhonov 正则化方法进行了深入研究。闵涛等[17]以函数逼近和Tikhonov正则化为基础，利用算子识别摄动法和线性化技术，建立了河流水质纵向弥散系数反问题的迭代算法，并进行了数值试验。闵涛等[18]利用有限元法求解了二维稳态对流-扩散方程,并利用迭代法对二维稳态对流-扩散方程参数反演进行了研究。闵涛等[19]利用遗传算法就对流-扩散方程的源项识别反问题进行了研究。潘军峰等[20]对一维对流-扩散方程的反问题利用Tikhonov正则化方法进行了研究。吴自库等[21]结合利用伴随同化方法和处理数学物理反问题的技巧就对流-扩散方程逆过程的反问题进行了数值研究。综上所述，由于对流-扩散方程反问题的不适定性，所以它的求解一般要采用特殊方法，如Tikhonov正则化方法、变分伴随方法和遗传算法等等。本文在贝叶斯理论的基础上，提出采用马尔科夫链蒙特卡罗(Markov Chain Monte Carlo,简称MCMC)方法[22,23]来识别对流-扩散方程中多个点源中的未知参数。 结合利用贝叶斯方法和MCMC算法求解反问题，具有以下优点：1) 能方便地将各种先验信息和误差信息高效地融合到问题求解过程中，减小问题的不确定性；2) 和确定性算法不同，反问题的不适定性不再是MCMC算法要考虑的问题，且计算获得的是全局最可能解，而通常的最优化算法可能陷入目标函数局部极小值；3) 能对定义在高维空间且无明确数学表达式的概率分布密度函数进行数值计算，而确定性方法无法解决此类问题；4) MCMC算法通过构造Markov链来进行随机模拟，是一种动态Monte Carlo方法，计算速度高于一般的Monte Carlo 方法和模拟退火算法，而且计算复杂度不依赖于计算空间的维数。 2 反问题模型 不失一般性，用对流-扩散方程来模拟污染物在河道中的扩散，考虑对流-扩散方程的初边值问题[19,21]，公式如下： 2 2 1 (), (,)(0,)(0,) (0,)0,(,)0,(0,) (,0)0,(0,) q i i i C C C u E kC s x x t x x x t L T C t C L t t T C x x L δ = ???? +=?+? ???? ?? ∈× ? ?==∈ ? =∈ ?? ∑ (1) 其中C为污染物的浓度，u为流速，E为扩散系数，k为污染物的降解率，L表示河道长度。δ是狄拉 克函数， i x和 i s，(1,2,) i q = 分别表示多个点污染源的位置和排放强度。假定(,) C x t在t T =时的分布已知，那么源项识别反问题就是根据这些已知 distribution, the Adaptive Metropolis algorithm was used to construct the Markov Chains of unknown parameters. And the converged samples were used to estimate the unknown parameters of source term. The results of numerical experiments show that the method has many virtues, such as high accuracy, quick convergent speed and easy to program and implement with computer． Key words: convection-diffusion equation; source term; inverse problem; Markov Chain Monte Carlo method