第 1章导论
1.1流体力学多介质大变形问题的研究背景
在内爆动力学、惯性约束聚变 (inertiAl con.ned fusion,ICF)、界面不稳定性、高速冲击、先进常规兵器等国防和高新技术领域 ,大规模、高置信度数值模拟是必不可少的工具 ,正在发挥越来越重要的作用 .数值模拟作为除了理论研究和实验研究之外的第三种研究方法 ,它既是连接理论研究和实验研究的桥梁 ,也是一种新的实验手段 (即数值试验 ).它对研究和发现物理过程中的规律 ,认识其中的机理具有重要的作用 .
上述这些领域所遇到的物理问题常常非常复杂 .从数学上来看 ,主要表现在需要求解可压缩流体力学方程组 ,其计算区域往往是三维复杂区域 ,问题具有动边界、多介质、大变形、强间断、强非线性以及多物理过程强耦合等特点 ,并且要求准确、清晰地刻画物质界面 .其中多介质大变形问题是流体力学计算的难点 ,迄今没有得到很好的解决 .
例如 ,为了探索实现氘氚核聚变 ,提供解决人类清洁能源问题的最终有效途径 ,中国和世界科技强国都开展了 ICF研究 .在激光间接驱动 ICF过程中 ,首先将多束高能激光注入黑腔 ,腔壁在吸收了入射的激光后将其转换成软 X射线辐射 (如图 1.1左图所示 ),之后再利用辐射驱动猛烈压缩置于腔体中间的靶丸 (氘氚燃料小球,外覆碳氢等多层材料 ),使其达到极高能量密度之聚变点火燃烧条件 .在靶丸压缩过程中 ,由于存在复杂波系的相互作用 ,其烧蚀层、氘氚冰层、氘氚气的各个界面上会出现流体力学界面不稳定性 (如图 1.1右图所示 ),包括 RAyleigh-TAylor不稳定性和 Richtmyer-Meshkov不稳定性 .当界面上的压强梯度和密度梯度方向相反时会出现 RAyleigh-TAylor不稳定性 ,当激波经过两种流体的界面时会出现 Richtmyer-Meshkov不稳定性 .界面不稳定性导致物质界面变形、翻转、破碎以及后期的湍流混合现象 ,直接影响靶丸的最终压缩状态 ,对聚变成功与否至为关键 .因此在 ICF中界面不稳定性以及与之相关的多介质大变形是最具挑战性的问题之一 .
1.2拉氏方法、欧拉方法和任意拉氏欧拉方法
基于网格类的流体力学计算方法主要分为采用物质网格的拉氏方法 (LAgrAng-iAn method)、采用固定空间网格的欧拉方法 (EuleriAn method)和采用运动网格的任意拉氏欧拉方法 (ArbitrAry LAgrAngiAn EuleriAn method,简称 ALE方法).
拉氏方法采用物质网格 ,即网格的运动速度取为物质的运动速度 ,从而避免了输运计算及由此造成的误差 .同时用网格边界清晰地刻画、描述物质界面 ,界面处理的精度较高 .其缺点是当流场中发生大变形时网格极度扭曲 ,往往导致计算终止.图 1.2是拉氏方法示意图 ,其中 tn和 tn+1分别表示第 n时间层和第 n +1时间层 .拉氏方法最早是 von NeumAnn和 Richtmyer在文献 [1]中针对一维情况提出来的 . Wilkins在文献 [2]中将该方法推广到二维情况 .这种方法采用人工粘性(Arti.ciAl viscosity)处理激波间断 .采用这种方法的早期的计算格式求解时 ,由于该方法求解内能方程 ,总能量是不守恒的 ,并且会产生虚假的数值振荡 .尽管有这些缺点 ,在过去的几十年中 ,该方法在多介质流动的数值模拟中被广泛应用 .为了解决上述这些问题 ,最近十几年来 ,针对该方法的缺点作了很多改进 . CArAmAnA和 ShAshkov在文献 [3]中提出采用子网格压力消除沙漏运动和虚假漩涡 .在文献 [4,5]中, ShAshkov等基于支撑算子的思想建立了保持总能量守恒性的相容拉氏动力学计算方法 .其主要特点是保证控制方程的主要物理性质和方程中各微分算子之间严格的数学关联在离散后得到保持 ,将控制方程的基本物理性质融入到数值计算中去 ,从而改善了数值模拟结果的可靠性和预测能力 .此外 ,人工粘性的计算方法也取得了明显的进步 [6,7].
欧拉方法采用固定空间网格进行计算 ,能够处理大变形问题 .其缺点是由于输运计算带来的误差 ,欧拉方法很难给出精确的物质界面 .图 1.3是欧拉方法示意图 .近几十年来欧拉方法取得了长足的进步 ,发展了基于 RiemAnn间断解的 Godunov方法 [8]和 Roe方法 [9]、积分平均型间断解方法 MUSCL[10]和 PPM[11],高分辨率方法 TVD[12]、ENO[13]、WENO[14]、紧致格式 [15],以及间断有限元方法 [16-19]等.其中 ENO方法采用逐次扩展的节点模板来构造高阶的逼近多项式 ,利用比较差商绝对值大小进行取舍的方法判定、选择扩展的节点模板 ,尽量避免在所选择的模板中包含间断 ,以此提高插值方法精度并实现高分辨率和无振荡的效果 .在 ENO方法的实施过程中 ,有许多中间计算结果被可惜的丢弃不用 .文献 [14]提出了 WENO方法 ,对这个不足作了弥补 ,因此得到了更为广泛的应用 .间断有限元方法是提高数值逼近精度的一个有效方法 .间断有限元方法的出现 ,最早可以追溯到 1973年 Reed和 Hill关于中子输运方程问题的论文 [20]. 20世纪 90年代以来 ,以 Cockburn和 Chi-WAng Shu为代表提出的 Runge-KuttA间断 GAlerkin有限元方法[16-19]在解决含有间断现象的问题中发挥着越来越大的作用 ,卓有成效地应用到了水动力学、气动力学和波传播等问题 .间断有限元方法的基本思想是用检验函数乘以原方程并在网格单元上积分 ,然后通过分步积分获得原问题的弱形式 .选取有限元逼近函数空间 (允许其在单元边界不连续 ),选择合适的数值流通量构造单元边界连接条件 ,形成可解的封闭线性代数方程组 .间断有限元方法具有很多优点：易于处理复杂边界和边值问题 ;具有灵活处理间断的能力 ,克服了一般有限元不适合间断问题的缺点 ;可以通过提高单元插值函数的次数来提高精度 ,不必像有限体积方法那样首先扩大节点模板 ,再做复杂的多项式重构 ;容易实施自适应策略 ;可以适当选取基函数 ,使得质量矩阵是分块对角的 ,容易求逆 ;为了求解给定单元的自由度,只需要相邻单元的自由度 ,易于实现并行计算 .用间断有限元等高阶格式计算含有强间断的问题会出现非物理振荡甚至非线性不稳定 (Gibbs现象 ). Qiu JiAnxiAn和 Shu ChiwAng提出了一个小模板的 HWENO(Hermite weighted essentiAlly non-oscillAtory)重构方法 [21,22]来抑制非物理振荡 ,并用作间断 GAlerkin有限元方法的限制器 .此外 , H. Luo等也提出了一个类似的方法 [23].
为了克服拉氏方法和欧拉方法的缺点并结合它们的优点 ,人们提出并发展了 ALE(ArbitrAry LAgrAngiAn-EuleriAn)方法 [24-29],其主要特点是采用网格速度可人为控制的运动网格 ,物质内部网格可以通过优化保持较好的几何品质 ,能够适应大变形问题的数值模拟 ,同时物质界面仍用网格边界显式刻画 ,界面处理清晰、准确 .图 1.4是 ALE方法示意图 .由于上述优点 ,近几十年来 ALE方法一直是流体力学计算方法研究的热点 ,具有重要的学术价值和工程应用价值 ,在相关基础研究领域和应用基础研究领域 ,特别是内爆动力学和 ICF等领域得到广泛应用 .
根据物理量定义位置的不同 , ALE方法大致可以分为两类：一类为运动量 (如位移、速度和加速度 )定义在网格节点、状态量 (如密度、压力、内能等 )定义在网格中心的交错网格 ALE方法 (stAggered grid ALE method),一类为所有物理量 (包括运动量和状态量 )均定义在网格中心的非交错网格 ALE方法 (cell-centered ALE method).
另一方面 ,根据物理量更新方式的不同 , ALE方法大致可以分为两种：比较常见的一种 ALE方法通常由三个步骤组成 ,即拉氏步计算、网格重分和物理量重映 (如文献 [24-27]).当网格变形不大、形状较好时 ,仅进行拉氏步计算 ,把 tn时刻的拉氏网格上的物理量更新为 tn+1时刻的拉氏网格上的物理量 ,高精度地追踪物质界面 ;当 tn+1时刻的拉氏网格变形较大、网格品质较差时 ,首先对该网格进行优化 ,生成品质较好的新网格 ,然后进行重映计算 ,把 tn+1时刻的拉氏网格上的物理量映射到 tn+1时刻的新网格上 .另一种 ALE方法首先确定 tn时刻的网格速度 ,得到 tn+1时刻的新网格 ,然后直接求解含有网格速度的 ALE形式的流体力学方程组 (如 [28,29]),把 tn时刻的拉氏网格上的物理量直接更新为 tn+1时刻的新网格上的物理量 .文献 [30, 31]中指出 ,这两种 ALE方法的精度是相当的 .
1.3运动界面追踪方法
前面已经提到 ,由于输运计算带来的误差 ,欧拉方法很难给出精确的物质界面 .例如 , TVD、Roe等方法在光滑区域可以达到三阶、四阶精度 ,而在关键的间断区域却只能是一阶 (物质界面也是一种间断 ).为此 ,近几十年来人们发展了一类界面追踪方法 ,用于模拟、描述和追踪自由面和运动界面的轨迹、特征和发展 .目前比较流行的是 VOF(volume of .uid)方法 [32-33]、等值面 (level set)函数方法 [34,35]等.其中 VOF方法是在整个流场中定义一个函数 (称之为流体体积函数 ),在每个网格中 ,这个函数定义为一种流体 (称之为目标流体 )的体积与网格体积的比值 .在任意时刻 ,知道了这个流体体积函数在每个网格上的值 ,就可以通过某种途径显式地构造出运动界面 .然后在求解物理方程时可以在界面附近作特殊的精细处理 ,以提高分辨率和精度 .另一种较理想的做法是利用所谓的等值面函数 (level set function). (x,t)代替 VOF方法中的流体体积函数 (x是空间变量 , t表示时间 ).让 .以适当的速度移动 ,使其零等值面就是物质界面 .在任意时刻 ,只要知道 .,然后求出其零等值面,就知道了此时的运动界面 .等值面函数法不需要显式地追踪运动界面 ,从而可以较容易地处理复杂的物质界面及其拓扑结构发生变化的情形 ;而且界面的一些特征 (如法向、曲率等 )直接隐含在 level set函数中 ,便于精细地描述界面 ;此外 ,它还易于向高维推广 .由于这些优点 ,等值面函数法已经被用来处理几何、流体力学、工艺过程等许多方面的问题 .
1.4 MMALE方法
传统的 ALE方法在一个网格中只允许含有一种介质 ,在计算多介质问题时 ,一般将物质界面取为计算子区域的边界 ,与边界上的网格边重合 ,界面在法向保持拉氏运动 .这种情况下界面至少应该保持连续 ,最好不要发生较大的变形 ,只有这样才能保证在进行网格重分时能够在子区域中生成质量较好的新网格 .但是在内爆压缩等一些实际问题的数值模拟中 ,不仅网格随流体运动发生大变形 ,而且界面也会发生较大的变形 ,甚至出现破裂、合并等拓扑结构改变的情形 ,这时传统的拉氏方法和 ALE方法将无法继续应用 .美国 SAndiA国家实验室的 Peery等在 2000年提出的 MMALE方法 (multi-mAteriAl ArbitrAry LAgrAngiAn EuleriAn method)[36]是解决这个问题的一个比较有效的方法 ,其基本思想是把传统的 ALE方法和界面追踪方法相结合 ,处理包含网格扭曲以及界面破裂、合并等复杂现象的多介质大变形问题.其实现途径是在传统的 ALE方法的基础上 ,通过引进混合网格 (即允许一个网格内含有多种介质 ,物质界面可以穿过网格 ),用界面追踪方法在混合网格中刻画、重构运动界面 ,来进一步提高 ALE方法处理多介质大变形问题的能力 .他们将该方法应用于 ALEGRA程序中 ,计算 ICF研究中遇到的含有强剪切变形的多介质大变形问题 .图 1.5是 MMALE方法示意图 .
图 1.5 MMALE方法示意图
粗实线表示用网格边刻画的物质界面 ,虚线表示在混合网格中用界面追踪方法刻画、重构物质界面
文献 [36]中的 MMALE方法包含以下三个关键技术 .一是对混合网格提出一个热力学封闭模型 ,这个模型用来计算经过一个拉氏步以后混合网格中每种物质的体积分数和物理量的变化 .文献 [36]采用平均应变率混合模型 [37].在这个模型中 ,混合网格中各种介质的应变率就取为网格的平均应变率 ,这实际上就是假设经过一个拉氏步以后各种介质的体积份额保持不变 .当混合网格中界面两侧密度比或压力比很大时 ,应用这个模型会得到非物理的结果 [38].二是用重分重映处理网格的大变形.文献 [36]采用对流重映方法 [39-41],这种方法通过计算网格边上的通量来得到新网格上的积分量 ,因而要求新网格充分接近老网格 .但是这种方法不需要计算老网格和新网格的交点 ,计算效率较高 .三是用 VOF方法处理界面的大变形 .
近年来 ,由于内爆动力学和 ICF等领域中多介质大变形问题研究的需求牵引 , MMALE方法受到学术界的广泛关注 ,成为计算流体力学的一个研究热点 .
在混合网格的封闭模型方面 ,为了克服平均应变率混合模型的缺点 ,提出了几个新的封闭模型 ,例如： 1压力平衡模型 [42].这种模型假设混合网格中各种介质的压力在瞬间达到平衡 .压力平衡是非线性问题 ,需要设计复杂的迭代格式进行计算 .
2压力松弛模型 [43-46].这种模型引进一种类似于粘性的松弛机制使得混合网格中的压力趋于平衡 . 3一维子网格力学模型 [45,47-49].这个模型用声波 RiemAnn解估计界面速度 ,并由此计算每种物质的体积分数和物理量的改变量 . 4质量分数模型[50,51].这个模型也被称为浓度法 .就计算效率和有效性而言 ,压力松弛模型由于
展开
