移动网格上的二阶godunov型格式在一维爆轰波模拟中的应用.pdf

移动网格上的二阶G o d u n o v 型格式在一维爆轰 波模拟中的应用· 陈国贤，汤华中，张平文 ( 北京大学数学科学学院北京) 誓要= 本文考虑一维爆轰波的模拟我们的算法分为两部分，网格的移动和物理量的时问发展网格 的移动是基于变分的迭代网格分布物理量的时间发展是通过求解在逻辑区域上的控制方程使用的格 式是对原项没有分裂的基于M U S C L 型初始数据重构的二阶G o d u n o v 型格式数值结果显示我们设计的 算法的有效性 关键词：自适应移动网格方法爆轰波G o d u n o v 格式初始数据重构 1 引言 燃烧是很复杂的物理过程，涉及到如流体流动和化学反应等在通常情况下，火焰以几米每秒的 低速传播但是在特定情况下，传播速度可以达到很快( 如每秒2 0 0 0 米，甚至更快) ，这就是爆轰波 在爆轰过程中火焰具有一定的厚度，被称为反应层，反应层的厚度很小，与宏观流体的特征长度相差 好几个数量级正是由于这种巨大差别，燃烧的数值模拟非常困难 对爆轰波进行数值模拟的研究始于1 9 4 0 年代1 6 J 对此问题的数值方法大致分为两类：一类是低 分辨捕捉化学反应区方法，例如投影法1 2 , 1 2 ] 能在较粗的网格上正确捕捉间断的位置；另一类是高分辨 捕捉化学反应区方法，即在反应区布局足够的网格来达到对反应区高分辨捕捉的目的【4 ’州。

这两类方 法都需要在反应区进行特殊的计算处理 自适应移动网格方法是一种新的技术，它在科学与工程计算( 如在产生奇性的固体和流体计算) 中发挥着越来越重要的作用当前，已有一些成熟的实现网格自适应移动的方法，如基于变分原理的 移动网格方法[ 4 1 ，移动网格P D E s 方法【4 】，和基于D v i n s k y [ 5 1 的调和映射的移动网格方、法t 8 , 9 1 等A z a r e n o k 等⋯成功的把移动网格方法应用到了爆轰方程的计算，使用的格式是在移动网格上用G o d o n u v 型格式 解流体方程 本文对一维的爆轰方程进行模拟，采用的也是在移动网格上的G o d o n u v 型格式与【1 】不同的是， 第一，网格的移动使用的是【9 】中的基于变分的移动网格方法；第二，网格变动之后，我们将控制方程 映射到逻辑区域上，然后再用二阶G o d u n o v 格式进行求解 本文是这样组织的：第二节给出爆轰波的控制方程；第三节给出移动网格技术；第四节推导控制 方程在逻辑区域上的形式并给出它的离散方法；最后在第五节中我们给出一个数值实验证明格式的有 效性 2 控镧方程 一维理想爆轰波计算可以用如下的非线性偏微分方程组来描述 罢+ 辈爷( u ) ， ( 2 ．1 ) 研氖 、7 ‘’ 式中U = ( p ，p u ，E ，p Z ) r ，F ( U ) = ( ∥，倒2 + p ，甜( E + p ) ，m z ) r ，S ( U ) = ( 0 ，0 ，- p K ( T ) Z ) r ’部分得到了9 7 3 基金( 2 0 0 5 C B 3 2 1 7 0 3 ) ，N S F C 基金( 1 0 4 3 1 0 5 0 ，1 0 5 7 6 0 0 1 ) ，和计算物理实验室的资助． P 表示密度，U 表示气体的流动速度，E 表示气体的全能量，包括内能，动能和分子结构能三项，p 表示压强，丁表示温度，Z ∈【0 ，1 】表示未燃气体的质量百分比。

总能量E = p [ e + u 21 2 ] + p Z q o ，这 里，e 表示内能质量比，常数q 表示束缚能气体状态方程设为p 2 p e ( r 一1 ) ，其中y 表示绝热指数 我们这里考虑的是简单的模型，假定7 在已燃和未燃区域为相同的常数温度由丁7 ／p R 给出，其 中R 表示摩尔气体常数作为温度T 的函数K 仃) 0 表示化学反应速度，不同的模型有不同的表达 式，常见的有点火模型和A r r h e n i u s 模型在本文中的数值实验中，我们用的是点火模型 砸，= ∽麓黧P ， 亿2 ， 其中％称为点火温度，％是化学反应的时间尺度这里还有一个很重要的量％澍，它表示 激波速度，我们在移动网格时要用到这个量，其详细的定义可以参考【6 】 3 白适应移动网格方法 设逻辑区域Q ，= 【0 ，1 ] 上的坐标用孝表示，物理区域Q [ 而，■】上的坐标用x 表示这一节我们 沿用【9 】中设计的网格移动策略，但不对物理量进行守恒型插值 我们建立的从逻辑区域Q ，到物理区域Q 的一一映射 x = x ( f ) ， f ∈Q ，， ( 3 ．1 ) 是基于如下泛函的极小化 阶去t 文警砖 @ 2 ， 式中，国 0 称为控制函数，它一般情况下是非线性的，它的取值反应网格的集中。

利用W i n s l o w [ 1 1 1 的方法，泛函( 3 ．2 ) 的极小给出如下E u l e r - L a g r a n g e 方程 ( 嬲善) f = 0 ， ( 3 ．3 ) 加上边界条件x ( O ) = X 1 ，x ( 1 ) = x ，在一维情况，( 3 ．3 ) 等价于等分布原理我们给出逻辑区域的 均匀剖分{ 参I 磊= O 一1 ) ／Ⅳ，i = 1 ，2 ，．．．，N + I ) ，记△孝= 1 ／Ⅳ表示逻辑区域上的网格步长那么，( 3 ．3 ) 可以用一些迭代方法求解，比如J a c o b i 迭代， c o ( u ( x , + l ，2 ) ) ( 薯+ l 一置) ～缈( 甜( ■一l ，2 ) ) ( 霉一薯一1 ) = 0 ， ( 3 ．4 ) 其中‘+ l ，2 = ( ‘+ l + t ) ／2 网格边界条件由X 1 = x f ，x Ⅳ+ l = x ，给出从而解( 3 ．4 ) 可得到物理区域 7 4 上的非均匀剖分{ 薯Ix ，= x l z 2 ⋯一 h h “= x ，，i = 1 ，2 ，．．．，N + 1 ) ， 这里我们记 ( A x ) f + l ，2 = x f + l —t 表示物理区域上的网格步长。

一般情况下，物理区域Q 是固定的，用( 3 ．4 ) 来实现网格的移动就足够了但是，我们这里要模 拟的是爆轰波，它是～种超音速波，于是要模拟它，需要给出足够大的物理区域，使得在整个计算过 程中，我们关心的流场都在区域中，这会给计算量带来很大的浪费，所以我们这里让我们计算的物理 区域Q 也是随着流体运动而变化，这样可以大大地节省计算量为此，我们定义 啤( f ) = 【西( f ) ，‘( f ) ] ，啤( f ) 的移动速度为％离散形式为 霉= t + 出％ ( 3 ．5 ) 其中△f 是在计算流体方程时用的时间步长，墨表示新的网格位置，芦l 或，． 3 ．5 ) 对于很多计算 都是有效的，特别是考虑高速流，如超音速流，或者关心流体的长时间性态的模拟，只需要把上b 改 成相应的宏观速度从而我们的网格移动技术由两部分组成，先解( 3 ．5 ) ，再解( 3 ．4 ) 4 流体在逻辑区域上的表示以及相应的离散 在上一节中，我们给出了网格的移动技术但是我们并没有同时进行物理量在新的网格上的投影 | 8 , 9 1 本节，我们将流体方程映射到逻辑区域上，然后在逻辑区域上求解，从而回避新旧网格之间解 的插值。

每一个时间层上，我们已经给出逻辑区域Q ，到物理区域Q 的一一映射( 3 ．1 ) 显然在不同的时间， 这个映射是不一定相同的于是在整个的时间空间来看，Q ，到Q 之间的映射可以表示成 { x 翟力’ @ ·， 其中f 表示的是当空间取为逻辑区域时的时间变量在不引起混淆的情况下，变量在逻辑区域上 所采用的记号取成和物理区域上的记号相同，比如U ( 孝，f ) = ：U @ ( 孝，f ) ，f ) 对于任何变量g ，有 鱼：亟+ 一O q —O x ，鱼：一O q —O x ， ( 4 ．2 ) 8 t戳a x8 f j8 芒0 x8 芒j 、1 从而 型：型+ 型曼：一—O F ( —U ) + s ( u ) + 型戈 a f斑苏舐 、7 舐 ：一．O F ( U ) ／O x + s ( u ) + 型戈／鱼， = 一· 十·、¨，I 十一J ，一． a 乏a 芒、8 芒a 芒： 7 5 其中戈= 缸／a f ，两边同时乘以反／a 善得 记F ( U ) = ：F ( U ) 一Y o U ，则有 型鱼：一—a F ( —U ) + s ( ∽鱼+ 一a u 戈． 8 { 8 芒8 乏、。

8 芒8 芒 昙( u 妻) + 昙户( u ) = s ( ∽妻． ( 4 - 3 ) 8 T8 芒8 芒 、8 芒 、 方程( 4 ．3 ) 便是在逻辑区域上的流体方程下面对这个方程进行离散 设‰∽= 专p 似鳓帜= 击e 坝纠挚= 忐e 坝列进，同时我 们记霹l ，2 = 谚+ l ，2 ( 毛) 在区域噱，专+ 1 】×[ L ，L + l 】上对( 4 ．3 ) 积分z ⋯1 0 ， ( U A x ／A ：，3 + 1 ，2 一( u 血／△孝) ：1 ／2 一 △r 一亏∥2 一霉肘“2 △f+ @ ( u ) 簧／i + 1 ／2 1 ， 其中亏掣7 2 表示户( u ) 在r ∈【％，‘+ 上的积分平均值的一种近似消去△善，将△r 改记作肛， 等式两边同时除以缸茄j ：，则上式可以化简为 u j + U i n + 1 ／2 l - - 霹1 ，2 嘁1 ，2 ，一A v ㈨n + l 一一 苜露n+l／2“n+l／2怕(％_彳’磁Ax!'+112， ( 4 ．4 ) 其中童在网格边界由膏= @ 一x ) ／出给出为了给出一个二阶离散，我们需要( 一．7 j + l i n + l ，／2 2 ⋯霉7 7 2 的 一个二阶近似。

这里我们采用预估校正方法 第一步是预估步首先我们利用v a nL e e r t 加1 的重构技术在f 时刻的逻辑区域上重构分片线性函数 得到单元边界左右值W ’L 和W 一那么霉”1 7 2 ≈F ( y ( o ；W ⋯，W ’R ) ) 作为通量在预估步的近似，其 中y ( f ／f ；W 一，W ’R ) 是如下黎曼问题的解 f 昙cu 责，+ 壶户c u ) _ 0 ， a f a 毒a f 、 1u c 孝，f 【叼U 7 一L ，喜三： ‘4 ’5 ’ 7 6 x 0 ( 4 ．6 ) 除了有一个角度( 该角度的斜率为圣) 变化外，是仿射等价的设( 4 ．6 ) 的解为矽@ ／f ；u y ，W ’只) ， 则沙( o ；W ⋯，W ’尺) 等于≯( 戈；叫n ’，W ’R ) 霹学近似取作吸】，2 E h ( 4 ．4 ) 可以得到既j ：的预估值， 记作既；： 第二步是校正步利用预估值既；：，刚“ J ，7 i 川+ 1 ／门2 = ：( 既；：3 t ．．i 7 n + ，l ：) ／2 是二阶近似再一次利用 v a i l L e e r [ 1 0 1 的重构技术在t n 州2 逻辑区域上重构分片线性函数得到单元边界左右值W + 1 他’L 和 W “他一。

除了把W ’￡和W ’R 换成了W + 1 佗’上和W “他一，剩下的工作和预估步是完全一样，得到得 解嘴；2 就是一个二阶精度的近似 时间步长的取法仍然取成满足C F L 稳定性要求， A t = 印m i n 缸+ 1 ，2 ， 这里O c f l l ，每个网格上的时间步长△t + l ，2 由下式决定， △‰2 盂面瓦j ＆五( 7 i + 面1 ／2 瓦面， 其中砭1 ，2 和疋l ／2 表示U + l ／2 的左右特征值，毫+ l ／2 = @ j + 戈m ) ／2 5 数值算例 我们这里选取文献【1 ，7 ] h b 考。