偏微分方程的移动网格方法

第1章自适应方法 18世纪的拉普拉斯是一位法国的机械决定论者，被称为法国的“牛顿”，他把牛顿的质点运动确定论扩展到了无穷质点系统的确定论.拉普拉斯在《概率论的哲学试验》著作中写道：“我们可以把宇宙现在的状态看作是它历史的果和未来的因.如果存在这么一个智者，它在某一时刻，能够获知驱动这个自然运动的所有的力，以及组成这个世界的所有物体的位置，并且这个智者有足够强大的能力，可以把这些数据进行分析，那么宇宙之中从*宏大的天体到*渺小的原子都将包含在一个运动方程之中；对这个智者而言，未来将无一不确定，恰如历史一样，在它眼前一览无遗.” 拉普拉斯的这段名言，在科学和哲学界引起了轩然大波，余波至今未消.拉普拉斯这里所说的“智者”便是后人所称的“拉普拉斯妖”（Laplace’s　demon).事实上，拉普拉斯希望找到一个独立的公式，把宇宙的万物运动描述清楚.他提到：公式中要包含力、位置和原子状态等的描述.这样，宇宙的前因后果都确定了，也都能回溯过去和预测未来了.直到现在，人们还在不断完善公式，发展探测工具，获取高分辨率的资料，努力实现拉普拉斯的理想目标. 实际上，在自然科学与工程技术中，很多运动发展过程与平衡现象会遵循一定的规律.这些规律的定量表述一般地呈现为含有未知函数及其导数的方程.我们将只含有未知多元函数及其偏导数的方程，称为偏微分方程，初始条件和边界条件称为定解条件;而偏微分方程和定解条件作为一个整体，称为定解问题.我们居住的地球，表面上空被一厚度十几公里到二十多公里的大气层所环绕.我们每天感受到的阴晴雨霜、冷暖风雪天气就发生在这十几到二十几公里厚的大气层里.大气环绕着地球每天都在运动变化，它遵循牛顿运动定理、质量守恒定理、大气状态方程、热力学定理和水汽守恒定理.数值天气预报，就是将描述大气的流体动力学纳维-斯托克斯（Navier-Stokes)方程、热力学方程组，根据某一时刻观测到的大气状态，用数学方法求解，得到未来某一个时间的大气状态. 欧拉方程组可以理解为纳维-斯托克斯方程的简化形式.欧拉方程组适用的地方很多，可以描述飞行中的流场变化，比如飞行器在空中速度过快丨马赫数大于1)的情况下，会使周围的压力、空气密度、温度等产生不连续分布，即产生了激波.这些都可以通过欧拉方程组描述出来. 广义相对论是爱因斯坦在1915年提出来描述引力现象的几何理论，其基本观点是时空结构取决于物质的运动及分布.爱因斯坦提出的引力场方程，体现了运动的物质及其分布决定周围的时空性质，对于任意坐标变换，场方程的形式不变.求解爱因斯坦方程是人们了解宇宙运行规律的前提条件，但是爱因斯坦场方程是一个强非线性偏微分方程组，是自然科学中*复杂的偏微分方程之一，因此想要求得其精确解十分困难.尽管如此，仍有相当数量的精确解被求得，但仅有少数具有物理上的直接应用. 1948年，柯朗（Courant)和弗里德里希斯（Friedrichs)合作出版了《超音速流和激波》(Supersonic　Flow　and　Shock　Waves)，这是一本经典性的理论著作.这本书问世之后，刻画激波的守恒律方程为20世纪50年代兴起的一个主要研究领域.此类型方程的特征之一就是：即使初始数据是充分光滑的，守恒律的解在有限时刻也可能会发生间断，形成激波、切向间断和稀疏波.由这一理论形成的空气动力学偏微分方程组，成了研究高速飞行器、核武器等的重要工具. 在上述例子中，要想找到描述大气的流体动力学纳维-斯托克斯方程的解，找到描述激波的欧拉方程组的解，找到描述引力现象的爱因斯坦方程的解，或者是描述激波的双曲型方程的解，大多数情况下是不可能的.换句话说，要把这众多的数学方程（组）求解出来，给出一个数学公式，是十分困难或几乎不可能的事情.在这种情况下，很长时间以来，人们只能通过很多简化，找到有限个解析解，这对复杂问题的理解有很大的局限性. 自从20世纪50年代电子计算机逐渐出现在科学计算中以来，以上问题得到了根本解决.人们开始结合数学和计算机，通过数值计算方法得到满足精度要求的微分方程近似解.这样，纳维-斯托克斯方程、欧拉方程组、爱因斯坦方程等，都可以得到令人满意的近似解. 数值计算的基本思路就是用简单问题近似复杂问题（两种问题具有相同或非常接近的解)，用有限空间代替无限维空间，用有限过程代替无限过程，用代数方程代替微分方程，用线性问题代替非线性问题.而微分方程的数值方法，无论是常微分方程还是偏微分方程，都是将连续的、无限未知数的问题近似为离散的、有限未知数的问题，并进一步数值求解.经典数值分析通常会关心如下一些问题：相容性、稳定性、收敛性、收敛阶、计算量等.相容性是说格式在局部是不是做出了正确的近似，有一定的“精度”；稳定性是说局部的近似误差会不会随着计算而积累放大；收敛性是说当离散尺度无穷小的时候数值解是否会趋向于真实解；收敛阶则刻画了收敛的速度，高阶的格式可以用较大的离散尺度获得较好的数值结果.因此，数值方法的*终表现需要在精度、稳定性和计算量之间找到一个平衡. 偏微分方程的数值方法常见的有有限差分方法、有限元方法、有限体积法和谱方法等.本书根据不同的问题，分别采用前三种方法，即有限差分方法、有限元方法和有限体积法. **类方法是有限差分法，主要推导工具是泰勒展开，它是*早用来求解偏微分方程定解问题的数值方法，也是应用*广泛的方法之一.其基本思想是：**步，对求解区域作网格剖分（二维一般是正方形网格或长方形网格；三维就是立方体或长方体网格)，使得自变量的连续变化区域被有限离散点（网格点）集代替；第二步，将问题中出现的连续变量的函数用定义在网格点上的离散变量代替，通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组.如果差分格式有解，且当网格无限变小时其解收敛于微分方程定解问题的解，则差分格式的解就作为原问题的近似解丨数值解).这个方法主要基于泰勒展开，其优势是简单易行，容易理解、容易实现，但局限性就是对求解区域的要求相对苛刻，复杂区域的突现变得比较烦琐. 第二类方法是有限元方法，主要基于变分原理.有限元方法的**步是对整个求解区域进行分解，使每个子区域都成为简单的部分（比如在平面区域可以是一个个小三角形，这些小三角形内部互不相交，但所有小三角形又充满了整个给定的区域).这种简单剖分被称作有限元，而它形成的数值方法则被称作有限元方法.有限元方法的第二步就是对求解的偏微分方程采用“广义函数”，把偏微分方程转换为在“更弱”的函数空间（通常是分片多项式空间）上成立的积分形式.很多实际问题可以化为数学上的“泛函”，然后要找到*小化泛函，而达到这一目的就需要采用变分方法.有限元的核心思想就是假定未知函数在“更弱”的函数空间具有简单的表达式，比如在每个单元上都是一阶多项式，这样就会由变分形式推出每个单元上简单的代数方程组，而代数方程组问题可以交由计算机求解.有限元方法可以达到很高的计算精度，并且可以应对复杂求解区域，因此在科学和工程界非常受欢迎.另外，有限元方法的求解步骤可以系统化、标准化，能够开发出灵活通用的计算机程序，广泛应用于很多实际问题. 著名数学家冯康曾用简单形象的比喻形容有限元方法：分整为零、裁弯取直、以简驭繁，化难为易.他还形象地总结了有限元方法的巨大作用：求解微分方程的定解问题好像是大海搜针，成功的可能是微乎其微；但有限元离散后，寻求近似解就好像是碗里捞针，显而易见容易多了. 举个例子说明有限元计算的重要性.1991年8月23日，在挪威北海“SleipnerA”石油钻井平台的*后建造期间，发生了一次灾难性的故障.原始船体:塌，造成7亿美元的损失和里氏3.0级地震.这个钻井平台设计高度是82米，有24个格室，底座建筑面积有16000平方米.斯堪的纳维亚独立研究机构SINTEF主持的调查结论显示，基础结构24个格室中的一个格室壁破裂，导致泄漏量超出泵机的处理能力.根据SINTEF的结论，基于线性弹性模型的有限元计算不够准确，导致剪切应力被低估47%.事故发生之后，更精确的有限元计算结果显示原始设计将会在62米深度发生故障，与实际发生故障的65米深度基本匹配.也就是说，粗糖的有限元计算是这次灾难的罪魁祸首，如有相关的有限元计算结果足够准确，这个灾难将会避免. 第三类方法是有限体积法，也叫控制体积法.有限体积法是在有限差分基础上发展起来的，同时它又吸收了有限元方法的一些优点.有限体积法易于人们理解和使用，并且可以得到更合理的物理解释.它*大的意义在于，使用有限体积法得到的离散方程，完美地体现了守恒性.其具体的步骤如下：**，在计算过程当中，将需要计算的区域分割成一连串的具有不重复的控制体积，使得各个得以控制的体积都可以有一个作为代表的节点；第二，通过对未知函数做出时空分片函数假设，在任意的控制体积内对微分方程作积分，并对积分量作合理近似，从而得出离散方程组，其中的未知量是网格点上的因变量.有限体积法获得的离散方程，物理上表示的是控制容积的通量平衡，方程中各项有明确的物理意义，这也是有限体积法与有限差分法和有限元方法相比更具优势的地方.有限体积法的解可以达到高精度，且能很好地保持守恒性，因此在很多物理和工程问题计算上非常受欢迎. 1.1自适应方法综述 在过去的数十年里，尽管计算机的速度和内存都有大幅度提高，但很多实际问题的数值模拟仍需在很多简化下才可完成.也就是说，现有的计算机能力对付过多的自由度仍然有很多困难，特别是三维空间问题.在这种情形下，自适应算法应运而生，并在实际应用和理论研究上受到了广泛的重视.如果一个偏微分方程的解有足够的光滑度，则一致网格（有时也叫均匀网格）就可以给出满意的解.但是，也有一些很重要的问题，其解的光滑性并非很好，比如说解间断或解有大梯度情形.局部的奇性会导致求解区域上的网格过细，会造成不必要的计算时间和数据储存上的浪费. 在数值计算中，自适应网格所起的作用就是：在物理解变化剧烈的区域，通过某些数学方法，让网格在迭代过程中不断调节，使得网格点分布与物理解特性相耦合，从而提高解的精度和分辨率.自适应网格希望在物理解变动大的区域自动聚集网格，而在物理解变化平缓区域稀疏网格，这样可以兼顾计算效率和解的精度. 网格自适应方法主要分为三种类型，分别叫做方法、方法和方法(图1.1).其中方法是对网格进行自适应的局部加密和稀疏化，这里h代表有限元离散的网格特征长度.也就是说，/i-方法表示网格的大小可以改变丨通过增加或减少网格数量)，但计算方法的格式精度（比如说线性元）保持不变.而^方法是在网格的不同位置根据解的光滑性质采用不同的基函数，p是多项式(polynomial)的缩写，浐方法表示局部地方不改变网格大小和位置，但通过提高局部的 1.1自适应方法综述 格式精度，比如说从一次元到二次元或高次元，来提髙数值解的分辨率.方法就是把/I-方法和方法有机地结合起来.简单地说，在解的光滑度比较髙的地方，可用高阶元来逼近，而在光滑度比较低的地方，则通过减少网格大小来提高解的精度.r-方法是进行网格点的重新分布，又叫做移动网格方法，这里的r是重新分布（redistribution)的意思.这个方法的特点是网格点数固定，相邻点的排列顺序不变，但在不同的时间层或迭代层上，格点的位置要根据某些准则来重新分布. 网格自适应方法*根本的目标在于使用*少的计算资源来解决复杂问题，从而可以在现有的硬件资源条件下扩大计算的规模和提高计算的精度. 进一步地，/I-方法的基本想法是在原有网格的基础上，通过近似解的某种后验估计和解在局部的误差表现，在有需要的地方做局部加密和稀疏.这个方法可以对解的误差加以控制，也就是说，当计算停止时，近似解可以达到所需要的精度.为了能够使得局部加密方法得到比较好的效果，需要两个方面的准备工作：一个是非常精细