徐州工程学院 课程设计报告 课程名称偏微分方程数值解 课题名称对流扩散方程 的迎风格式的推导和求解专业信息与计算科学 班级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. 无限大物体中的扩散
对流扩散方程的定解问题是指物质输运与分子扩散的物理过程和黏性流体流动的数学模型，它可以用来描述河流污染、大气污染、核污染中污染物质的分布，流体的流动和流体中热传导等众多物理现象。关于对流扩散方程的求解很也备受关注，因此寻找一种稳定实用的数值方法有着重要的理论与实际意义。 求解对流扩散方程的数值方法有多种，尤其是对流占优扩散方程，这些方法有迎风有限元法，有限体积法，特征有限体积法，特征有限差分法和特征有限元法，广义差分法，流线扩散法，以及这些方法与传统方法相结合的方法如迎风广义差分法，迎风有限体积法有限体积——有限元法等这些方法数值求解效果较好，及有效的避免了数值震荡，有减少了数值扩散，但是一般计算量偏大 近年，许多研究者进行了更加深入的研究，文献提出了对流扩散方程的特征混合元法，再次基础上，陈掌引入了特征间断混合元方法，还有一些学者将特征线和有限体积法相结合，提出了特征有限体积元方法（非线性和半线性），于此同时迎风有限元也得到较大的发展，胡建伟等人研究了对流扩散问题的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年就多层渗流方程组合系统提出并分析了迎风分数步长差分方法，以上均是讨论的线性对流扩散问题，胡建伟等通过引入质量集中算子，构造并分析了一类基于三角网格的质量集中型的部分有限元方法处理线性和非线性对流扩散问
2013—2014学年第二学期 《Matlab 编程技术》作业 专业班级 石工博13-02 研究方向 油气田开发 姓 名 王壮壮 学 号 B13020075
结合自己研究方向，运用Matlab编写科学计算及可视化或其它相关程序。要求： 1）将要解决的问题交代清楚（数学模型、目标等）； 2）编写的程序的关键语句要有注释说明； 3）程序能顺利运行，运行结果和编写的m文件一并提交； 4）独立完成。
非线性对流扩散方程不同解法稳定性比较 流体力学基本方程组本身就是非线性的对流扩散方程，非线性Burgers 方程就是N-S 方程很好的模型方程，它的一维形式如下： L x x u x u u t u ≤≤??=??+??022μ (1) 边界条件为 ? ? ?====0,,00 u L x u u x (2) 初始条件是任意可以给出的。 我们知道，遇到对流项，我们用迎风格式是绝对没有问题，无论是一阶迎风还是二阶迎风格式都是能够解决非线性对流方程的，如果网格Peclet 数允许的话，中心差分也是可以考虑的。 不过，对于非线性对流，我们来看看另外两个G-S 格式，一个是G-S 型迎风半隐格式，另一个是G-S 型Samarskii 半隐格式，对于任何类型的对流扩散方程，可以收敛到定常解，并且是绝对稳定的，特别适合于解决定常问题。 对于式(1)这两个格式分别为 () 2 11 111111212h u u u R h u u u u u n i n i n i n i n i n i n i n i n i +-+++-+++-+=-+-μτ (3) 21 1 1111112112h u u u R R h u u u u u n i n i n i n i n i n i n i n i n i n i +-+++-+++-???? ??++=-+-μτ (4) 其中 μ 2h u R n i n i = 式(3)就是G-S 型迎风半隐格式，它具有一阶精度，是从一阶迎风格式发展而来的；式(4)是G-S 型Samarskii 半隐格式，具有二阶精度，它是从Samarskii 格式发展而来的。上面说过，它们只适用于求解定常解，因此上标的时层n 可以看作是迭代步，可以说它们没有时间精度，如果想用这两个格式求解非定常解，那可是徒劳的。 对于上两式，我们可以采用迭代法求解，把它们写成迭代式，分别为 ()[]( )( ) ( ) n i n i n i n i n i n i n i n i n i R h u h u u R u u u h u ++++++--= +-++-++142212*********τμτμτ (5)
对流扩散方程有限差分方法 求解对流扩散方程的差分格式有很多种，在本节中将介绍以下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 ?-++????????-??? ?????+????????=τ
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
Subjects ： A property φ is transported by means of convection and diffusion through the one-dimensional domain. The governing equation is ()()d d d u dx dx dx φ ρφ=Γ; boundary conditions are 01φ= at x=0 and 0L φ= at x =L. Using five equally spaced cells and the central differencing scheme for convection and diffusion calculate the distribution of φas a function of x for (i)Case 1:
对流扩散方程背景 提出一种隐格式用于解决二维时间依赖的Burgers型方程。迎风单边差分格式被用于对流项离散；对扩散项用二阶中心差分格式离散。我们建立了全隐的数值有限差分格式，分析了无条件稳定性和严格推导了收敛性，在空间是二阶收敛的和时间一阶收敛的。给出数值结果验证理论正确性。 关键词：隐格式，单边差分逼近，Burgers方程，稳定性，收敛阶 对流扩散方程背景 对流扩散方程描述黏性流体的动力学行为，这在许多工程应用中发挥了重要作用。对流占优型扩散方程一般具有对流比扩散的系数大得多的特点，通常数值模拟具有一定难度，因为一方面，扩散系数比传输速度小，并且在另一方面，由于数值扰动容易出现边界层现象。许多格式已用于这些问题的模拟，并有大量成功的数值方法[1-3]。通过离散方法来解决对流扩散问题时，一般运用标准Galerkin有限元方法求解，但此方法会导致非物理特征扰动。为了解决这类缺陷，几种稳定的有限元方法已经在[4]中被提出了。 我们感兴趣的是建立非耗散方法来克服数值扰动，并有鲁棒性和二阶精度，尤其是对Burgers问题。Burgers问题通常被认为非线性流体的流动和扰动的经典模型。在二维非线性的情况下，可以描述对流和扩散的现象，Burgers方程代表一种最基本的非线性模型方程。从一个数值格式出发研究是相当有趣的，因为Burgers已出现在众多的流体方程中[5-7]。并已经由霍普夫-科尔计算出多种组合的初始条件和边界条件下的结果[10,11]。此外，对于非线性Burgers方程的解析解也可以通过Homotop Perturbation法[12]得到。 众所周知，单独的选择一种基本差分格式如中心差分或者迎风格式，来计算纯对流式的方程，扩散项通常只是中心近似。而解决问题的关键在于对流方面构造稳定的离散结构来克服数值扰动。虽然单边差分近似格式已经被提出了30年之久[13]，人们却很少关注他们在计算流动问题。一阶或者二阶单边迎风有限差分
A
对流扩散方程的求解 对流扩散问题的有效数值解法一直是计算数学中重要的研究内容，求解对流扩散方程的数值方法主要是有限差分法(FDM)、有限元法(FEM)、有限体积法(FVM)、有限解析法(FAM)、边界元法(BEM)、谱方法(SM) 等多种方法。但是对于对流占优问题，用通常的差分法或有限元法进行求解将出现数值震荡。 为了克服数值震荡，80年代，J.Douglas,Jr.和T.F.Russell 等提出特征修正技术求解对流扩散占优的对流扩散问题，与其它方法相结合，提出了特征有限元方法、特征有限差分方法、特征混合元方法；T.J.Hughes和A.Brooks提出过一种沿流线方向附加人工黏性的间断有限元法，称为流线扩散方法(SDM)。有限差分法、有限元法、有限体积法是工程应用中的主要方法。 对流扩散方程的特点 对流扩散方程右端第一项为扩散项，左端第二项则是对流项。由于其方程本身的特点，给建立准确有效的数值求解方法带来一定的困难。对流和扩散给流体中由流体携带的某种物理量的变化过程，可以通过一个无量纲的特征参数(Peclet数)来描述，Peclet数Pe的定义为：Pe=|ν|L/D。这里v是来流速度，L是特征长度，D是物质的扩散系数。如果Pe数较小，即对流效应相对较弱，这类问题中，扩散占主导地位，方程是椭圆型或抛物线型；
如果Pe数较大，即溶质分子的扩散相对于流体速度而言是缓慢的，这类问题中，对流占优，方程具有双曲型方程的特点。 对于对流占优问题的求解，采用常规的Galerkin有限元方法，为了避免求解结果产生数值振荡，获得稳定解，则应使每个单元的局部Peclet数，Peh=|ν|h/D≤2，这里h为单元的最大尺寸，|v|为单元中的最大速度分量值。因此，用本文方法求解对流占优对流扩散问题，要得到稳定解，则要通过加密有限元网格来实现。