有限差分求解扩散的数值解

C2
Jeff
C2
无孔结构
J0 =-D ∂C ∂z
多孔介质
ε ∂C Jeff =- D τ ∂z
其中ε为孔隙率
曲折因子
τ= J0 ε Jeff
3.扩散模拟表征多孔介质物质输运性能（相转换流延制备） SRCT构建孔结构数组
灰度图
二值图
黑色：孔 白色：实体 孔隙率：0.35
三维结构重构图A
3.扩散模拟表征多孔介质物质输运性能
1.有限差分的基本原理（ 一维扩散例子） 计算结果
稳态下
∂C ∂2 C =D 2 = 0 ∂t ∂x
2.边界条件处理（ 二维柱坐标扩散例子）
扩散进口
扩散出口
绝缘 环
圆柱上底面和侧面为高浓度面（红色） 圆柱下底面同心小圆为低浓度面（白色） 圆柱下底面圆环为绝缘面 （蓝色）
2.边界条件处理（ 二维柱坐标扩散例子） 0 Z
r
扩散微分方程： ∂C ∂2 C 1 ∂C ∂2 C =D( 2 + r + 2) ∂t ∂r ∂r ∂Z
圆心对称
柱坐标下取体积微元由质量衡算推导扩散偏微分方程
r+ △r
△Z
r
体积微元质量变化量
∆m=-2π r+∆r ∆Z∆t*J -π r+∆r 2 -r2 ∆t*J
移项化简得
r+∆r+2πr∆Z∆t*J r Z
2 2 + π r+ ∆r -r ∆t*J Z+∆Z
∆m 1 -J = * 2πr∆Z∆r∆z r
r+∆r *
即
r+∆r -J r *r J Z+∆Z -J Z ∆r ∆Z ∂C ，将 J=-D 代入可得 ∂x
∂C ∂2 C 1 ∂C ∂2 C =D( 2 + + ) ∂t ∂r r ∂r ∂Z2
2.边界条件处理（ 二维柱坐标扩散例子） 0
方向六个点Ai ±1,j±1,k±1的结构， 如果周围某个点为实体，则令其 浓度等于计算点Ai,j,k的浓度Ci,j,k 例如： Ai +1,j,k ， Ai -1,j,k ， Ai,j+1,k， Ai,j-1,k四个点为实体，则令其 Ci +1,j,k = Ci -1,j,k = Ci,j+1,k= Ci,j-1,k = Ci,j,k 则迭代时简化为
L
X轴
设法在x—t平面网格点 上逼近浓度C(x,t)值
1.有限差分的基本原理（ 将微分方程转化为代数方程）
Cik,
j+1 h -Cik,jh
h Ci,j+1 = D*h k
2
=D*(
C i+1
k,jh -2Cik,jh +C i-1 k,jh ) 2
k
Ci+1,j +Ci-1,j + 1-2*
D*h k
1.有限差分的基本原理（ 一维扩散例子）
高浓度端 10-5 低浓度端 0
设定参数
扩散系数D： 10-7cm2/s 透氧距离L： 0.01cm 高浓度C1： 10-5 mol/ml 低浓度C2： 0 距离步长k： 0.0005cm 距离等分数： 20 时间步长h： 1s 初始浓度： 0
迭代代数式，Matlab编程计算 Ci,j+1 =s Ci+1,j +Ci-1,j + 1-2*s *Ci,j
Ckm,j Ckm,j+1
Cki,j-1
Cki,j Cki,j+1
Cki+1,j
虚拟点：Ckm+1,j=Ckm,j
2.边界四
圆柱轴中心线上格点
Cki-1,1
Cki,1 Cki,2
Cki+1,1
无法用
2.边界四
圆柱轴中心线上格点微元体积质量衡算
Cki-1,1 Cki,1 Cki+1,1 Cki,2
孔 A
0<Cs<1
B
Cs=1
固相
Cs=初始值
C
D
Cs=0
下表面浓度：0
通过判断稳态时各网格点浓度值区分四类孔 初始浓度设为非0非1的数，以区分孤立孔 死孔包含B、C和D
3.扩散模拟表征多孔介质物质输运性能
孔弯曲度（曲折因子）
曲折因子τ=1
曲折因子τ>1
气体扩散测定曲折因子
C1 C1>C2 C1
J0
k k k k Ck+1 i-1,1 +Ci+1,1 +2Ci,2 -4Ci,1 +Ci,1
2.Matlab计算结果
50秒浓度及其等高线
2.Matlab计算结果
200秒浓度及其等高线
2.Matlab计算结果
500秒浓度及其等高线
2.Matlab计算结果
1200秒浓度及其等高线
3.扩散模拟表征多孔介质物质输运性能
2
*Ci,j
令
D*h k
2
=s
，得代数迭代格式
Ci,j+1 =s Ci+1,j +Ci-1,j + 1-2*s *Ci,j
Ci,j+1
由上式可知，某一时刻的浓度值，可以由前一 个时刻的浓度值计算得出，即可从初始零时刻 逐步迭代计算任一时刻任一离散点的浓度值
Ci-1,j
Ci,j
Ci+1,j
注意：要使迭代式收敛稳定，须使0<s ≤0.5
浓度不变的网格点构成的孔为死孔
3.扩散模拟表征多孔介质物质输运性能
内部孔结构
3.扩散模拟表征多孔介质物质输运性能（相转换流延制备）
纵向剖面结构图
纵向剖面浓度分布
3.扩散模拟表征多孔介质物质输运性能（孔随机分布）
孔随机分布三维结构 黑色：孔 白色：实体 孔隙率：0.35
3.扩散模拟表征多孔介质物质输运性能（孔随机分布）
体积微元质量变化量
∆m=Dπ h ∆t*
移项化简得
2
Ci-1,1 -2Ci,1 +Ci+1,1 h
+2πh ∆t
2
Ci,2 -Ci,1 h
∂C D = = 2 Ci-1,1 -4Ci,1 +Ci+1,1 +2Ci,2 3 πh ∆t ∂t h D∆t h2
∆m
得到轴中心线上网格点迭代格式：
Ck+1 i,1 =
模拟计算原理 原理：三维扩散方程
∂C ∂2 C ∂2 C ∂2 C =D + + ∂t ∂x2 ∂y2 ∂z2
Ci-1,j,k Ci,j+1,k
y

x z
Ci,j,k-1 Ci,j-1,k
Ci+1,j,k
计算方法：有限差分
Ci,j,k+1
n n n n n n n n n Cn+1 Cn i,j,k -Ci,j,k i+1,j,k -2Ci,j,k +Ci-1,j,k Ci,j+1,k-2Ci,j,k +Ci,j-1,k Ci,j,k+1 -2Ci,j,k +Ci,j,k-1 =D + + 2 2 ∆t ∆x ∆y ∆z2
D∆t ∆x2
n n Cn i,j,k+1 +Ci,j,k-1 -2Ci,j,k
Ci-1,j,k
Ci+1,j,k Ci,j+1,k
Ci,j,k+1
n Cn+1 i,j,k =Ci,j,k +
为一维扩散方程离散代数表达式
3.扩散模拟表征多孔介质物质输运性能（相转换流延制备）
根据以上迭代式逐个求每个时间点网格点的浓度，直到进口 流量和出口流量相等，即到达稳态。 稳态时浓度分布图
Cki-1,j Cki,j-1 Cki,j Cki+1,j Cki,j+1
Cki-1,j
Cki,j-1
Cki,j Cki,j+1 Cki+1,j
变量离散，微分方程化为代数方程：
k+1 D∆t Ci,j = 2 k k k k Ck +C +C +C -4C i,j+1 i,j-1 i+1,j i-1,j i,j +
稳态时浓度分布图
3.扩散模拟表征多孔介质物质输运性能（孔随机分布）
纵向剖面结构图
纵向剖面浓度分布
3.扩散模拟表征多孔介质物质输运性能
两种孔结构输运性能对比
样品 相转换制备 孔随机分布 孔隙率 0.35 0.35 总曲折因子 2.15 70.67 连通孔/总孔 95.12% 65.69% 封闭孔/总孔 1.68% 26.73%
h
1 k Ck -C j-1 i+1,j i,j
+Ck i,j
2.边界一
圆柱上表面和侧面边界
Ck1,1——Ck1,n=C0 Ck1,n——Ckm,n=C0
2.边界二
扩散出口边界条件
扩散出口边界
Ckm,1——Ckm,n1=0
2.边界三
下表面绝缘边界条件：
Ckm-1,j
Cki-1,j
Ckm,j-1
实验室制备多孔陶瓷膜 A B
C A：样品实物图 B：SRCT三维重构图 C：孔结构图
3.扩散模拟表征多孔介质物质输运性能
孔连通性
孔 A
B
固相
C
D
A：与样品上下表面皆连通的孔 C：封闭在样品内的孤立孔 B：与样品上表面连通的孔 D：与样品下表面连通的孔
3.扩散模拟表征多孔介质物质输运性能
上表面浓度：1
C t+∆t -C t (C x+∆x -2C x +C =D* ∆t ∆x2
x-∆x )
1.有限差分的基本原理 （变量离散）
t轴
在x—t平面将距离和时间离 散化：时间步长△t取h，距 离步长△x取k，步数分别为n 和m。
Ci,j+1
Ci-1,j Ci,j Ci+1,j Ci,j-1
tj=j*h
Xi=i*k m=L/k，n=T/h,平面内有m×n个网格点
令x,y,z方向步长值相等， ∆x=∆y=∆z 简化上式得
n Cn+1 i,j,k =Ci,j,k +
D∆t ∆x2
n n n n n n Cn i+1,j,k +Ci-1,j,k +Ci,j+1,k +Ci,j-1,k +Ci,j,k+1 +Ci,j,k-1 -6Ci,j,k
含孔结构计算处理技巧
n Cn+1 i,j,k =Ci,j,k +
D∆t ∆x2
n n n n n n Cn i+1,j,k +Ci-1,j,k +Ci,j+1,k +Ci,j-1,k +Ci,j,k+1 +Ci,j,k-1 -6Ci,j,k
x y
z
Ci,j,k-1
计算 Cn+1 i,j,k 前，先判断该点三维
Ci,j-1,k
谢谢大家！
有限差分法求解扩散问题的 数值解
王智文 2013.5.15
提 纲
1.有限差分原理（一维坐标系）
2.边界条件处理（二维柱坐标系） 3.实验室应用举例（三维直角坐标系）
1.有限差分的基本原理（一维坐标系）
菲克第二定律的一维方程：
∂C ∂2 C =D 2 ∂t ∂x
∆x → 0, ∆t → 0 时 微分方程转化为代数方程：