剖面二维非恒定悬移质泥沙扩散方程的数值方法<
1 引言
数学模拟方法正在成为研究河流泥沙问题的重要手段。目前，一维数学模型发展较成熟，已广泛应用于模拟长河段的长期变形，但它只能给出河段平均冲淤深度的沿程变化，如需了解短河段的河床变形细节，则要采用二维以至三维数学模型。不论是一维数学模型还是平面二维维数学模型，都不能反映含沙量沿垂线的分布状况，并忽略了含沙量沿垂线分布对垂线平均含沙量变化过程的影响。要解决这类问题，必须建立剖面二维数学模型。这种模型主要通过解剖面二维泥沙扩散方程来研究悬移质泥沙沿水深的分布及含沙量的变化过程，对水电站进口和其它引水工程的引水口高程的确定都能提供较好的数值模拟。
泥沙扩散方程实际上是一个变系数的二阶线性偏微分方程，这样的方程在各种复杂边界条件下求解是极为困难的。求扩散方程的解析解在数学上存在着难以克服的困难，往往只能通过对方程的简化，才能得到一些简单边界条件下的解析解，在这方面， A.A.Kalinske 、野满隆治、 W.E.Dobbins 、俞维强、张启舜、韦直林 [ ２ ] 等都做了有益的尝试；求扩散方程的数值解曾经因为缺乏高效率的计算工具而难以实现，直到 60 年代后，随着计算机的广泛应用，在各种复杂边界条件下求扩散方程的数值解不但成为可能，而且得到迅速的发展，在这方面，曹志先、崔侠 [ ４ ] 等做了大量工作，取得了很多成果。
数值方法相对于解析方法在求解偏微分方程上有着明显的优势，即简单灵活、计算方便快捷，但要寻找一种精度高、稳定性好、计算方便的差分格式也并非易事。本文拟在前人研究的基础上着重讨论剖面二维泥沙扩散方程的数值解问题，希望能提供一种精度高、稳定性好、计算方便的数值解。
2 基本方程
剖面二维泥沙扩散方程的形式为
(1) 式中 x,y 为水流方向和铅直方向的维轴； u,v 分别为沿水流方向和铅直方向
的时均流速；ε
sx , ε
sy
分别为水流方向和铅直方向的泥沙扩散系数；ω ,S
分别为泥沙静水沉速和含沙量。
对于式 (1) 的求解，研究者一般会对它进行不同程度的简化，为此引入以下假定中的一种或几种：
A. 非恒定流可以概化为梯级式恒定流，即
;
B. 在一个时段内，认为泥沙运动可以概化为处于恒定状态，即
;
C. 在二维流动中，纵向扩散系数与方程其他项相比，可以忽略不计，即认为方程右端第一项可以忽略；
D. 认为悬移质泥沙粒径均一，即ω =const;
E. 认为水流为二维均匀流，即 v=0 。
为简单起见，我们讨论的范围限于水流条件为二维非恒定均匀流，悬移质泥沙粒径均匀，为此引入假定 C 、 D 、 E 。这时，泥沙扩散方程为
(2) 目前，对ε
的变化规律研究得不很充分，一般假定
s
ε s = βε m
(3) 其中ε
为动量传递系数，β为修正值。由勃兰特尔掺长理论可得 m
ε m = κ u*y/(1-y/h)
(4)
式中κ为卡门常数， u
*
为摩阻流速。
对于 u ，我们取卡曼 - 勃兰特尔对数流速分布公式
(u
max -u)/u
*
=(1/ κ )ln(h/y)
(5)
令 W= ω + β ( κ u
*
/h)(1-则式 (2) 可变形为
(6) 3 差分方程
3.1 网格的剖分
为建立差分方程，首先必须剖分网格。我们取时间步长Δ t= τ， X 方
向的空间步长Δ x=h
1， Y 方向的空间步长为Δ y=h
2
, 这样形成如下网格
D h ={(x
j
,y
l
,t
n
) ｜ x
j
=jh
1
,y
l
=lh
2
,t
n
=n τ｜ } 其中j=0,...,N; l=0,...,M;n ≥ 0;
3.2 构造差分格式
通过对流方程和扩散方程的差分格式的构造，我们可以得到对流扩散方程的差分格式。由于隐式格式稳定性好，考虑 Crank-Nicholson 型隐式格式。为此，引入差分算子记号
δ 2 y S n j,l=S n j,l+1-2S n j,l+S n j,l-1
L x S n
j,l
=S n
j+1,l
-S n
j-1,l
L y S n
j,l
=S n
j,l+1
-S n
j,l-1
为了看得更清楚，暂且取 h
1=h
2
=h. 对式 (6) 离散，则 C-N 格式为
S n+1
j,l -S n
j,l
/ τ +u/4hL
x
(S n
j,l
+S n
+1j,l
)= ε
s
/2h2δ 2
y (S n+1
j,l
+S n
j,l
)+W/4hL
y
(S n+1
j,l
+S n
j,l
)
(7)
C-N 格式的精度是二阶的，绝对稳定。但对于二维问题，由 (7) 导出的方程组，其系数矩阵不是三对角矩阵，不能用追赶法求解。因此，考虑构造交替方
向的隐式格式 ( 命名为 Z-C 格式 )
(S n+1/2
j,l -S n
j,l
)/( τ /2)+u/2hL
x
S n
j,l
= ε
s
/h2δ 2
y
S n+1/2 j,l +W/2hL
y
S n+1/2 j,l
(8)
(S n+1
j,l -S n+1/2
j,l
)/ τ /2+u/2hL
x
S n+1
j,l
= ε
s
/h2δ 2
y
S n+1/2
j,l
+W/2hL
y
S n+1/2
j,l
(9)
可以看出，计算 S n+1
j,l
是由两步组成的，每一步仅是一个方向的隐式，故用两次追赶法即可。
3.3 精度分析
现在，我们考虑 Z-C 格式的精度。先设法消去过渡值 S n+1/2
j,l
，为此，将(8) 和 (9) 两式相加，可得
S n+1
j,l -S n
j,l
+ τ u/2hL
x
(S n+1
j,l
+S n
j,l
)= τε
s
/2h2δ 2
y
S n+1/2
j,l
+ τ W/4hL
x
S n+1/2
j,l
(10) 将 (8) 和 (9) 两式相减，可得
S n+1/2
j,l =(S n+1
j,l
+S n
j,l
)/2+ τε
s
/4h2δ 2
y
(S n
j,l
-S n+1
j,l
)+ τ W/8hL
y
(S n
j,l
-S n+1
j,l
)
(11)
把式 (11) 代入 (10) ，变形整理，可得
(S n+1
j,l -S n
j,l
)(1- τ u ε
s
/8h3L
x
δ 2
y
- τ 2 uW/16h2L x L
y
)=(S n
j,l
+S n+1
j,l
)( τ u ε s /2h2δ 2 y + τ W/4hL y- τ u/4hL x)
(12)
设 S(x,y,t) 是 (12) 的精确解，并假定 S(x,y,t) 关于 t 三次连续可微，关于 x,y 四次连续可微，那么利用 Taylor 级数展开可得
S(x
j ,y
l
,t
n+1
)-S(x
j
,y
l
,t
n
)/ τ (1- τ 2 ε s /8h3L xδ 2 y - τ 2 uW/16h2L x L y)-
S(x
j ,y
l
,tn+1)-S(x
j
,y
l
,t
n
)/ τ ( τε
s
/2h2) δ 2
y
+ τ W/4hL
y
- τ
u/4hL
x
)=O( τ 2 +h2)
(13) 由此可见， Z-C 格式具有二阶精度。
3.4 稳定性分析
现在，我们来讨论 Z-C 格式的稳定性。为此，把式 (12) 变形整理得
(1- τε
s
/2h 2 δ 2 y - τ W/4hL y)(1+ τ u/4hL x)S n+1j,l=(1+ τε s /2h 2 δ y 2
+ τ W/4hL
y )(1- τ u/4hL
x
)S n
j,l
(14)
由式 (14) 可得出过渡因子为
G( τ ,k)=[(1-2 τε
s /h 2 sin2k2h/2+i τ W/2hsink
2
h)(1-i τ
u/2hsink
1h)]/(1+2 τε
s
/h2sin2k2h/2-i τ W/2hsink
2
h)(1+i τ u/2hsink
1
h)
(15)
令 a=2 τε
s /h2sin2k2h/2,b= τ W/2hsink2h,c= τ u/2hsink
1
h, 则
|G|2=|1-a2-b2+2bi/(1+a)2+b2|2· |1-c2-2ci/1+c2|2=1 · (1-a2-b2)2+4b2/[(1+a)2+b2]2=1-4a2+4a+4b2+4a3/[(1+a)2+b2]2
(16) 显然，对于任意的τ ,h,|G( τ ,k)|2所以 Z-C 格式是绝对稳定的。
我们考虑初边值问题。
(1) 初始条件
用 Rouse 公式给出含沙量沿垂线分布
S(x,y,0)=S
a
· 5(a/h-a)z(h-y/y)z
(17)
式中 z= ω /ku
*为悬浮指标， S
a
为近底含沙量， h 为水深，一般取
a=0.01 ～ 0.05h 。
(2) 水面条件
(y=h)
(18) (3) 底部边界条件
(y=a)
(19)
式中 S
a * 为近底挟沙力，即输沙平衡时的近底售沙量 S
a
.
(4) 近底含沙量计算近底含沙量在求解泥沙扩散方程时具有边界条件性质，它
选取的正确与否，意味着所给边界条件是否正确。实际工程中一般缺乏实测资料，近底含沙量不易测定。这里，我们利用水流挟沙力和含沙量沿垂线分布公
式来反求近底含沙量 [ ４ ] 。
已知断面平均挟沙力为
S
*
=k(u3/gh ω )m
(20)
假定
S
a
*= α S*
(21)
输沙平衡时，含沙量沿垂线分布用 Rouse 公式 (17) 表示，用 (17) 表达挟沙力的垂线分布，然后沿垂线积分得断面平均挟沙力为
(22)
将 (22) 与 (21) 比较，可得
(23) 4.2 计算步骤
为方便计算，将式 (8) 和 (9) 式变形整理，并对 X ， Y 方向取不同的
空间步长
( τε
s /2h2
2
- τ W/4h
2
)Sn+1/2
j,l+1
-( τε
s
/h2
2
+1)S n+1/2
j,l
+( τε
s
/2h2
2
- τ
W/4h
2
)Sn+1/2
j,l+1
= τ u/4h
1
L
x
S n
j,l
-S n
j,l
(24)
- τ u/4h
1S n
j-1,l
+S n
j,l
+ τ u/4h
1
S n
j+1,l
= τε
s
/2h2
2
ζ 2
y
S n+1/2
j,l
+ τ
W/4h
2
L
y
S n+1/2
j,l
+S n+1/2
j,l
(25)
在一个时间层 ( 第 n 层 ) 内，计算分两步进行：
第一步，对式 (24) 用追赶法求第 n+1/2 层的过渡值。
令 C
1=- τ u/4h
1
,C
2
=1,C
3
= τ u/4h
1
,E
1
= τε
s
/2h2
2
,E
2
= τ W/4h
2
D
1= τε
s
/2h2
2
- τ W/4h
2
,D
2
=- τε
s
/h2
2
-1,D
3
= τε
s
/2h2
2
+ τ W/4h
2
l=1 时， D
1S n+1/2
j,0
+D
2
S n+1/2
j,1
+D
3
S n+1/2
j,2
=C
3
L
x
S n
j,1
-S n
j,1
l=2 时， D
1S n+1/2
j,0
+D
2
S n+1/2
j,1
+D
3
S n+1/2
j,2
=C
3
L
x
S n
j,1
-S n
j,1
……………
l=M 时 ,D 1S n+1/2j,M-1+D 2S n+1/2j,M +D 3S n+1/2j,M+1=C 3L x S n j,M -S n j,M
令 H l =-S n j,l +C 3L x S n j,l (1 H 1=-S n j,1+C 3L x S n j,1-D 1S n+1/2j,0 H M =-S n j,M +C 3L x S n j,M -D 3S n+1/2j,M+1 其中 S n+1/2j,0 和 S n+1/2j,M+1 由边界条件给出，则用矩阵形式表示为 （ 26 ） 第二步，再对 (25) 式用追赶法求第 n+1 层的值 令 F j =E 1 δ 2 y S n+1/2 j,l +E 2L y S n+1/2j,l +S n+1/2j,l F 1=E 1 δ 2 y S n+1/2 1,l +E 2 L y S n+1/2 1,l +S n+1/2 1,l -C 1 S n+1 0,l F N =E 1 δ 2 y S n+1/2 N,l +E 2 L y S n+1/2 N,l +S n+1/2 N,l -C 3 S n+1 N+1,l 其中 Sn+10,l 和 Sn+1N+1 ， l 由边界条件给出，则同理可得矩阵方程 (27) 这样，按此步骤一层层地计算。 4.3 数值模拟合理性分析 受所掌握的实测资料的限制，目前尚无法对本文提出的算法与含沙量沿垂线分布的实测值进行对比。我们用库里·阿雷克沉沙池 [ ５ ] 的实测资料作了垂线平均值沿程变化的比较。该沉沙池的主要数据为：池深 h=1.53m; 平均流速 u=0.12m/s; 泥沙沉速ω =0.0176cm/s; 悬浮指标 Z=0.01 。计算时取卡门常数κ =0.4 ， a=0.05h 。表 1 给出了计算值和实测值，结果表明，计算值和实测值比较符合。 表 1 断面平均含沙量验证 ( 单位： kg/m3) Verifications of cross-sectional average sediment concentrations 距离 (m) 200 400 600 800 1000 1200 计算值3.012 2.573 2.071 1.639 1.209 0.849 0.539 实测值3.012 2.650 1.810 1.520 1.141 0.822 0.561 为了进一步分析含沙量垂线分布计算结果的合理性，我们对另一组较粗的泥沙 ( ω =0.616cm/s,Z=0.4) 进行了对比计算。图 1\, 图 2 分别为两组沙的计算结果。从图中可以看出，计算结果符合含沙量沿垂线分布的一般规 律，粗沙分布不均匀，细沙分布较均匀；近底浓度相对较大，水面浓度相对较小，不存在 Rouse 公式中水面含沙量为 0 的缺陷；含沙量沿程衰减的特性较为明显。图 3 为较粗一组泥沙的相对含沙量沿垂线分布的沿程变化情况。图 3 表明，尽管进口断面按 Rouse 公式给出了含沙量沿垂线的分布，但由于该断面实际处于不平衡输沙状态，这种分布并不是稳定的。在距进口 200m 处，泥沙的分布调整到一种不平衡输沙状态，随着泥沙的沿程淤积，水流输沙向平衡方向发展，垂线平均含沙量趋向于水流挟沙力，而含沙量沿垂线分布向平衡时的分布状态 (Rouse 公式 ) 发展。由于这种发展是趋向于稳定状态，因此愈接近下游，分布愈靠近 Rouse 公式。计算结果表明，本文提出的计算方法是合理可行的。 图 1 含沙量垂线分布的沿程变化 (Z=0.01) 图 2 含沙量垂线分布的沿程变化 (Z=0.4) concentrations(Z=0.01) 5 结语 本文建立了求解二维非恒定泥沙扩散方程的一种差分格式 (Z-C 格式 ) 。这种格式具有如下特点： 1. 精度较高 ( 具有二阶精度 ) 。 2. 稳定性好 ( 无条件稳定 ) 。 3. 计算较方便 ( 每一时段利用两次追赶法即可 ) 。 图 3 相对含沙量的垂线分布变化 (Z=0.4) Changes of vertical distributions of relative sediment concentrations(Z=0.4) 徐州工程学院 课程设计报告 课程名称偏微分方程数值解 课题名称对流扩散方程 的迎风格式的推导和求解专业信息与计算科学 班级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>式，求出其他层的函数值。误差转化成图的形式，并输出最大值。 针对三种不同的输入对应输出结果 ： 扩散方程稳态扩散与非稳态扩散 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.1抛物型扩散方程 抛物型偏微分方程是一类重要的偏微分方程。考虑一维热传导方程： 22(),0u u a f x t T t x ??=+<≤?? (1.1.1) 其中a 是常数，()f x 是给定的连续函数。按照初边值条件的不同给法，可将(1.1.1)的定解分为两类： 第一，初值问题（Cauchy 问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ?=0,, ∞<<∞-x (1.1.2) 第二，初边值问题（也称混合问题）：求足够光滑的函数()t x u ,，满足方程(1.1.1)和初始条件： ()()x x u ?=0,, 0x l << (1.1.3) 及边值条件 ()()0,,0==t l u t u ， T t ≤≤0 (1.1.4) 假定()x f 和()x ?在相应的区域光滑，并且于()0,0，()0,l 两点满足相容条件,则上述问题有唯一的充分光滑的解。 1.2抛物线扩散方程的求解 下面考虑如下热传导方程 22()(0.)(,)0(,0)()u u a f x t x u t u L t u x x ????=+????? ==??=??? (1.2.1) 其中，0x l <<，T t ≤≤0,a (常数)是扩散系数。 取N l h = 为空间步长，M T =τ为时间步长，其中N ，M 是自然数，用两族 平行直线jh x x j ==, ()N j ,,1,0 =和k t t k τ ==, ()M k ,,1,0 =将矩形域 G {}T t l x ≤≤≤≤=0;0分割成矩形网格。其中 (),j k x t 表示网格节点；h G 表示 网格内点（位于开矩形G 中的网格节点）的集合；h G 表示位于闭矩形G 中的网格节点的集合；h Γ表示h G -h G 网格边界点的集合。 k j u 表示定义在网点(),j k x t 处的待求近似解，N j ≤≤0，M k ≤≤0。 现在对方程进行差分近似： (一) 向前差分格式 =-+τ k j k j u u 111 2 2(())k k k j j j j j j u u u a f f f x h +--++= (1.2.2) ()j j j x u ??==0， k u 0=k N u =0 (1.2.3) 计算后得： 111(12)k k k k j j j j j u ru r u ru f τ++-=+-++ (1.2.4) 其中，2 a r h τ = ，1,,1,0-=N j ，1,,1,0-=M k 。 显然，这是一个四点显示格式，每一层各个节点上的值是通过一个方程组求解到的。方程组如下： 1000 121011000 232121000 3432310001121(12)(12)(12)(12)N N N N N u ru r u ru f u ru r u ru f u ru r u ru f u ru r u ru f ττττ----?=+-++?=+-++??=+-++? ???=+-++? (1.2.5) 若记 () T k N k k k u u u 1 21,,,-= u ，()()()()T N x x x 121,,,-=???? ，()()()()T N x f x f x f 121,,,-=τττ f 则显格式(1.2.4)可写成向量形式 10 ,0,1,,1 k k k M φ +?=+=-?=? u Au f u (1.2.6) 其中 《微分方程数值解法》 【摘要】自然界与工程技术中的很多现象，可以归结为微分方程定解问题。其中，常微分方程求解是微分方程的重要基础内容。但是，对于许多的微分方程，往往很难得到甚至不存在精确的解析表达式，这时候，数值解提供了一个很好的解决思路。，针对于此，本文对常微分方程数值解法进行了简单研究，主要讨论了一些常用的数值解法，如欧拉法、改进的欧拉法、Runge —Kutta 方法、Adams 预估校正法以及勒让德谱方法等，通过具体的算例，结合MA TLAB 求解画图，初步给出了一般常微分方程数值解法的求解过程。同时，通过对各种方法的误差分析，让大家对各种方法的特点和适用范围有一个直观的感受。 【关键词】 常微分方程 数值解法 MA TLAB 误差分析 引言 在我国高校，《微分方程数值解法》作为对数学基础知识要求较高且应用非常广泛的一门课程，不仅 在数学专业，其他的理工科专业的本科及研究生教育中开设这门课程．近四十年来，《微分方程数值解法》不论在理论上还是在方法上都获得了很大的发展．同时，由于微分方程是描述物理、化学和生物现象的数学模型基础，且它的一些最新应用已经扩展到经济、金融预测、图像处理及其他领域 在实际应用中，通过相应的微分方程模型解决具体问题，采用数值方法求得方程的近似解，使具体问题迎刃而解。 2 欧拉法和改进的欧拉法 2.1 欧拉法 2.1.1 欧拉法介绍 首先，我们考虑如下的一阶常微分方程初值问题 ???==0 0)() ,('y x y y x f y （2--1） 事实上，对于更复杂的常微分方程组或者高阶常微分方程，只需要将x 看做向量，（2--1）就成了一个一阶常微分方程组，而高阶常微分方程也可以通过降阶化成一个一阶常微分方程组。 欧拉方法是解常微分方程初值问题最简单最古老的一种数值方法，其基本思路就是把（2--1）中的导数项'y 用差商逼近，从而将一个微分方程转化为一个代数方程，以便求解。 设在[]b a ,中取等距节点h ，因为在节点n x 点上，由（2--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年就多层渗流方程组合系统提出并分析了迎风分数步长差分方法，以上均是讨论的线性对流扩散问题，胡建伟等通过引入质量集中算子，构造并分析了一类基于三角网格的质量集中型的部分有限元方法处理线性和非线性对流扩散问 扩散方程 扩散方程稳态扩散与非稳态扩散 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第一定律） ，，， （即第二个面的扩散通量为第一个面注入的溶质与在这一段距离内溶质浓度变化引起的扩散通 1. 课本2p 有证明 2. 课本812,p p 有说明 3. 课本1520,p p 有说明 4. Rit2法，设n u 是u 的n 维子空间，12,...n ???是n u 的一组基底，n u 中的任一元素n u 可 表为1n n i i i u c ?==∑ ，则,11 11()(,)(,)(,)(,)22j n n n n n n i j i j j i j j J u a u u f u a c c c f ???=== -=-∑∑是12,...n c c c 的二次函数，(,)(,)i j j i a a ????=，令 () 0n j J u c ?=?，从而得到12,...n c c c 满足1 (,)(,),1,2...n i j i j i a c f j n ???===∑，通过解线性方程组，求的i c ，代入1 n n i i i u c ?==∑， 从而得到近似解n u 的过程称为Rit2法 简而言之，Rit2法：为得到偏微分方程的有穷维解，构造了一个近似解，1 n n i i i u c ?== ∑， 利用,11 11()(,)(,)(,)(,)22j n n n n n n i j i j j i j j J u a u u f u a c c c f ???===-=-∑∑确定i c ，求得近似解n u 的过程 Galerkin 法：为求得1 n n i i i u c ? == ∑形式的近似解，在系数i c 使n u 关于n V u ∈，满足(,)(,) n a u V f V =，对任 意 n V u ∈或（取 ,1j V j n ?=≤≤） 1 (,)(,),1,2...n i j i j i a c f j n ???===∑的情况下确定i c ，从而得到近似解1 n n i i i u c ?==∑的过程称 Galerkin 法为 Rit2-Galerkin 法方程： 1 (,)(,)n i j i j i a c f ???==∑ 5. 有限元法：将偏微分方程转化为变分形式，选定单元的形状，对求解域作剖分，进而构 造基函数或单元形状函数，形成有限元空间，将偏微分方程转化成了有限元方程，利用 有效的有限元方程的解法，给出偏微分方程近似解的过程称为有限元法。 6. 解：对求解区间进行网格剖分，节点01......i n a x x x x b =<<<<=得到相邻节点1,i i x x - 对流扩散方程有限差分方法 求解对流扩散方程的差分格式有很多种，在本节中将介绍以下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 数学模型 本作业从最简单的模型方程，即一维、稳态、无源项的对流扩散方程出发，方程如下： 22, 02f f f U D x t x x ???+=≤≤??? (1) 初始条件 (),0sin(2)f x t A kx π== (2) 解析解 ()()()224,sin 2Dk t f x t e A k x Ut ππ-=- (3) 式中，1,0.05,0.5,1U D A k ==== 函数(3)描述的是一个衰减波的图像，如图1所示 t=0 t=0.5 t=1 图1 函数()()()224,sin 2Dk t f x t e k x Ut ππ-=- 的图像（U=1，D=0.05，k=1） 2 数值解法 2.1 数值误差分析 在网格点(),i n 上差分方程的数值解n i f 偏离该点上相应的偏微分方程的精确解 (),f i n 的值，称为网格节点上的数值误差。 当取定网格节点数21N =时，观察差分方程的解与微分方程的解在不同时间步长下的趋近程度，其中时间步长分别取值0.05,0.025,0.0125,0.0005t ?=。 （a ）21,0.05N t =?= （b ）21,0.025N t =?= （c ）21,0.0125N t =?= （d ）201,0.0005N t =?= 图2 数值误差随步长的变化情况 从图2的(a)~(d)可以定性的看出，数值误差与步长的大小有关。在满足稳定性条件的前提下，数值误差随着时间步长的减小而减小，同时，图（d ）表示增大网格的分辨率也有助于减小网格误差。 为了对数值误差有一个定量的认识，接下来取定时间步长为0.0005t ?=，分别算出 11,21,41,61,81,101,121,161N =时，指标E =1所示。 表1 不同网格节点数下指标E 的值 基于Peclet 数判别法的一维对流扩散方程分类研究 摘要：采用Peclet 数的绝对值大小来判别一维对流扩散方程为对流占优型或是扩散占优型方程，运用三种隐式差分格式—中心隐式格式、对流C-N 型格式和扩散C-N 格式，对不同Peclet 数的算例进行离散和求解。然后，将计算区域中所有节点的解析解与数值解表示成矩阵形式，并求解出它们的矩阵2范数之后作比较，两者越接近则代表差分格式精度越高。通过比较得出了当方程Peclet 数的绝对值小于等于0.5时，方程为扩散占优型方程。在离散方法选取方面，针对扩散项的离散可以采用更高精度的差分格式，如扩散C-N 格式；当Peclet 数的绝对值大于等于20时，方程为对流占优型方程。此时，针对对流项可以采用更高精度的差分格式，如对流C-N 格式；当Peclet 数的绝对值介于0.5与20之间时，无法用Peclet 数判断方程类型,不过可以选择折衷的离散格式减小误差，如中心隐式格式。 关键字：一维对流扩散方程 Peclet 数判别法 有限差分方法 数值模拟 MR(2010)主题分类号：39A14；65M06 中图分类号：O242.2 文献标识码： A 1.引言 一维对流扩散方程是描述流体流动和传热问题的一类线性化模型方程。土壤、大气等多孔介质中水、盐分、温度以及污染物质的对流扩散问题都会遇到此类方程。在一维对流扩散方程的求解过程中，反映流体对流和扩散两种物理作用的分别是对流项和扩散项。所以，根据方程中对流项还是扩散项占主导作用，通常可将方程分为对流占优型和扩散占优型两类方程。然而，要想得到精确度较高的数值结果，这两种类型方程的离散方法不能采用相同的离散格式。因此，需要有一种判别方法来判断方程的类型，关于对流占优型和扩散占优型方程的判别方法一直是近年来研究的热点问题。这对研究不同类型的方程使用合适的差分格式进行离散具有实际的意义。由于Peclet 数的绝对值表示了对流作用相对扩散作用的大小，即绝 大，扩散所起的作用就可以忽略。反之，当Peclet 数为零时，方程就为纯扩散方程。本文选用一维定解非稳态对流扩散方程为例，通过考察Peclet 数的绝对值大小来对方程进行分类，方程一般形式如下： 2(,),,0 122(1)(,0)()(,)(),(,)()12(,) u u u a f x t x x x t t x x u x g x u x t t u x t t u u x t υ?φ???? ?? ?? ????+=+≤≤≥???==== 其中a 和υ分别代表对流项系数和扩散项系数。假定求解区间长度为s ， Peclet 数的绝对值 剖面二维非恒定悬移质泥沙扩散方程的数值方法< 1 引言 数学模拟方法正在成为研究河流泥沙问题的重要手段。目前，一维数学模型发展较成熟，已广泛应用于模拟长河段的长期变形，但它只能给出河段平均冲淤深度的沿程变化，如需了解短河段的河床变形细节，则要采用二维以至三维数学模型。不论是一维数学模型还是平面二维维数学模型，都不能反映含沙量沿垂线的分布状况，并忽略了含沙量沿垂线分布对垂线平均含沙量变化过程的影响。要解决这类问题，必须建立剖面二维数学模型。这种模型主要通过解剖面二维泥沙扩散方程来研究悬移质泥沙沿水深的分布及含沙量的变化过程，对水电站进口和其它引水工程的引水口高程的确定都能提供较好的数值模拟。 泥沙扩散方程实际上是一个变系数的二阶线性偏微分方程，这样的方程在各种复杂边界条件下求解是极为困难的。求扩散方程的解析解在数学上存在着难以克服的困难，往往只能通过对方程的简化，才能得到一些简单边界条件下的解析解，在这方面， A.A.Kalinske 、野满隆治、 W.E.Dobbins 、俞维强、张启舜、韦直林 [ ２ ] 等都做了有益的尝试；求扩散方程的数值解曾经因为缺乏高效率的计算工具而难以实现，直到 60 年代后，随着计算机的广泛应用，在各种复杂边界条件下求扩散方程的数值解不但成为可能，而且得到迅速的发展，在这方面，曹志先、崔侠 [ ４ ] 等做了大量工作，取得了很多成果。 数值方法相对于解析方法在求解偏微分方程上有着明显的优势，即简单灵活、计算方便快捷，但要寻找一种精度高、稳定性好、计算方便的差分格式也并非易事。本文拟在前人研究的基础上着重讨论剖面二维泥沙扩散方程的数值解问题，希望能提供一种精度高、稳定性好、计算方便的数值解。 2 基本方程 剖面二维泥沙扩散方程的形式为 微分方程的分类及其数值解法 微分方程的分类： 含有未知函数的导数，如dy/dx=2x 、ds/dt=0.4都是微分方程。 一般的凡是表示未知函数、未知函数的导数与自变量之间的关系的方程，叫做微分方程。未知函数是一元函数的，叫常微分方程；未知函数是多元函数的叫做偏微分方程。微分方程有时也简称方程。 一、常微分方程的数值解法： 1、Euler 法： 00d (,), (1.1)d (), (1.2) y f x y x y x y ?=???=? 001 (),(,),0,1,,1n n n n y y x y y hf x y n N +=??=+=-? (1.4) 其中0,n b a x x nh h N -=+=. 用(1.4)求解(1.1)的方法称为Euler 方法。 后退Euler 公式???+==+++),,(),(111 00n n n n y x hf y y x y y 梯形方法公式 )].,(),([2 111+++++=n n n n n n y x f y x f h y y 改进的Euler 方法11(,),(,),1().2p n n n c n n p n p c y y hf x y y y hf x y y y y ++?=+??=+???=+??? 2、Runge-Kutta 方法： p 阶方法 ： 1()O h -=?总体截断误差局部截断误差 二阶Runge-Kutta 方法 ??? ????++==++=+),,(),,(,2212 1211hk y h x f k y x f k k h k h y y n n n n n n 一维稳态对流扩散方程稳定性条件的推导 姓名： 班级：硕5015 学号： 2015/12/15 证明： 一维稳态对流扩散方程： 22u x x φφρ??=Γ?? 采用控制容积积分法，对上图P 控制的容积作积分，取分段线性型线，对均分网格可得下列离散方程： ()()()()()()()()11112222e w e w P E W e w e w w w e e u u u u x x x x φρρφρφρδδδδ??????ΓΓΓΓ+-+=-++????????????????记：()()()()1122e w P e w w e a u u x x ρρδδΓΓ=+-+ ()()12 e E e e a u x ρδΓ=- ()()12w W w w a u x ρδΓ= + 定义通过界面的流量u ρ记为F ,界面上单位面积扩散阻力的倒数x δΓ记为D ,则原式简化为： P P E E W W a a a φφφ=+ 12 E e e a D F =- 12 W w w a D F =+ ()P E W e w a a a F F =++- 令 u x F Pe D ρδ==Γ 则 1111222 E W P Pe Pe φφφ????-++ ? ?????= 当Pe 大于2以后，数值解出现了异常；P φ小于其左右邻点之值，在无源项情 况下是不可能的。因为当2Pe >时系数12 E e e a D F =-小于零，即右边点的通过对流及扩散作用对中间点所产生的影响是负的，这会导致物理上产生不真实的解，所以2u x Pe ρδ=≤Γ 证毕。 一、用有限差分法求解二维扩散方程的初边值问题 该问题的精确解为1 ()2 (,,)x y t u x y t e +-= 二、用下列差分格式编程计算，并比较计算速度、精度、稳定性。 1． 古典显式格式： 1,,1,,1,,1,,1 2 2 22n n n n n n n n j l j l j l j l j l j l j l j l u u u u u u u u h h τ ++-+---+-+= + 2． 交替方向隐式格式（P-R 格式）： 1111 2222 ,,1,,1,,1,,1 22 1111 1111 2222,,1,,1,,1,,122 222 222 n n n n n n n n j l j l j l j l j l j l j l j l n n n n n n n n j l j l j l j l j l j l j l j l u u u u u u u u h h u u u u u u u u h h ττ+++++-+-+++++++++-+-?--+-+?=+?????--+-+?=+??? 3． 局部一维格式： 1111 2222 ,,1,,1,1,,1,221111 11112222 ,,,1,,1,1,,122 2211()222 2211()222 n n n n n n n n j l j l j l j l j l j l j l j l n n n n n n n n j l j l j l j l j l j l j l j l u u u u u u u u h h u u u u u u u u h h ττ+++++-+-++++ +++++-+-?--+-+?=+?????--+-+?=+??? 4． 预测-校正格式： 第一步：1111 4444 ,,1,,1,2 11111 24222 ,,,1,,1222 22 n n n n n j l j l j l j l j l n n n n n j l j l j l j l j l u u u u u h u u u u u h ττ+++++-++++++-?--+?=?????--+?=??? ， 122 ()2221 ()2 11 (1)2211 (1)223,(0,1,01) 2(,,0)(0,1)(0,,),(1,,)(01,01) (,0,),(,1,)(01,01)x y t x y y t y t x t x t u u u e x y t t x y u x y e x y u y t e u y t e y t u x t e u x t e x t +-+-+--+-?????-+=-<<<≤??????=< i.常微分方程初值问题数值解法 常微分方程初值问题的真解可以看成是从给定初始点出发的一条连续曲线。差分法是常微分方程初值问题的主要数值解法，其目的是得到若干个离散点来逼近这条解曲线。有两个基本途径。一个是用离散点上的差商近似替代微商。另一个是先对微分方程积分得到积分方程，再利用离散点作数值积分。 i.1 常微分方程差分法 考虑常微分方程初值问题：求函数()u t 满足 (,), 0du f t u t T dt =<≤ （i.1a ） 0(0)u u = (i.1b) 其中(,)f t u 是定义在区域G : 0t T ≤≤, u <∞上的连续函数，0u 和T 是给定的常数。我们假设(,)f t u 对u 满足Lipschitz 条件，即存在常数L 使得 121212(,)(,), [0,]; ,(,)f t u f t u L u u t T u u -≤-?∈∈-∞∞ （i.2） 这一条件保证了（i.1）的解是适定的，即存在，唯一，而且连续依赖于初值0u 。 通常情况下，（i.1）的精确解不可能用简单的解析表达式给出，只能求近似解。本章讨论常微分方程最常用的近似数值解法-差分方法。先来讨论最简单的Euler 法。为此，首先将求解区域[0,]T 离散化为若干个离散点： 0110N N t t t t T -=<< <<= （i.3） 其中n t hn =，0h >称为步长。 在微积分课程中我们熟知，微商（即导数）是差商的极限。反过来，差商就是微商的近似。在0t t =处，在（i.1a ）中用向前差商 10()()u t u t h -代替微商du dt ，便得 10000()()(,())u t u t hf t u t ε=++ 如果忽略误差项0ε，再换个记号，用i u 代替()i u t 便得到 1000(,)u u hf t u -= 一般地，我们有 1Euler (,), 0,1, ,1n n n n u u hf t u n N +=+=-方法： （i.4） 从(i.1b) 给出的初始值0u 出发，由上式可以依次算出1,,N t t 上的差分解1,,N u u 。 对流扩散方程背景 提出一种隐格式用于解决二维时间依赖的Burgers型方程。迎风单边差分格式被用于对流项离散；对扩散项用二阶中心差分格式离散。我们建立了全隐的数值有限差分格式，分析了无条件稳定性和严格推导了收敛性，在空间是二阶收敛的和时间一阶收敛的。给出数值结果验证理论正确性。 关键词：隐格式，单边差分逼近，Burgers方程，稳定性，收敛阶 对流扩散方程背景 对流扩散方程描述黏性流体的动力学行为，这在许多工程应用中发挥了重要作用。对流占优型扩散方程一般具有对流比扩散的系数大得多的特点，通常数值模拟具有一定难度，因为一方面，扩散系数比传输速度小，并且在另一方面，由于数值扰动容易出现边界层现象。许多格式已用于这些问题的模拟，并有大量成功的数值方法[1-3]。通过离散方法来解决对流扩散问题时，一般运用标准Galerkin有限元方法求解，但此方法会导致非物理特征扰动。为了解决这类缺陷，几种稳定的有限元方法已经在[4]中被提出了。 我们感兴趣的是建立非耗散方法来克服数值扰动，并有鲁棒性和二阶精度，尤其是对Burgers问题。Burgers问题通常被认为非线性流体的流动和扰动的经典模型。在二维非线性的情况下，可以描述对流和扩散的现象，Burgers方程代表一种最基本的非线性模型方程。从一个数值格式出发研究是相当有趣的，因为Burgers已出现在众多的流体方程中[5-7]。并已经由霍普夫-科尔计算出多种组合的初始条件和边界条件下的结果[10,11]。此外，对于非线性Burgers方程的解析解也可以通过Homotop Perturbation法[12]得到。 众所周知，单独的选择一种基本差分格式如中心差分或者迎风格式，来计算纯对流式的方程，扩散项通常只是中心近似。而解决问题的关键在于对流方面构造稳定的离散结构来克服数值扰动。虽然单边差分近似格式已经被提出了30年之久[13]，人们却很少关注他们在计算流动问题。一阶或者二阶单边迎风有限差分 1 引言 2 对流方程及算法介绍 2.1对流方程的概述 对流：是指由于流体的宏观运动，从而使流体各部分之间发生相对位移，冷热流体相互掺混所引起的热量传递过程。 对流仅发生在流体中，对流的同时必伴随有导热现象。 人们研究对流扩散方程，主要的研究对象是流体在流动过程中，流体所携带的某种物质的物理量的变化规律，例如传热过程中温度的变化规律或者溶解于流体中溶质的物质浓度等物理量的变化规律。这些变化通常包括对流、扩散以及由于某种物理或者化学的因素而引起的物理量的自身衰减或增长。 最简单的一维对流扩散方程形如(2-1)式: （2-1） 其中C 是常数，它属于双曲型方程，可以被用来描述流体的运动等物理现象。 2.2水对流现象的简易演示 2.2.1 基本步骤 用两只相同的小烧杯，各装上冷水，再如图1所示插入长短两根吸管，虹吸管由普通化学实验用玻璃管在酒精灯上加热弯成，一根查到被子底部，一根只插入水的表面，再在右杯中滴入几滴墨水并搅拌均匀，现在开始用酒精灯加热左边的烧杯，一段时间后就可以明显的看到染了颜色的水从右杯源源不断的流入左杯，左杯的水源源不断的流入右杯，最后两杯水都变成了墨水的颜色，与此同时用手摸右边的杯子，右边的水也热了起来，这就是冷热水发生了对流的缘故。 2.2.2 实验注意事项 短虹吸管只插入水的表面，不能过深。 玻璃管宜选壁较厚一些的，这样绝热性好一些，效果也好一些。 0=??+??x u C t u 2.2.3 实验原理分析 对左边的水杯用酒精灯加热，水受热密度变小开始上升，右边水杯的 冷水从下边的吸管流向左边的水杯进行补充，左边水杯的热水从上边的吸管流向右边的水杯，这样一会儿两杯水都变成墨水的颜色了[1]。 在冷水里面掺热水也是一样的道理，在不搅拌的情况下，最后水温基本都是一个温度，这就是水的对流，除了水的对流还有刮风是空气的对流，气压高的一方向气压低的一方补充空气，这就形成了对流，就会产生风；还有冬天在家里开空调，形成空气对流，最后整个房间的温度都升了起来。 2.3对流方程及其现有算法 1.针对常系数对流扩散方程，我们利用指数变换， 构造四阶紧致差分格式。 2.针对一维变系数对流扩散方程，将其转化为扩散方程，并构造四阶紧致差分格式。 3.对于常系数二维对流扩散方程，构造出四阶紧致差分方程，以及特殊的变系数 对流扩散方程的四阶紧致差分格式。 4.针对一维常系数对流扩散方程 和一维变系数对流扩散方程，分别构造了几种基于线性和双线性插值 的特征差分格式。 5.针对二维对流扩散方程 ，构造了几种基于线性和双线性插值的特征差分格式。 2.3影响物理量?的三个过程 用),,,(t z y x ??=来表示流体中单位体积的流体所携带的某种物理量，它可以是流体的质量或温度。流体的温度可以用?来表示，流体的密度ρ也可以用? ),(22t x f x u x u a t u +??=??+??ε2 2),(x u x u t x a t u ??=??+??εf y u x u a y u q x u p t u =??+??-??+??+??)(2222) ,()()()(2222y x f y u y q x u x p y u x u =??+??+??+??22x u x u v t u ??=??+??ε) ,()()()(22t x f x x a x u x b x u x c u =??-??+??) ,())(,(),(),(),(222221y x f y u x u y x a y u y x b x u y x b t u y x c =??+??-??+??+?? 包括基本概念，差分格式的构造、截断误差和稳定性，这些内容是贯穿整个教材的主线。解答问题关键在过程，能够显示出你已经掌握了书上的内容，知道了解题方法。这次考试题目的类型：20分的选择题，主要是基本概念的理解，后面有五个大题，包括差分格式的构造、截断误差和稳定性。 习题一 1． 略 2． y y x f -=),(，梯形公式：n n n n n n y h h y y y h y y )121(),(2111+-+=+- =+++，所以0122)1(01])121[()121()121(y h h y h h y h h y h h n h h n n n +--+--+-+=+-+==+-+= ，当0→h 时， x n e y -→。 同理可以证明预报-校正法收敛到微分方程的解. 3． 局部截断误差的推导同欧拉公式; 整体截断误差: ? ++++++-++≤1 ),())(,(11111n n x x n n n n n n n dx y x f x y x f R εε 11)(++-++≤n n n y x y Lh R ε，这里R R n ≤ 而111)(+++-=n n n y x y ε，所以 R Lh n n += -+εε1)1(，不妨设1对流扩散方程
扩散方程稳态扩散与非稳态扩散
抛物形扩散方程的有限差分法及数值实例
微分方程数值解法
对流扩散方程引言
fick定律扩散方程
偏微分方程数值解法答案
对流扩散方程有限差分方法.
一维对流扩散方程的数值解法
基于Peclet数判别法的一维对流扩散方程分类研究
剖面二维非恒定悬移质泥沙扩散方程的数值方法.
微分方程的分类及其数值解法
一维对流扩散方程的稳定性条件推导
二维扩散方程
常微分方程数值解法
对流扩散方程背景
对流方程及算法介绍
微分方程数值解法答案
相关主题
文本预览