浅水方程的二阶广义迎风差分方法

南京师范大学 硕士学位论文 浅水方程的二阶广义迎风差分方法 姓名：李春兰 申请学位级别：硕士 专业：数学；计算数学 指导教师：张志跃 2009-05-15 摘 要 摘摘摘 要要要 广义迎风差分方法, 结合了有限差分方法和有限元方法的特 点,与当前求解计算流体力学常用的有限体积数值解法相接近． 本文第一部分即引言主要介绍了浅水方程的相关及其发展状 况,并对广义迎风差分方法作简略介绍．第二部分主要介绍了一 维浅水方程的广义迎风差分格式并对数值格式给出 L2误差估计, 误差阶为 O(∆t + h3/2),并给出数值例子验证该格式的有效性;第 三部分将一维广义迎风差分格式推广到二维浅水方程,给出了二 维浅水方程的广义迎风差分格式及其 L2误差估计,得到误差阶为 O(∆t+ h3/2);最后给出数值例子说明该格式的有效性．在一维和 二维的数值实验中,都采取了二次多项式函数构造有限元函数空 间． 关键词: 浅水方程；广义迎风差分方法；L2误差估计 – iii – Abstract Abstract In this thesis, we apply the generalized upwind difference method to analyze the problem of shallow water model numerically. The gen- eralized upwind difference method has the advantage of finite differ- ence methods and finite element methods, and it is similar to the fi- nite volume method. There are three parts in this thesis. In part 1, we give some background information about the shallow water model and the generalized difference method; In part 2, we present a full dis- crete generalized upwind difference scheme and its L2 error estimates for the 1D shallow water equations and we showed the error estimate is O(∆t+ h3/2). At last, we give numerical experiments that demonstrate that our method is efficient; In part 3, we present some preliminary re- search about the 2D shallow water equations, and show that the error estimates is O(∆t + h3/2). In this paper, the polynomial family of two degree is used to construct a finite element space in the 1D and 2D nu- merical tests. Keywords: Shallow water equations; generalized upwind difference method; L2-error estimate. – iv – 学位论文独创性声明 本人郑重声明: 1、坚持以 “求实、创新”的科学精神从事研究工作。 2、本论文是我个人在导师指导下进行的研究工作和取得的研究成果。 3、本论文中除引文外，所有实验、数据和有关材料均是真实的。 4、本论文中除引文和致谢的内容外，不包含其他人或其它机构已经发表 或撰写过的研究成果。 5、其他同志对本研究所做的贡献均已在论文中作了声明并表示了谢意。 作者签名： 日 期： 学位论文使用授权声明 本人完全了解南京师范大学有关保留、使用学位论文的，学校有权 保留学位论文并向国家主管部门或其指定机构送交论文的电子版和纸质版； 有权将学位论文用于非赢利目的的少量复制并允许论文进入学校图书馆被查 阅；有权将学位论文的内容编入有关数据库进行检索；有权将学位论文的标 题和摘要汇编出版。保密的学位论文在解密后适用本规定。 作者签名： 日 期： 前 言 前前前 言言言 许多工程中的流动问题常常可以用浅水流动模型来描述,如近海,湖泊,水 库及河口中的流动 [2,5,6]. 其中溃坝问题是大家最关心的问题.大坝溃决时，库 区蓄水将突然下泄,造成下游水位陡涨和库水位陡降,这种非恒定流,包括坝 址上、下游河道的水流波动称为溃坝波 [1],它是一种断波,水流的一些物理参 数在断面上产生间断. 溃坝洪水与其他类型的洪水相比,具有其特殊性,影响 因素复杂,产生机理目前仍不完全清楚,溃决过程中溃决模式 (瞬溃、渐溃)、 溃口形状、决口泥沙输运、边坡稳定性、决口的扩展过程、下游河床的水文 水力因素等都会影响溃坝洪水波的传播.溃坝洪水理论上属于洪水水力学的 范畴,但溃坝洪水波往往发生于大坝和堤防的突然溃决时段,其流量、水位在 极短的时间内急剧变化,波形在传播过程中发生间断,流态急剧变化 [2],其分 析计算需要采取特殊的方法作特殊的处理. 浅水方程的数值模拟在计算机出现的早期就开始得到应用,在本世纪四 十年代后期, Charney等人就利用它做了大气的数值模拟 [5]. Hansen数值模拟 了海洋流动 [6]. 至今,它的数值研究已取得了很大的发展 [7]. 浅水方程在数学 上属双曲方程组,浅水方程的求解 [12,19−25] 是计算数学,计算流体力学的一个 重要研究内容.现今发展起来的方法主要有有限差分法 [8,9],特征法 [10,19],有限 元法 [11]等. 特征线法最初源自Massau对浅水方程的图解积分,在 x− t平面上绘制特 征线,在其交点上确定因变量来求解. Haetree发展了特征有限差分法 [26],林秉 南在 20世纪 40-50年代曾对浅水方程的特征理论有过深入的研究 [1]. 特征线 法不足之处在于求解格式复杂,尤其对高维问题更为烦琐. 因此,目前很少用 于实际的数值计算,多作为理解其他数值方法的基础. 有限差分法是数值计算 最早采用且至今仍被广泛应用的基本方法. Chow等 1973年提出了基于有限 差分法、可处理间断与支流入流的非恒定流模型 [27]. Sakkas等 [28] 于 1978年 建立的模型将特征线法用于移动网格,可模拟缓流与急流、连续流与间断流, 且可处理复杂地形. 槐文信等 [29] 用贴体坐标变换与剖开算子法,在计算平面 的交错网格上对控制方程进行隐格式离散,进行了洪潮遭遇情况下的水动力 学计算.胡四一等 [30] 对高性能差分格式在明渠非恒定流计算中的应用进行了 有益的探讨. 流体力学的有限元法借鉴于固体力学的有限元法,在 20世纪 60 – 1 – 前 言 年代开始应用于流体数值计算,其理论基础是变分原理和剖分、插值技术. 常 见的有限元计算方法有直接法、变分法、加权余量法及能量平衡法等,有限 元法适应性强、精度较高,原来存在的计算格式复杂、计算量较大、大型系 数矩阵求解困难等问题,都随着算法的优化和计算机性能的不断改进已不再 是难以解决的问题.但是有限元法对连续性要求限制了其处理流体参数大变 形的能力,捕捉锐利波形比较困难,针对这种情况,发展了对间断问题处理能 力较强的间断有限元法 (discontinuous finite element method)[31−32]. 有限体积法的基本思路是把计算区域离散为若干点,以这些点为中心,将 整个计算区域划分为若干互相连接互不重叠的控制体.在有限体积法的计算 中,将基本方程对每个控制体积分,得到一组以计算节点上的物理量为未知数 的代数方程组,同时方程组离散沿坐标方向进行,形成的离散网格和离散方程 与有限差分法相似,由于采用守恒型的微分方程并对每一计算体进行质量和 动量守恒形式的离散,使得微分方程包含的守恒性质在每一个控制体都得到 满足,保持各单元两侧相邻控制体的计算数值通量相等,整个计算区域上也能 保持守恒性质. 有限体积法兼有有限差分法的物理概念清晰和有限元法的适 应不规则网格、复杂边界条件及计算精度较高的特点,因此在流体动力学的 数值计算领域有着很大的发展潜力. 国外对有限体积法的研究较早,也较完善 [17,18,33]. 一些大型流体计算软件 也运用了有限体积法 (如 Fluent, Star-CD等)[34] 进行计算,计算结果较为理想. 同时,国外将有限差分法的 TVD概念引入有限体积法,构造了许多高分辨率 的有限体积法 TVD格式 [35]. Zhou等 [36]1994年建立了以无结构三角形或四边 形网格、采用 Osher格式的有限体积法模拟二维非恒定流的模型,用于 Florida 的 Kissimmee河流域的水流模拟. 国内,王船海等 [27] 于 1996年将流域径流概 化为调蓄单元的零维模型、河道水流的一维模型、洪泛区的二维模型,采用 直接坐标系下非均匀网格的有限体积法进行流域的溃坝洪水模拟;胡四一等 [24] 在研究溃坝波的数值模拟上参考了空气动力学的某些概念和计算方法,并 根据圣维南方程组在物理和数学意义上与空气动力学中的不可压缩气体的无 黏运动方程组类似以及两者之间的差别发展了一系列高性能、适合非恒定流 (溃坝流)模拟的有限体积法的计算格式;谭维炎 [16]于 1996年用有限体积法进 行了长江中下游及洞庭湖防洪系统的洪水演进模拟;谭维炎等 [37]提出一种有 限体积法高性能格式 (Osher),在非结构网格单元引入逆风概念,从而实现了跨 单元界面法向数值通量的逆风分解,并对长江口二维浅水流动进行了模拟. 同 时还有不少学者提出了诸如运动界面追踪法、有限分析法、边界元法等方法. – 2 – 前 言 本文对浅水方程运用广义差分法 [11](又称有限体积元法)进行求解. 此方 法是由李荣华教授于 1982年提出,综合了 FDM和 FEM的优点,被认为是两者 之间的桥梁. 本文在格式构造上采用迎风格式. 本文下面的安排是分别对一维 和二维浅水方程构造广义迎风差分格式,给出误差估计,并通过数值例子说明 该格式的有效性. – 3 – 第 1章 一维浅水方程的二阶广义迎风差分方法 第第第 1章章章 一一一维维维浅浅浅水水水方方方程程程的的的二二二阶阶阶广广广义义义迎迎迎风风风差差差分分分方方方法法法 在本章里研究一维浅水方程,将一维浅水方程转化为特征方程组,给出其 广义迎风差分格式,并对所给的格式进行误差估计,最后用数值例子说明所给 格式的有效性. 1.1 控控控制制制方方方程程程 考虑如下方程 [14,15]: ∂H ∂t + ∂(Hu) ∂x = s, (1.1.1a) ∂(Hu) ∂t + ∂(Hu2 + gH2/2) ∂x = −gHB′(x). (1.1.1b) 其中 x ∈ [a, b], t ∈ [0, T ]. H 为水深; u为 x方向的深度平均速度分量; g为 重力加速度; s为源项. B 代表底部地形 (B = 0时为平底情形). 浅水方程属于 双曲方程,解的基本特点是具有特征构造. 令 p = u + 2c, q = u − 2c, a1 = u + c, a2 = u − c. 则 u = 12(p + q), H = 1 16g (p− q)2,方程 (1.1.1)的特征方程形式为： ∂p ∂t + a1 ∂p ∂x = s1, (1.1.2a) ∂q ∂t + a2 ∂q ∂x = s2, (1.1.2b) 其中 c = √gH . 1.2 广广广义义义迎迎迎风风风差差差分分分格格格式式式 设 I = [a, b],下面对区间 I = [a, b]做剖分 Th,节点为 a = x0 < x 1 2 < x1 < x 3 2 < · · · < xn− 1 2 < xn = b. – 4 – 第 1章 一维浅水方程的二阶广义迎风差分方法 其中 xi− 1 2 = (xi+xi−1)2 .单元的长度为 hi = xi − xi−1,记 h = max1≤i≤nhi. 并且假定 剖分满足正则性条件: hi ≥ µh,其中 µ为正常数, i = 1, . . . , n. 再做 Th的对偶剖分 T ∗h ,节点为 a = x0 < x 1 4 < x 3 4 < · · · < xn− 3 4 < xn− 1 4 < xn = b. T ∗h = {I∗i , I∗i− 1 2 : I∗i = [xi− 1 4 , xi+ 1 4 ]; I∗ i− 1 2 = [xi− 3 4 , xi− 1 4 ], i = 1, 2, . . . , n − 1, I∗0 = [x0, x 1 4 ], I∗n = [xn− 1 4 , xn]},其中 xi− k 4 = xi − k4hi(k = 1, 3, i = 1, 2, . . . , n). 取试探函数空间为 Uh = {φi, φi− 1 2 , i = 1, 2, . . . , n}. 其中: φi(x) =  (2|x−xi|hi − 1)( |x−xi| hi − 1), xi−1 ≤ x ≤ xi, (2|x−xi|hi+1 − 1)( |x−xi| hi+1 − 1), xi ≤ x ≤ xi+1, 0, 其他. φi− 1 2 (x) = { 4(1− x−xi−1hi )( x−xi−1 hi ), xi−1 ≤ x ≤ xi, 0, 其他. 对应于 T ∗h 的检验函数空间取作分片常数空间 Vh = {ψi(x), ψi− 1 2 (x), i = 1, 2, . . . , n},其中 ψi(x) = { 1, xi− 1 4 ≤ x ≤ xi+ 1 4 , 0, 其他. ψi− 1 2 (x) = { 1, xi− 3 4 ≤ x ≤ xi− 1 4 , 0, 其他. 由于 Vh 的对偶单元边界不连续,不能在整个区域 I 上应用 Galerkin有限 元法,但可以在单个的对偶单元 I∗j , I∗j− 1 2 上应用. 因此,求 ph ∈ Uh满足:  ∫ I∗j [∂ph∂t + a1 ∂ph ∂x ]vhdx = ∫ I∗j s1vhdx,∫ I∗ j− 12 [∂ph∂t + a1 ∂ph ∂x ]vhdx = ∫ I∗ j− 12 s1vhdx, ∀vh ∈ Vh. (1.2.1) 应用 Green公式 – 5 – 第 1章 一维浅水方程的二阶广义迎风差分方法  ∫ I∗j a1 ∂ph ∂x vhdx = (a1phvh)∂I∗j − ∫ I∗j ph ∂(a1vh) ∂x dx,∫ I∗ j− 12 a1 ∂ph ∂x vhdx = (a1phvh)∂I∗j− 12 − ∫I∗ j− 12 ph ∂(a1vh) ∂x dx, (1.2.2) 因此我们可以把 (1.2.1)写为: ∫ I∗j ∂ph ∂t vhdx− ∫ I∗j ph ∂(a1vh) ∂x dx+ (a1phvh)∂I∗j = ∫ I∗j s1vhdx,∫ I∗ j− 12 ∂ph ∂t vhdx− ∫ I∗ j− 12 ph ∂(a1vh) ∂x dx+ (a1phvh)∂I∗j− 12 = ∫ I∗ j− 12 s1vhdx, (1.2.3) 区间 I 的边界为 ∂I ,定义: (∂I)− = a, (∂I)+ = b, ai ≥ 0; (∂I)+ = a, (∂I)− = b, ai < 0; i = 1, 2. 同样定义 (∂I∗i )±.对 ∀x ∈ ∂I∗i ,记 p+h (x) =  limx′→x ph(x ′), x ∈ (∂I∗i )−, x′∈¯I∗i , lim x′→x ph(x ′), x ∈ (∂I∗i )+, x′ ∈ I∗i . p−h (x) =  limx′→x ph(x ′), x ∈ (∂I∗i )−, x′ ∈ I∗i , lim x′→x ph(x ′), x ∈ (∂I∗i )+, x′∈¯I∗i . 他们分别对应 x ∈ ∂I∗i 时 ph(x)的迎风和顺风值.类似于经典迎风格式,我 们用 p+h 来代替 (1.2.3)式左边的边界上的值 ph得:  ∫ I∗j ∂ph ∂t vhdx− ∫ I∗j ph ∂(a1vh) ∂x dx+ (a1p + h vh)∂I∗j = ∫ I∗j s1vhdx,∫ I∗ j− 12 ∂ph ∂t vhdx− ∫ I∗ j− 12 ph ∂(a1vh) ∂x dx+ (a1p + h vh)∂I∗j− 12 = ∫ I∗ j− 12 s1vhdx. ∀vh ∈ Vh. (1.2.4) – 6 – 第 1章 一维浅水方程的二阶广义迎风差分方法 由 (1.2.2)有 − ∫ I∗j ph ∂(a1vh) ∂x dx = ∫ I∗j a1 ∂ph ∂x vhdx− (a1phvh)∂I∗j , − ∫ I∗ j− 12 ph ∂(a1vh) ∂x dx = ∫ I∗ j− 12 a1 ∂ph ∂x vhdx− (a1phvh)∂I∗ j− 12 . 将上式代入 (1.2.4)就得到 (1.1.2a)的半离散广义迎风差分格式: ∫ I∗j ∂ph ∂t vhdx+ ∫ I∗j a1 ∂ph ∂x vhdx+ (a1[ph]vh)(∂I∗j )− = ∫ I∗j s1vhdx,∫ I∗ j− 12 ∂ph ∂t vhdx+ ∫ I∗ j− 12 a1 ∂ph ∂x vhdx+ (a1[ph]vh)(∂I∗j− 12 )− = ∫ I∗ j− 12 s1vhdx. ∀vh ∈ Vh. (1.2.5a) 同理得到 (1.1.2b)的半离散广义迎风差分格式： ∫ I∗j ∂qh ∂t vhdx+ ∫ I∗j a2 ∂qh ∂x vhdx+ (a2[qh]vh)(∂I∗j )− = ∫ I∗j s2vhdx,∫ I∗ j− 12 ∂qh ∂t vhdx+ ∫ I∗ j− 12 a2 ∂qh ∂x vhdx+ (a2[qh]vh)(∂I∗j− 12 )− = ∫ I∗ j− 12 s2vhdx. ∀vh ∈ Vh. (1.2.5b) 其中 [ph] = p+h − p−h 是经过对偶剖分边界的跳跃. 同理 [qh] = q+h − q−h 是经 过对偶剖分边界的跳跃. 初边值条件为{ (ph, qh)(x, t) = (0, 0), x ∈ (∂I)−, (ph, qh)(x, 0) = Φh(x), x ∈ I. (1.2.5c) 对时间导数可以采用不同的差商方法进一步离散,可以得到相应的全离 散格式,下面用向后差分对时间离散得:  ∫ I∗j pnh−pn−1h 4t vhdx+ ∫ I∗j a1 ∂ph ∂x vhdx+ (a1[ph]vh)(∂I∗j )− = ∫ I∗j s1vhdx,∫ I∗ j− 12 pnh−pn−1h 4t vhdx+ ∫ I∗ j− 12 a1 ∂ph ∂x vhdx+ (a1[ph]vh)(∂I∗j− 12 )− = ∫ I∗ j− 12 s1dx. ∀vh ∈ Vh. (1.2.6a) – 7 – 第 1章 一维浅水方程的二阶广义迎风差分方法  ∫ I∗j qnh−qn−1h 4t vhdx+ ∫ I∗j a2 ∂qh ∂x vhdx+ (a2[qh]vh)(∂I∗j )− = ∫ I∗j s2vhdx,∫ I∗ j− 12 qnh−qn−1h 4t vhdx+ ∫ I∗ j− 12 a2 ∂qh ∂x vhdx+ (a2[qh]vh)(∂I∗j− 12 )− = ∫ I∗ j− 12 s2vhdx. ∀vh ∈ Vh. (1.2.6b) 其中4t = tn − tn−1为时间步长. 令 p¯ = pe−c1t,其中c1 = c0 + 1 2 sup x∈Ω |diva1(x)|, c0 > 0. q¯ = qe−c2t,其中c2 = c0 + 1 2 sup x∈Ω |diva2(x)|, c0 > 0. 定义双线性形式 A1(p¯, ψj) = n∑ j=1 {∫I∗j a1 ∂p¯∂xdx+ (a1[p¯])(∂I∗j )− + ∫I∗j c1p¯dx}, A1(p¯, ψj− 1 2 ) = n∑ j=1 {∫I∗ j− 12 a1 ∂p¯ ∂xdx+ (a1[p¯])(∂I∗j− 12 )− + ∫ I∗ j− 12 c1p¯dx}. (1.2.7)  A2(q¯, ψj) = n∑ j=1 {∫I∗j a2 ∂q¯∂xdx+ (a2[q¯])(∂I∗j )− + ∫I∗j c2q¯dx}, A2(q¯, ψj− 1 2 ) = n∑ j=1 {∫I∗ j− 12 a2 ∂q¯ ∂xdx+ (a2[q¯])(∂I∗j− 12 )− + ∫ I∗ j− 12 c2q¯dx}. (1.2.8) 因此,我们可以把 (1.2.5)写为等价形式:{ (p¯h,t, ψj) + A1(p¯h, ψj) = (s¯1, ψj), (p¯h,t, ψj− 1 2 ) + A1(p¯h, ψj− 1 2 ) = (s¯1, ψj− 1 2 ). j = 1, 2, · · · , n. (1.2.9a) { (q¯h,t, ψj) + A2(q¯h, ψj) = (s¯2, ψj), (q¯h,t, ψj− 1 2 ) + A2(q¯h, ψj− 1 2 ) = (s¯2, ψj− 1 2 ). j = 1, 2, · · · , n. (1.2.9b) – 8 – 第 1章 一维浅水方程的二阶广义迎风差分方法 记 p¯nh = p¯h(x, tn), s¯n1 = s¯1(x, tn), tn = n4t.则全离散格式的等价形式为: (p¯nh, vh) +4tA1(p¯nh, vh) = (p¯n−1h +4ts¯n1 , vh), ∀vh ∈ Vh. (1.2.10a) (q¯nh , vh) +4tA2(q¯nh , vh) = (q¯n−1h +4ts¯n2 , vh), ∀vh ∈ Vh. (1.2.10b) 1.3 误误误差差差估估估计计计 引引引理理理 1.1: 对上述定义的双线性形式,有 A1(p¯h, p¯h) ≥ γ0(‖p¯h‖2∂I + ‖p¯h‖20). 其中: γ0 = min(σ1, 12), σ1 = c1 − 12diva1h, ‖p¯h‖20 = (p¯h, p¯h), ‖p¯h‖2∂I = n∑ j=0 {(|a1h|[p¯h]2)(∂I∗j )− + (|a1h|[p¯h]2)(∂I∗j+12 )−}+ 1 2 (|a1h|(p¯+h )2)(∂I∗0 )+ . 即 A1(p¯h, p¯h)是正定的. 证证证明明明: 因 div(avh) = vhdiva+ a · ∇vh,因此我们有 (av2h)∂I∗j = 2 ∫ I∗j vh(a · ∇vh)dx+ ∫ I∗j v2hdivadx. 因此, A1(p¯h, p¯h) = n∑ j=1 { ∫ I∗j a1h ∂p¯h ∂x p¯hdx+ (a1h[p¯h]p¯h)(∂I∗j )− + ∫ I∗j c1p¯ 2 hdx} + n∑ j=1 { ∫ I∗ j− 12 a1h ∂p¯h ∂x p¯hdx+ (a1h[p¯h]p¯h)(∂I∗ j− 12 )− + ∫ I∗ j− 12 c1p¯ 2 hdx} – 9 – 第 1章 一维浅水方程的二阶广义迎风差分方法 = n∑ j=1 {1 2 (a1hp¯ 2 h)∂I∗j − 1 2 ∫ I∗j p¯2hdiva1hdx+ (a1h[p¯h]p¯h)(∂I∗j )− + ∫ I∗j c1p¯ 2 hdx} + n∑ j=1 {1 2 (a1hp¯ 2 h)∂I∗ j− 12 − 1 2 ∫ I∗ j− 12 p¯2hdiva1hdx+ (a1h[p¯h]p¯h)(∂I∗ j− 12 )− + ∫ I∗ j− 12 c1p¯ 2 hdx} = n∑ j=1 {1 2 (a1hp¯ 2 h)(∂I∗j )+ + 1 2 (a1hp¯ 2 h)(∂I∗j )− + [a1h(p¯ + h − p¯−h )p¯h](∂I∗j )− + 1 2 [a1h((p¯ + h ) 2 − 2p¯+h p¯−h + (p¯−h )2)](∂I∗j )− − 1 2 (a1h[p¯h] 2)(∂I∗j )− + ∫ I∗j σ1p¯ 2 hdx} + n∑ j=1 {1 2 (a1hp¯ 2 h)(∂I∗ j− 12 )+ + 1 2 (a1hp¯ 2 h)(∂I∗ j− 12 )− + [a1h(p¯ + h − p¯−h )p¯h](∂I∗ j− 12 )− + 1 2 [a1h((p¯ + h ) 2 − 2p¯+h p¯−h + (p¯−h )2)](∂I∗ j− 12 )− − 1 2 (a1h[p¯h] 2)(∂I∗ j− 12 )− + ∫ I∗ j− 12 σ1p¯ 2 hdx} = 1 2 n∑ j=1 {(a1hp¯2h)(∂I∗j )+ − (a1h[p¯h]2)(∂I∗j )− + (a1h(p¯+h )2)(∂I∗j )− + 2 ∫ I∗j σ1p¯ 2 hdx} + 1 2 n∑ j=1 {(a1hp¯2h)(∂I∗ j− 12 )+ − (a1h[p¯h]2)(∂I∗ j− 12 )− + (a1h(p¯ + h ) 2)(∂I∗ j− 12 )− + 2 ∫ I∗ j− 12 σ1p¯ 2 hdx}. 因 (∂I∗ j− 1 2 )+ = (∂I ∗ j )−, (∂I ∗ j )+ = (∂I ∗ j+ 1 2 )−, A1(p¯h, p¯h) = 1 2 n∑ j=0 {(|a1h|[p¯h]2)(∂I∗j )− + (|a1h|[p¯h]2)(∂I∗j+12 )−}+ ∫ I∗j σ1p¯ 2 hdx+ ∫ I∗ j− 12 σ1p¯ 2 hdx + 1 2 (|a1h|(p¯+h )2)(∂I∗0 )+ − (|a1h|(p¯+h )2)(∂I∗n)− , – 10 – 第 1章 一维浅水方程的二阶广义迎风差分方法 而 p¯+h |(∂I∗n)−= 0. 故 A1(p¯h, p¯h) ≥ γ0(‖p¯h‖2∂I + ‖p¯h‖20). 证毕. 所以,对 p¯ ∈ H2(I), Ritz投影 Rhp¯ ∈ Vh存在且唯一,且有: A1(Rhp¯, vh) = A1(p¯, vh), ∀vh ∈ Vh, (1.3.1) 同理: A2(q¯h, q¯h) ≥ γ0(‖q¯h‖2∂I + ‖q¯h‖20). 其中: γ0 = min(σ1, 12), ‖q¯h‖20 = (q¯h, q¯h), ‖q¯h‖2∂I = n∑ j=0 {(|a2h|[q¯h]2)(∂I∗j )− + (|a2h|[q¯h]2)(∂I∗j+12 )−}+ 1 2 (|a2h|(q¯+h )2)(∂I∗0 )+ . 即 A2(q¯h, q¯h)是正定的. 所以,对 q¯ ∈ H2(I), Ritz投影 Rhq¯ ∈ Vh存在且唯一,且 有: A2(Rhq¯, vh) = A2(q¯, vh), ∀vh ∈ Vh, (1.3.2) 引引引理理理 1.2: 若 p¯, q¯ ∈ H2(I),则有下列估计： |‖p¯−Rhp¯|‖ ≤ Ch 32‖p¯‖2, |‖q¯ −Rhq¯|‖ ≤ Ch 32‖q¯‖2, (1.3.3) 其中 |‖v|‖2 = ‖v‖20 + ‖v‖2∂I + h ∑ xi ∫ I∗xi (a1 ∂v ∂xi )2dx,则 |‖v|‖和 ‖v‖在 Vh中等价. 对全离散广义迎风差分格式,下列误差估计成立: 定定定理理理 1.3: 设 p, q, ph, qh 分别是 (1.2.1), (1.2.6) 的解, 并且满足 pt, qt ∈ H2(I), ptt, qtt ∈ L2(I),则下列误差估计成立: ‖p(tn)− pnh‖0 ≤ ‖p0− p0h‖0 +Ch 3 2‖p0‖2 +4t ∫ tn 0 ‖ptt‖0dt+Ch 32 ∫ tn 0 ‖pt(t)‖2dt, (1.3.4) – 11 – 第 1章 一维浅水方程的二阶广义迎风差分方法 和 ‖q(tn)− qnh‖0 ≤ ‖q0 − q0h‖0 + Ch 3 2‖q0‖2 +4t ∫ tn 0 ‖qtt‖0dt+ Ch 32 ∫ tn 0 ‖qt(t)‖2dt. (1.3.5) 证证证明明明：：：首先证明 (1.3.4). 记 pnh − p(tn) = ξn1 + ηn1 ,其中 ξn1 = Rhp(tn)− p(tn), ηn1 = pnh −Rhp(tn). 由 (1.3.3)有 ‖ξn1 ‖ ≤ Ch 3 2‖p(tn)‖2. 又 p(tn) = p(0) + ∫ tn 0 pt(t)dt, ‖p(tn)‖2 ≤ ‖p(0)‖2 + ∫ tn 0 ‖pt(t)‖2dt. 则 ‖ξn1 ‖ ≤ Ch 3 2 (‖p0‖2 + ∫ tn 0 ‖pt(t)‖2dt). (1.3.6) 下面对 ηn1 进行估计.记 ∂tpnh = (pnh − pn−1h )/4t.由 (1.2.9)和 Rh的定义: A1(η n 1 , vh) = A1(p n h −Rhp(tn), vh) = −(∂tpnh, vh) + (sn1 , vh)− A1(p(tn), vh) = (pt(tn)− ∂tpnh, vh), 因此 (∂tη n 1 , η n 1 ) + A1(η n 1 , η n 1 ) = (pt(tn)−Rh∂tp(tn), ηn1 ) = (wn1 + w n 2 , η n 1 ), (1.3.7) – 12 – 第 1章 一维浅水方程的二阶广义迎风差分方法 其中 wn1 = pt(tn)− ∂tp(tn), wn2 = ∂tp(tn)−Rh∂tp(tn). 注意到 (ηn1 ) +|(∂I)− = 0, A1(ηn1 , ηn1 ) ≥ 0. 因此由 (1.3.7)我们有 ‖ηn1 ‖0 ≤ ‖ηn−11 ‖0 +4t‖wn1 + wn2‖0 ≤ ‖η01‖0 +4t n∑ j=1 ‖wj1 + wj2‖0. (1.3.8) 显然 ‖η01‖0 = ‖p0h −Rhp(x, 0)‖0 ≤ ‖p0h − p(x, 0)‖0 + ‖p(x, 0)−Rhp(x, 0)‖0 ≤ ‖p0h − p(x, 0)‖0 + Ch 3 2‖p(x, 0)‖2. (1.3.9) wj1 = pt(tj)−4t−1(p(tj)− p(tj−1)) = 4t−1 ∫ tj tj−1 (t− tj)ptt(t)dt, wj2 = (I −Rh)∂tp(tj) = 4t−1 ∫ tj tj−1 (I −Rh)pt(t)dt. 因此我们有: 4t n∑ j=1 ‖wj1 + wj2‖0 ≤ 4t ∫ tn 0 ‖ptt(t)‖0dt+ Ch 32 ∫ tn 0 ‖pt(t)‖2dt. (1.3.10) 把 (1.3.9), (1.3.10)代入 (1.3.8)得 ‖ηn1 ‖0 ≤ ‖p0 − p0h‖0 + Ch 3 2‖p0‖2 +4t ∫ tn 0 ‖ptt(t)‖0dt+ Ch 32 ∫ tn 0 ‖pt(t)‖2dt. (1.3.11) – 13 – 第 1章 一维浅水方程的二阶广义迎风差分方法 最后联合 (1.3.6), (1.3.11)得到误差估计 ‖pnh − p(tn)‖0 ≤ ‖ηn1 ‖0 + ‖ξn1 ‖0 ≤ ‖p0 − p0h‖0 + Ch 3 2‖p0‖2 +4t ∫ tn 0 ‖ptt(t)‖0dt +Ch 3 2 ∫ tn 0 ‖pt(t)‖2dt. 则 (1.3.4)得证. 同理 (1.3.5)成立. 证毕. 1.4 数数数值值值实实实验验验 算例一:静水流动,初始流速为零,水深为 1米,计算区域为 [−100, 100] [单 位: 米],空间步长为 h = 2米,时间步长为 ∆t = 0.2秒. 计算结果表明水深保持 不变.见图 1− 1. 算例二: 用经典算例理想溃坝水流来验证上面算法的有效性. 理想溃坝问 题是典型的非线性流的例子,经常用来检验数值格式,设右端项 s = 0,初始条 件为: H(x, 0) = { 1.0, x ≤ 0, 0.5, x > 0. u(x, 0) = 0. 下面用简单中心格式 [14] 和本文的方法做比较来说明格式的有效性. 在例 子中我们取空间区域 I = [−100, 100](单位: 米),空间步长为 h = 2米,时间步 长为 ∆t = 0.2秒. 图 1− 2为水深在时刻 T = 5秒时简单中心格式与广义差分 的比较,图 1− 3为水深在 T = 10秒时的变化图,图 1− 4为水深在 T = 20秒 时的变化图. 定义克朗数 Cr1 = a1∆t/h, Cr2 = a2∆t/h,方程 (1.1.2)的简单中 心格式为: pn+1i = p n i − Cr1 2 (pni+1 − pni−1) + Cr21 2 (pni+1 − 2pni + pni−1), qn+1i = q n i − Cr2 2 (qni+1 − qni−1) + Cr22 2 (qni+1 − 2qni + qni−1). – 14 – 第 1章 一维浅水方程的二阶广义迎风差分方法 −100 −80 −60 −40 −20 0 20 40 60 80 100 1 1 1 1 1 1 1 图 1-1静水流动 −100 −50 0 50 100 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 −100 −50 0 50 100 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 图 1-2 (左)简单中心格式, (右)GUDM,T=5s – 15 – 第 1章 一维浅水方程的二阶广义迎风差分方法 −100 −50 0 50 100 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 −100 −50 0 50 100 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 图 1-3 (左)简单中心格式, (右)GUDM,T=10s −100 −50 0 50 100 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 −100 −50 0 50 100 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 图 1-4 (左)简单中心格式, (右)GUDM,T=20s – 16 – 第 1章 一维浅水方程的二阶广义迎风差分方法 从图 1− 2, 1− 3, 1− 4可以看出在同一时刻,简单中心格式在间断处产生 了假震现象，而本文的广义迎风差分的结果消除了假震. 从图 1 − 3右与图 1− 4右可以看出随着时间的变化,水流左边的水深不断下降,这与实际相符. – 17 – 第 2章 二维浅水方程的二阶广义迎风差分方法 第第第 2章章章 二二二维维维浅浅浅水水水方方方程程程的的的二二二阶阶阶广广广义义义迎迎迎风风风差差差分分分方方方法法法 在本章里将研究二维浅水方程 [15],给出其广义迎风差分格式,并对所给的 格式进行 L2范数误差估计,最后用数值例子说明所给格式的有效性． 2.1 控控控制制制方方方程程程 考虑如下二维浅水方程: ∂H ∂t + ∂(Hu) ∂x + ∂(Hv) ∂y = 0, (2.1.1a) ∂(Hu) ∂t + ∂(Hu2 + gH 2 2 ) ∂x + ∂(Huv) ∂y = 0, (2.1.1b) ∂(Hv) ∂t + ∂(Huv) ∂x + ∂(Hv2 + gH 2 2 ) ∂y = 0. (2.1.1c) 其中 (x, y) ∈ Ω为具有光滑边界的有界区域, t ∈ [0, T ], H 为水深, u, v 分别为 x, y方向的平均水流速度, g为重力加速度. 将上述方程改写成下式: Ht + (u, v) · ∇H +Hux +Hvy = 0, ut + (u, v) · ∇u+ gHx = 0, vt + (u, v) · ∇v + gHy = 0. (2.1.2) 2.2 广广广义义义迎迎迎风风风差差差分分分格格格式式式 设 Ω ⊂ R2 是有界区域,并且 ∂Ω是其分片光滑边界, Th = {kQ}是区域 Ω 的拟一致三角剖分, Q为三角形重心,三角形的顶点和边的中点是剖分的结点. Ωh 为所有三角形单元顶点集, Mh 为所有单元边的中心集, Ω∗h 为所有单元重 心的集合. Ω˙h = Ω¯h\∂Ω, M˙h = M¯h\∂Ω. – 18 – 第 2章 二维浅水方程的二阶广义迎风差分方法 Th的对偶剖分记为 T ∗h ,包含围绕 P0 ∈ Ωh的区域K∗P0 和围绕M0 ∈Mh的 区域 K∗M0 . 这些小区域叫做对偶单元,他们的构造如下: 1) K∗P0 的构造.假设 P0 ∈ Ωh，Pi(i = 1, 2, · · · , 7)为 P0 的相邻顶点, P0i 为 P0Pi 上的点且 P0P0i = 13P0Pi. 顺次连接 P0i(i = 1, 2, · · · , 7)得到一个围绕 P0 的区域 K∗P0 . 见图 2− 1. P1 P2 P3 P4 P5 P6 P7 P0 P01 P02 P03 P04 P05 P06 P07 图 2-1 K∗P0 的构造 2)K∗M0 的构造.令 M0 ∈ Mh 为相邻三角单元 KQ1 = 4P0P1P2 和 KQ2 = 4P0P1P3 的公共边上的中点, 记 Q12, Q13, Q02, Q03 分别为 P01P02, P01P03, P10P12 和 P10P13 的中点,依次连接 P10, Q03, Q2, Q13, P01, Q12, Q1,Q02, P10,则 得到围绕M0的区域 K∗M0 . 见图 2− 2. 对 Th 选 Lagrange二次元空间. 对 ∀P0 ∈ Ω˙h,M0 ∈ M˙h,相应的基函数为 分片二次多项式满足下列插值条件: φP0(P ) = { 1, P = P0, 0, P ∈ Ωh ⋃ Mh \ {P0}. – 19 – 第 2章 二维浅水方程的二阶广义迎风差分方法 p0 p3 p2 p1 Q2 Q1 Q12 Q13 Q02 Q03 M0 p01 p10 p12 p13 p03 p02 图 2-2 K∗M0 的构造 φM0(P ) = { 1, P =M0, 0, P ∈ Ωh ⋃ Mh \ {M0}. 因此试探函数空间为 Uh = span{φP0 , φM0 , P0 ∈ Ω˙h,M0 ∈ M˙h}. 对 T ∗h 选分片常数空间, ∀P0 ∈ Ω˙h,M0 ∈ M˙h,相应的基函数为 K∗P0，K∗M0 的特征函数. ψP0(P ) = { 1, P ∈ K∗P0 , 0, P ∈¯K∗P0 . ψM0(P ) = { 1, P ∈ K∗M0 , 0, P ∈¯K∗M0 . 因此检验函数空间为 Vh = span{ψP0 , ψM0 , P0 ∈ Ω˙h,M0 ∈ M˙h}. 令 a = (u, v),将边界 ∂Ω分成两部分:{ (∂Ω)− = {(x, y) ∈ ∂Ω : a · γ ≤ 0}, (∂Ω)+ = {(x, y) ∈ ∂Ω : a · γ > 0}. (2.2.1) – 20 – 第 2章 二维浅水方程的二阶广义迎风差分方法 其中 γ 为 ∂Ω的外法向量. 由于空间在对偶单元边界的不连续性, 不能在整个区域 Ω 上应用 Galerkin方法,但可以在单个的对偶单元 K∗P0 和 K∗M0 上应用. 首先我们考虑 (2.1.2)中的第一式. 设 w1h ∈ Vh满足:∫ K∗P0 [Hh,t + ah · ∇Hh +Hhuh,x +Hhvh,y]w1hdxdy = 0, (2.2.2) ∫ K∗M0 [Hh,t + ah · ∇Hh +Hhuh,x +Hhvh,y]w1hdxdy = 0. (2.2.3) 记 γ, ν 分别为 ∂K∗P0 , ∂K∗M0 的外法向量,对 (2.2.2)式中的第二项应用 Green公 式得∫ K∗P0 (ah · ∇Hh)w1hdxdy = − ∫ K∗P0 Hh · div(ahw1h)dxdy + ∫ ∂K∗P0 (ah · γ)Hhw1hds. (2.2.4) 代入 (2.2.2)得:∫ K∗P0 Hh,tw1hdxdy − ∫ K∗P0 Hh · div(ahw1h)dxdy + ∫ ∂K∗P0 (ah · γ)Hhw1hds, + ∫ K∗P0 Hhuh,xw1hdxdy + ∫ K∗P0 Hhvh,yw1hdxdy = 0. (2.2.5) 类似 (2.2.1)定义 (∂K∗P0)+和 (∂K∗P0)−,对 ∀(x, y) ∈ ∂K∗P0 ,设 H+h (x, y) =  lim (x′,y′)−→(x,y) Hh(x ′, y′),当(x, y) ∈ (∂K∗P0)−, (x′, y′)∈¯K∗P0时, lim (x′,y′)−→(x,y) Hh(x ′, y′),当(x, y) ∈ (∂K∗P0)+, (x′, y′) ∈ K∗P0时. H−h (x, y) =  lim (x′,y′)−→(x,y) Hh(x ′, y′),当(x, y) ∈ (∂K∗P0)−, (x′, y′) ∈ K∗P0时, lim (x′,y′)−→(x,y) Hh(x ′, y′),当(x, y) ∈ (∂K∗P0)+, (x′, y′)∈¯K∗P0时. – 21 – 第 2章 二维浅水方程的二阶广义迎风差分方法 用 H+h 代替 (2.2.5)左边的线积分得:∫ K∗P0 Hh,tw1hdxdy − ∫ K∗P0 Hh · div(ahw1h)dxdy + ∫ ∂K∗P0 (ah · γ)H+h w1hds + ∫ K∗P0 Hhuh,xw1hdxdy + ∫ K∗P0 Hhvh,yw1hdxdy = 0. (2.2.6) 由 (2.2.4)有 − ∫ K∗P0 Hh · div(ahw1h)dxdy = ∫ K∗P0 (ah · ∇Hh)w1hdxdy − ∫ ∂K∗P0 (ah · γ)Hhw1hds = ∫ K∗P0 (ah · ∇Hh)w1hdxdy − ∫ (∂K∗P0)+ (ah · γ)Hhw1hds− ∫ (∂K∗P0)− (ah · γ)Hhw1hds. 则 (2.2.6)式写为∫ K∗P0 Hh,tw1hdxdy + ∫ K∗P0 (ah · ∇Hh)w1hdxdy + ∫ (∂K∗P0)− (ah · γ)[Hh]w1hds + ∫ K∗P0 Hhuh,xw1hdxdy + ∫ K∗P0 Hhvh,yw1hdxdy = 0. (2.2.7) 同理对 (2.2.2)式的第三第四项应用 Green公式,得到:∫ K∗P0 Hhuh,xw1hdxdy = ∫ K∗P0 Hhuh,xw1hdxdy + ∫ (∂K∗P0)− (Hhnx)[uh]w1hds, ∫ K∗P0 Hhvh,yw1hdxdy = ∫ K∗P0 Hhvh,yw1hdxdy + ∫ (∂K∗P0)− (Hhny)[vh]w1hds. – 22 – 第 2章 二维浅水方程的二阶广义迎风差分方法 把上式代入 (2.2.7)即得 (2.2.2)的半离散广义迎风差分格式:∫ K∗P0 Hh,tw1hdxdy + ∫ K∗P0 (ah · ∇Hh)w1hdxdy + ∫ (∂K∗P0)− (ah · γ)[Hh]w1hds + ∫ K∗P0 Hhuh,xw1hdxdy + ∫ (∂K∗P0)− (Hhnx)[uh]w1hds + ∫ K∗P0 Hhvh,yw1hdxdy + ∫ (∂K∗P0)− (Hhny)[vh]w1hds = 0. (2.2.8) 这里 [zh] = z+h − z−h 是指 zh通过 ∂(K∗P0)−的跳跃,其中 z = H, u, v. nx, ny 分别 为 x, y方向的方向导数. 同理对 (2.2.3)应用 Green公式得半离散广义迎风差分格式:∫ K∗M0 Hh,tw1hdxdy + ∫ K∗M0 (ah · ∇Hh)w1hdxdy + ∫ (∂K∗M0)− (ah · ν)[Hh]w1hds + ∫ K∗M0 Hhuh,xw1hdxdy + ∫ (∂K∗M0)− (Hhnx)[uh]w1hds + ∫ K∗M0 Hhvh,yw1hdxdy + ∫ (∂K∗M0)− (Hhny)[vh]w1hds = 0. (2.2.9) 同理, (2.1.2)的第二式和第三式的半离散广义迎风差分格式为:∫ K∗P0 uh,tw1hdxdy + ∫ K∗P0 (ah · ∇uh)w1hdxdy + ∫ (∂K∗P0)− (ah · γ)[uh]w1hds + ∫ K∗P0 gHh,xw1hdxdy + ∫ (∂K∗P0)− (gnx)[Hh]w1hds =