（流体力学专业论文）三维非理想爆轰波阵面传播过程的数值方法和数值模拟.pdf

三维非理想爆轰波阵面传播过程的数值方法和数值模拟 摘要 本文基于爆轰冲击波动力学( d e t o n a t i o ns h o c kd y n a m i c s ) 方法，研究模拟 三维非理想爆轰波阵面传播过程的l e v e ls e t 方法。为了简化直角坐标系下弯曲 边界条件处理繁琐的缺点，采用了贴体坐标系，同时给出边界条件的处理方 法。根据h a m i l t o n j a c o b i 方程的g o d u n o v 差分格式，本文研究了非正交的贴体 坐标系中非理想爆轰波传播方程的差分格式并提出了相应的数值方法。我们编 写了适用于二维和三维非理想爆轰波传播计算的程序，给出了几个测试算例。 计算结果表明，在贴体坐标系中，本文的计算结果与解析计算结果符合，对于 三维旋转体模型，本文的二维和三维计算结果一致，验证了差分格式，计算方 法和计算程序的正确性。 本文的研究成果可用于二维和三维复杂构形中的非理想爆轰波传播及其相 互作用的数值计算，为进一步耦合流体动力学计算打下基础。 关键词：非理想爆轰，l e v e ls e t 方法，贴体坐标系，爆速曲率关系 三维非理想爆轰波阵面传播过程的数值方法和数值模拟 a b s t r a c t t h i sp a p e rs t u d yt h el e v e ls e tm e t h o df o rp r o p a g a t i o no fn o n i d e a ld e t o n a t i o n b a s e do nd s d ( d e t o n a t i o ns h o c kd y n a m i c s ) m o d e l i no r d e rt o t os i r e p l i f yt h e b o u n d a r yc o n d i t i o nt r e a t m e n t i nc a r t e s i a nc o o r d i n a t e ，b o d y - f i t t e dc o o r d i n a t ea r e e m p l o y e d w eg i v e a na l g o r i t h mo ft h eb o u n d a r yc o n d i t i o n a c c o r d i n gt ot h e h a m i h o n - j a c o b if o r m u l m i o no ft h eg o d u n o v ss c h e m e ，w es t u d yf i n i t ed i f f e r e n c e m o t h e df o rt h ep r o p a g a t i o ne q u a t i o no fn o n - i d e a ld e t o n a t i o ni nn o n 。o r t h o g o n a lb o d y f i r e dc o o r d i n a t ea n dp r e s e n ta na l g o r i t h m w ed e v e l o pt h ec o d ef o rs i m u l a t i n g p r o p a g a t i o no fn o n - i d e a ld e t o n a t i o ni n2 da n d3 da n dg i v eaf e wc a l c u l a t i o n s t h e r e s u l t sa r ei na g r e e m e n tw i t ha n a l y t i cr e s u l t si nb o d y f i r e dc o o r d i n a t e f o rb o d yo f r e v o l u t i o n ，t h e2 dc a l c u l a t i o n sa r ei na g r e e m e n tw i t h t h e3 dc a l c u l a t i o n s t h e a l g o r i t h ma n dt h ec o d ei nt h i sp a p e r h a v eb e e nv a l i d a t e d t h ep r o d u c t i o ni n t h i sp a p e rmf o rn u m e r i c a ls i m u l a t i n gp r o p a g a t i o na n d i n t e r a c t i o no fn o n i d e a ld e t o n a t i o ni n2 da b d3 d i ti st h eb a s ef o ri n t e g r a t i n gf l u i d d y n a m i c sm e t h o d k e y w o r d ：n o n i d e a ld e t o n a t i o n ，l e v e ls e tm e t h o d ，b o d y f i t t e dc o o r d i n a e t ， d 。r e l a t i o n s 独创性声明 y8 0 9 9 12 本人声明所呈交的学位论文是本人在导师t - t 4 - 。导下进簪的研究工作及取得的 趼究成果。据我所知，除了文中特别加以标注和致谢的地万外，沦文中不包含其 他人已经发表或撰写过晌研究成果，也不包含为获得中蹦1 ：程物理研究院或其他 教育机构的学位或证书使用过的材料。与我同工作的同志对本研究所做的任何 贡献均已在论文中作了明确的说明并表示谢意。 学位论文作者签名 针欲 签字同期：上刃岁年占县仞同 学位论文版权使用授权书 本学位论文作者完全了解并接受中国工程物理研究院研究生音# 有关保存、使 用学位论文的规定，允许论文被查阅、借阅和送交国家有关部门或机构，同时授 权中国工程物理研究院研究生部可以将学位论文全部或部分内容编入有关数据 库进行检索，可以影印、缩印或扫描等复制手段保存、汇编学位论文。 学位论文作者签名导师签名 签字日期：2 口涉年d 月o m 签字日期时年二月加日 三维非理想爆轰波传播过程的数值方法和数值模拟 1 绪论 复杂几何形状炸药中爆轰波阵面传播过程的描述是炸药驱动装置设计中的 一个重要课题。对于非理想爆轰，d s d ( d e t o n a t i o ns h o c kd y n a m i c s ) 理论给出 的爆速曲率关系是其传播规律的最简洁描述，可以将非理想爆轰波阵面传播问 题转化为曲率相关的界面运动问题。l e v e ls e t 方法作为一种描述运动界面的数 值方法，在流体动力学计算等方面的到了广泛的应用，是适合于处理曲率相关 的界面运动问题的数学方法。a s l a m l l l 采用l e v e ls e t 方法结合爆速曲率关系研 究了直角坐标系中的爆轰波阵面传播的l e v e ls e t 函数方程，能模拟爆轰波阵面 的运动，且实现了三维计算。但在直角坐标系中，弯曲的边界难以对应到网格 点上，边界条件的处理非常困难。陈森华，张旭【2 】等提出了采用等提出采用贴 体坐标系简化边界条件处理的方法，研究了二维贴体坐标系中的爆轰波传播问 题，可以处理较为复杂的边界，证明了方法的可行性。但采用文【2 】中给出的贴 体坐标系中的l e v e ls e t 函数方程的差分形式计算中心一点起爆模型时，预期为 圆形的波阵面上出现不合理的尖突，分析原因为该差分形式用于非正交网格中 的l e v e ls e t 函数方程时，仅适用于方程中的交叉项影响较小的情形。并且文 2 】 中未引入重新初始化方法，在计算绕射模型时，由于波传播方向的改变，计算 初始时所给出的l e v e ls e t 函数值已经失去意义，无法完成计算。同时，文【2 中 给出的处理边界条件所采用的拟合二次曲线方法在所选取的6 点插值多项式线 性相关时无法工作。因此，我们的研究目的为研究和发展二维和三维爆轰波传 播问题适用于贴体坐标系中非正交网格的l e v e ls e t 方法。通过选择合适的差分 格式，研究适用于非正交网格的计算爆轰波阵面传播问题的数值方法，引入重 新初始化方法解决绕射计算的问题，研究适合贴体坐标系的边界条件处理方 法。我们的研究内容包括三维问题，而三维问题的计算量巨大，需要引入提高 计算效率的方法。最后编制适用于二维和三维贴体坐标系中爆轰波阵面传播方 程的数值计算程序，建立测试模型以验证方法。 1 1 小曲率弯曲爆轰波传播规律 非理想爆轰波具有较厚的反应区，使爆速明显依赖于波阵面曲率。对于弯 曲爆轰波传播问题，采用忽略爆轰波厚度的c - j 爆轰模型，已不能准确描述波 阵面传播规律。d s d 理论指出，对于小曲率弯曲爆轰波，爆速仅仅是波阵面曲 率的函数。因此，爆速曲率关系可以认为是小曲率弯曲爆轰波传播规律的数学 表述。采用爆速曲率关系描述的波阵面传播，可以认为是一个界面位置决定界 三维非理想爆轰波传播过程的数值方法和数值模拟 面速度的界面运动问题。通过实验给出的爆速曲率关系的各个参数，数值计算 小曲率弯曲爆轰波传播问题可以转化为计算曲率相关运动界面的问题。 1 2 非理想爆轰波传播问题的计算方法 运动界面处理是e u l e r 坐标系中多介质流体动力学数值计算问题中的关键 之一，数值计算对运动界面的处理最初主要有两种方法。其一为示踪粒子法， 即以无质量的示踪粒子跟随界面运动，利用这些示踪粒子来确定界面位置的方 法。另一类方法为v o f ( v o l u m eo ff l u i d ) 方法。v o f 方法多用于捕获多介质 流体的界面。数值计算非理想爆轰波传播问题，可以采用示踪粒子法，在小扰 动的运动情形下，该方法有很高的精确度，但在较为复杂的界面运动问题上， 采用示踪粒子法会遇上较大的问题。首先，随着界面的运动，示踪粒子的分布 会变得过于密集或过于稀疏，需要采用重分技术增加或减少示踪粒子，而重分 技术不可避免的带来了计算误差。更大的问题是，当界面的拓扑结构发生改变 时，例如，爆轰波相遇时，使用示踪粒子法处理界面的变化会有很大的困难。 在1 9 8 8 年，o s h e r t l 】提出了处理运动界面的适合计算非理想爆轰波传播问题一种 新方法，称为l e v e ls e t 方法( 位标函数法) ，近年得到了广泛的应用。 1 3 l e v e ls e t 方法及国内外研究现状 文 3 中给出的l e v e ls e t 方法的基本概念是用l e v e ls e t 函数的零等值面来表 示物质的界面。对于随时间运动的界面，采用随时间变化的函数( i ，f ) 的零等 值面来表示，函数y e ，f ) 满足定的方程。在每个时刻t ，我们求出函数妒( i ，r ) 的值，从而得到当前时刻界面的位置。 函数g ，f ) 称为l e v e ls e t 函数，如图1 , 3 1 所示，对封闭的运动界面r o ) ， 有 f 妒g ，r ) o 当i 位于r ( f ) 夕 v ( i ，t ) 的值在满足上式的基础上有一定的任意性，可根据计算问题的需要进行 定义。但在大部分的应用中，缈的初值妒g ，o ) 一般取为主点到r ( 0 ) 的符号距 离。 一一 三丝韭里望堡薹鎏堡塑垫堡盟塑堕立蓬塑塑篁堕型 图1 3 1二维的l e v e ls e t 函数( 引自文【1 】) t u ( x ，y ，f ) = o 代表界面，0 ，y ，r ) o 代表界面内部区域 在给出界面运动的法向速度“。后，根据初始条件求解下面的h a m i l t o n j a c o b i 偏微分方程来确定 + 面v y = 0( 1 3 2 ) 式中y ，为】f ，关于时间的微分，由妒= 0 就可确定随时间变化的界面r ( f ) 。采用 l e v e ls e t 方法计算运动界面问题，无需直接追踪物质界面，提高了追踪复杂界 面的能力。并且，采用l e v e ls e t 方法后界面的拓扑变化不需要特殊处理就可以 通过_ f ，( i ，t ) 的变化自然得出分离或融合后的界面，( 图1 3 2 ) 。同时，界面 r ( f ) 表示为妒i ，f ) = 0 的隐函数形式，可以方便地计算界面的曲率，法方向等几 何参数，使采用该方法数值计算弯曲爆轰波传播问题较为方便。 三维非理想爆轰波传播过程的数值方法和数值模拟 图1 3 2拓扑变化的自然反映( 引自文1 4 】) l e v e ls e t 方法在推进时可以自然的处理阵面的拓扑变化 为构造直角坐标系中l e v e ls e t 函数方程的差分格式，o s h e r l 2 等研究了如下 数学问题： 。 f 紫i ? ( 1 3 3 ) ( x ，o ) = ( 工) 、。 式中，d y = 帆1 一，y 。，日0 l ，_ 一，) = 一- ；+ + “：产 对于一维情形，”= 虮有 坼+ 旧0 ) 】，= 0 ( 1 3 4 ) 如果g ，+ 坂= g ( ，。，“m + 。) 满足l i p s c h i t z 条件和相容性要求 g ( ”一，“) = 日莓) ，并且能构成对所有变量“不减的差分格式 “= “? 一( g ，+ 一g ，必) ( 1 3 5 ) 则方程( 1 3 3 ) 的解由下述差分格式给出 矿? “= y ? 一a t g ( d _ l ，+ 1 ，d + y k ) ( 1 3 6 ) 式中d 一? = ? 一y 三。，d + j ；f ，? = 盘- 一? 方程( 1 3 2 ) 可写为 对于方程( 1 3 7 ) 形式。 + ”。i v 1 = 0 ( 1 3 7 ) o s h e r 等给出了直角坐标系中的高阶和多维情况下的差分 在自由面追踪或二相流问题中， 采用以下形式 考虑有r ( f ) = 每g ，f ) = o ) ，方程( 1 3 2 ) 可以 + “虬+ ”虬= 0 ( 1 - 3 8 ) 对于方程( 1 3 8 ) ，s u s s m a n 5 ，6 1 等人采用l e v e ls e t 方法计算直角坐标系下的 不可压两相流动问题，文献【6 中采用了高阶e n o 格式的计算方法，计算了二 维三维轴对称液气流动问题，并且考虑了表面张力和粘性的影响。在此基础 上，文献【7 】中s u s s m a i l 等人还将l e v e ls e t 方法和自适应计算方法相结合，提高 4 三维非理想爆轰波传播过程的数值方法和数值模拟 计算精度。在s u s s m a n 等的两相流动问题计算方法中，将l e v e ls e t 函数的零等 值面设为液气分界面，取符号距离函数为l e v e ls e t 函数。因为假定了流体的不 可压性，因此在气体和液体中，密度都为常数，在分界面上产生了间断。 s u s s m a n 等人利用在界面的邻近区域l e v e ls e t 函数可视为计算网格到物质分界 面的距离这一点将液，气分界面临近区域的密度间断处理为连续密度，使相关的 流体动力学方程可以求解。 f e d k i w 和a s l a m 8 】等将l e v e ls e t 方法和g h o s tf l u i d 方法结合，计算多物质 流动的界面运动问题。l e v e ls e t 函数确定了物质分界面，也就指出了r e a l 网格 与g h o s t 网格，在曲o s t 网格上的压力、速度等认为在界面两侧是连续的物理量 可以采用g h o s t 网格上真实流体的相应物理量，但界面被认为是一个接触间断， 熵在界面两侧是不连续的。为了给出曲o s t 网格上的熵值，文献f 7 】采用了 i s o b a r i c 血技术，界面的法线方向可由l e v e ls e t 函数确定，即 雨= 呙 ( 1 3 9 ) v 驯 、7 则界面附近的熵值立调整为常数，由以下方程确定 i ：n v i = 0( 1 3 1 0 ) a s l a r n 9 1 采用l e v e ls e t 方法处理冲击波传播问题，利用l e v e ls e t 函数确定 间断面区分波后区域和波前区域，在间断面的法线方向上求解r i e m a n n 问题来 确定冲击波传播速度，并为g h o s t 网格赋值，然后分别求解波前和波后区域的欧 拉方程，进而完成整个区域的计算。a s l a m 等还成功地利用l e v e ls e t 方法计算 了立方体炸药多点起爆问题，其控制方程使用的是直角坐标系中的基本方程。 国内对l e v e ls e t 方法的应用主要用在流体力学计算问题中处理多介质界面 问题。柏劲松【1 0 1 等提出了一种处理三种以上介质的情况下消除三重点空白区的 方法。将l e v e ls e t 方法用于爆轰波传播问题并不多见，文尚刚瞎研究了直角 坐标系中的爆轰波传播问题。陈森华，张心1 , 1 2 , 1 3 等研究了二维贴体坐标系下的 爆轰波传播问题。 l e v e ls e t 函数的取值在满足定义的条件下有一定的任意性，但符号距离函 数有较明确的意义。大多数问题中l e v e ls e t 函数采用的是符号距离函数，在计 算过程中，由于数值方法的影响，l e v e ls e t 函数会偏离符号距离函数。为使 l e v e ls e t 函数保持一个良好的性质，基于o s h e r 和s e t h i a n 的工作，p e n g 口4 等人 提出了重新初始化的方法，通过迭代求解下面的h a m i l t o n j a e o b i 偏微分方程， ( 1 3 1 1 ) 得到符号距离函数d 来实现重新初始化，然后令相应计算节点上的 v ：d 将l e v e ls c t 函数重置为保持波阵面位置的符号距离函数。 d rd - s o ) _ g n 然芝t ；) ， id b ，o ) = 玩( x ) = _ f ，b ，) ” 三维非理想爆轰波传播过程的数值方法和数值模拟 计算中s 毋z p ) 光滑化为 s i g n ( d ) = s 。0 ) = 了兰( 1 3 1 2 ) q d + s 文【1 4 等给出了直角坐标系中重新初始化( 1 - 3 1 1 ) 方程的差分方法。 采用计算l e v e ls e t 函数的方法来计算界面的运动，需要额外的计算量。在 计算模型较大的时候，如果在整个流场中计算l e v e ls e t 函数的值，消耗的计算 时间过长，导致计算效率低下，并且在距离界面过远的位置，难以给出合适的 法向速度。而从方程( 1 3 2 ) 可看出，在数值计算界面位置时，下一时刻的界面 位置由当前时刻界面附近区域的网格上的l e v e ls e t 函数值和界面的法向速度确 定。而在应用中，所关心的l e v e ls e t 函数值往往只在界面附近的一个狭窄的带 状区域内。可以想见如果采用窄带计算方法，即只计算界面附近的l e v e ls e t 函 数值并保持界面运动不受影响，应该可以极大地提高计算速度。在文 1 4 1 中， p e n g 等人引入常数0 卢 y ，通过截断函数c ) 1 当m 卢 c 眵) = 0 j i c ，i 一7 ) 2 ( 2 1 妒1 + y 一3 n ) p 一卢) 3 当卢 ， 将l e v e ls e t 函数y 的偏微分方程( 1 3 7 ) 变化为 + c 如。i v 纠= 0 ( 1 _ 3 1 4 ) 这样，将计算区域限制在一个由y 决定狭窄的范围内，经过一个时间步长的计 算后，l e v e ls e t 函数已经偏离了符号距离函数，这时需在一个较大一些的区域 内进行重新初始化，将l e v e ls e t 函数重置为符号距离函数d ，最后按( 1 3 1 5 ) 式 截断d 得到下一时间步计算需要的l e v e ls e t 函数初值妒( 图1 3 3 ) 。 l 一卢ys p y = d一卢 d 卢 三维非理想爆轰波传播过程的数值方法和数值模拟 a b 图1 3 3窄带计算的实现( 引自文【9 1 ) 实线为当前步骤的结果，虚线为前面步骤地结果a ) 计算初始的l e v e ls e t 函数。b ) 推进一个 时间步。c ) 经过重新初始化。d ) 截断后的l e v e ls e t 函数，下一时间步的初始值。 l e v e ls e t 方法作为一种能准确描述界面位置( 界面为无厚度的几何面，对 它的形状不做任何数学上的假定) ，并且可以方便地计算界面的诸如曲率等几 何性质，自出现以来在流体动力学计算中得到了广泛的应用。该方法还被用于 其他研究领域，如图象处理，机器人设计等。l e v e ls e t 方法在流体动力学计算 中主要用于运动界面的处理，一般取符号距离函数为l e v e ls e t 函数，在界面附 近可以表示和界面的距离。对l e v e ls e t 方法的研究首先是将该方法与其他流体 动力学计算计算方法相结合，扩展其应用范围，另外就是提高l e v e ls e t 方法的 计算精度和速度。 1 4 边界条件与贴体坐标系 爆轰波传播问题中，很多情况下不可避免会遇上炸药边界上的波阵面。在 波阵面边界点，边晁稀疏不能进入化学反应区的物理条件，在爆轰产物状态方 程为多方指数的情况，可归结为波阵面外法线与边界外法线的夹角小于依赖炸 药的角的几何边界条件【13 1 ，一般的形式为： 三维非理想爆轰波传播过程的数值方法和数值模拟 1 等一，v _ ) f ，d 。i 0 ( 1 4 1 ) o n 式中的y _ x k 为边界网格点上的沿边界外法线的方向导数。 d 门 a s l a m 【1 】等在计算立方体炸药多点起爆问题时，在炸药形状为规则的立方体 时得到了满意的结果。但对于较为复杂的炸药边界，如果在直角坐标系下进行 处理，一般的炸药边界不会位于网格边上，文 1 中处理方法非常繁琐。对于这 样的问题，为简化边界条件的处理，陈森华，张旭f 2 】等提出采用贴体坐标系的 方法，计算的网格边界即为炸药边界，炸药边界的几何性质蕴含于计算网格 中，可以处理较为复杂的边界条件。但文 2 】中仅考虑了二维情况，并且其中 l e v e ls e t 方程的差分格式尚不完善。文 2 中也没有采用重新初始化的方法将 l e v e ls e t 函数重置为符号距离函数，在计箅有绕射情况的爆轰波传播问题时会 由于波传播方向发生变化而无法得到合理的波阵面。 1 5 本文的主要研究内容及特色 本文主要研究将l e v e ls e t 方法用于数值模拟采用爆速曲率关系描述的小曲 率弯曲爆轰波传播问题。主要的研究内容及其特色如下。 ( 1 1 、发展了a s l a m 的研究工作。a s l a m 等已经成功地利用l e v e ls e t 方法 计算了立方体炸药多点起爆问题。但是他们使用直角坐标系中的基本方程，得 到的差分格式很难应用到曲面边界的炸药，对于三维曲面边界情况困难更大。 本文的方法简化了三维曲面边界的处理。 ( 2 ) 、推广了国内同行的研究工作i 2 1 。陈森华，张魍的方法计算了二维曲边 炸药中的爆轰波传播问题，但末推广到三维，并且其中的差分形式和数值计算 方法用于非正交网格时，仅适用于l e v e ls e t 函数方程中的交叉项影响较小的情 况。本文改进了二维的数值计算方法，并且推广到了三维情况。 ( 3 1 、给出二维和三维贴体坐标系下的适用于非正交网格的l e v e ls o t 计算方 法。本文从h a m i l t o n j a c o b i 偏微分方程的一般差分形式出发，给出贴体坐标系 下l e v e ls e t 函数方程和重新初始化方程的具体的差分形式。给出的差分形式适 用于二维和三维，没有非物理振荡，可以应用于非正交的网格。同时讨论了适 用于二维和三维贴体坐标系下具有复杂边界的爆轰波传播问题的几何边界条件 处理方法。 最后，本文讨论了采用给出的差分形式和边界条件处理方法的爆轰波传播 问题f o r t r a n 语言数值计算程序试算算例的计算结果，并探讨了与其他流体动力 学计算程序的鹚合应用。 三维非理想爆轰波传播过程的数值方法和数值模拟 2 小曲率爆轰波波阵面传播问题的爆速曲率关系【1 2 1 2 1 动力学基本方程 一维爆轰化学反应流可用以下e u l e r 方程组描述： 鱼+ d p u + z p u ：o 8tarr 要肛+ 要( 2 + p ) + 丛：o 瓦肛+ 万( 1 p ) + 芋。0 丢肛+ 鲁郇+ 争 孕+ u 孕：i ( 旯，舢) 一十= r t 几，p 口i atar ”一。 丝唑幽：o 2 川) 方程中p 、1 4 、p 、和r 分别为密度、质点速度、压力和e u l e r 坐标；五为 反应进程变量，未反应炸药中旯= 0 ，反应结束时五= 1 ；( 咒，p ，p ) 为化学反 应速率方程；是几何因子，a = o ，1 ，2 分别对应平面、柱面和球面情况；e 是 总能量e = q + “2 一坦，式中e t 是热力学能，o 是炸药的反应热。 对小曲率曲面爆轰波，反应区内流动是拟稳定的，设先导冲击波阵面坐标 为，= r s ( f ) ，爆速为硼2 = d s ( r ) 。可用波阵面中值曲率匿表示： = 2 k 。引进坐标变换x = r ；一r ，般 l 时，略去石的高阶量后方程( 1 ) 变 化为： d _ 哑p ，v ，p ，e ) = l c f ( p ，v ，p ，e ) 靠 m ( 2 1 2 ) v 拿：t ( 九膨p ) 甜 蛾唧，驴 菇三葛， ，聊，唧= 一：z - 2 p 嚣( d , 秽- v ) 二，) c z s ， 式中v = v s 一“。总能量仍表示为e ，但与( 2 1 1 ) 中的定义略有不同，此处的定 义为e = 巳+ 必v 2 一蚀。 9 三维非理想爆轰波传播过程的数值方法和数值模拟 塞= ( 舅 肼塞+ 滢l 害+ ( 鲁) 。警 c z _ 。， 内能e = e 。一旭。从方程( 2 1 2 ) 解得车、宰、! 拿、掣代入方程( 2 i ，4 ) 可得 a , x 靠甜a x 妾= 盟掣1 ) 出一m ： 、7 式中肘，2 为m a c h 数，c ，为冻结声速，以为热性系数，它们的定义为 弘鬻印一b 插， 社d 2 眷警 亿：。， 硪 r o ( z ，p ，p ) 由方程( 2 2 1 ) 积分得到积分方程 中( p ( 五) ，v ( z ) ，p ( 五) ，e ( 咒) ) = 中( p o ，d 。，p o ，e o ) + s ( 五) ( 2 2 2 ) 下标o 表示未反应炸药中的参数，源项s ( 丑) 是积分 鼢) - k i 燮蹉铲赢 ( 2 z 3 ) 假定在a = 1 附近一具有下形式 名= ( 1 一a ) 。f ( p j ，乃) + o ( 0 一五) 。) ( 国0 ) ( 2 2 4 ) j 表示c j 点的参数，国是化学反应阶数。 求解常微分方程( 2 ，2 1 ) 的初值问题时， = 0 处的流场参数初值可由任意给 定的冲击波速度d ，( 一。 d 。 0 是一个不依赖于k 的常数。不难看出，( 2 2 7 ) 给出的可阻保证上述两个 条件熊褥到满足。将( 2 2 7 ) 代入( 2 2 。6 ) 得 s f ) = 鑫( 2 2 。8 ) 三维非理想爆轰波传播过程的数值方法和数值模拟 式中 占= r + o ( j r ) 一芷i n k + 0 0 , i n 盯、 r 。+ 口。) ( o 1 1 7 鼍彰筹铲西( 0 1 时，的分 量都包含因子七p ，0 t 与t p 成正比，仪不能被唯一确定。 根据爆速偏离值的定义，将( 2 2 9 ) 代入( 2 2 1 4 ) 目p n 得到爆速曲率关系 d ，= k s d ，+ 口, 2 + o ( k )( 0 1 ) 爆速对于0 ) 的依赖是连续的，所以对于0 ) 1 应当有 厂 土 土、 1 i i n i 口( 国) 茁。+ o ( 茁。) l = 一口( 1 ) r l n r + o ( r l n 芷) ( 2 2 1 9 ) 一l 从( 2 2 1 0 ) 可知，当 1 时系数仪包含了一个因子( 1 一) 一。从( 2 2 1 9 ) 可以估计 出d ( k 坞) = ( p 一1 c + d ( 1 ( ) 。此处，p 与0 l 相似，是依赖于和其它炸药参数而 不依赖于k 的系数。对于f 0 l ，从( 2 2 1 0 ) 可知，当呻1 时系数值的渐近表达 式包含了一个因子( 1 ) 。可以按同样的办法估计爆速曲率关系中的高阶项。 最后得到爆速曲率关系 d s = 七。d o , + 口j f + ( 卢一口) k + o ( 茁) 一d r l n 茁+ 芦茁+ 0 0 c ) 1 觥。+ ( 一口) 盯+ d ( _ i r ) ( o 功 】) 三维非理想爆轰波传播过程的数值方法和数值模拟 3 贴体坐标系下非理想爆轰波传播方程的差分形式 3 1 贴体坐标系下的基本方程。 d s d 理论给出的爆速曲率关系是对弯曲爆轰波传播规律的最简洁描述，应 用爆速曲率关系数值计算小曲率弯曲爆轰波传播爆轰波阵面传播问题就转化为 计算界面运动速度与界面当地曲率相关的界面运动问题，计算所需的参数由实 验给出，计算波阵面位置时不涉及具体状态方程形式，具有较强的通用性。 o s h e r 等提出了计算曲率相关的界面运动问题的l e v e ls e t 方法，在直角坐标系 中波阵面满足方程y g ，f ) = 0 ，l e v e ls e t 方法的主要工作是根据初始条件求解确 定v 的偏微分方程，方程可写为 + 见仁】v j = 0 ( 3 1 1 ) 式中d 。体) 是阵面传播速度在法线方向的分量，k 是阵面的中值曲率。我们总 是取法线线方向和爆轰波传播的方向相同，因此对( 2 2 2 0 ) 式有一些变化， d 。伍) 可写为： 或( k ) = d 。+ 占( 茁) ( 3 1 2 ) 为简化计算，忽略高阶项，取6 ( x ) = 一a j f ，( 3 1 2 ) 成为 v - ，+ d 。i v 1 = a r i v p l ( 3 1 3 ) 为叙述简单起见，这里以二维为例，推导计算公式可推广到三维情况。 在二维情况下引进坐标变换【2 】 f 2 f ( x ，一) r 3 1 4 、 叩= g ( x ，y ) 、 其中变换应满足两个基本要求： ( 1 ) 为了保证逆变换存在，应当有j a c o b i 行列式j = d 皓g ，y ) o ； ( 2 ) 在计算空间皓，叩) 边界即炸药边界上，应当有f = 茧或7 7 = 叩，计算区域是矩 形。 于是有 中x2 妒； ? x + l f ，q g x vy = 中；f p 七l l ，q g y 只、g ，、g ，利用j a c o b i 矩阵可得出。此时 m ：再万瓦习q 话丽 展开合并同类项，可以写为 l v i = q p f ) 2 + d ：眠) 2 + 2 口。p ；p 。 ( 3 1 5 ) ( 3 1 6 ) ( 3 1 7 ) 1 4 三维非理想爆轰波传播过程的数值方法和数值模拟 式中 a l = ( t ) 2 + 眩) 20 2 = ( g ，) 2 + ( g ，) 2 口。= c g 。+ g ， 于是方程( 3 1 3 ) 在二维贴体坐标系下可以表示为 妒：+ h ( v f ，y 目) = 髓r i v _ ；f ，l ( 3 1 8 ) 式中 。日( “，v ) = p 。口。0 ) 2 + 日：( v ) 2 + 2 a 。v ( 3 1 9 ) 当网格为非正交网格时，( 3 1 9 ) 式中的a 。项不为零。( 3 1 8 ) 式的左边出现了交叉 项，不能实现变量分离，大部分文献提供的直角坐标系中的差分计算方法无法 应用。 3 2 。曲率部分戆计算 考虑方程( 3 1 8 ) ，科( ，) 仅含有对空间嫩标的一阶偏导数，方程的左 逮灸一个h a m i l t o n - j a c o b i 攘徽分方程，方穗熬右边窘弯謦骜窭豹鞭，戆率霉要 用对空闽嫩标的一酚釉二阶俯学数表示，不能交承成为h a m i l t o n - j a c o b i 偏徽 努方程。我们对方程戆溪逸分秀诗舞，方羧靛砉逸可以采瘸中，差分谤箕，嚣 w 为l e v e ls e t 函数，计弊所用的曲率可采用隐函数曲率计算公式，即 一( 尚 b z - ， 三维情况下的曲率计算公式由( 3 2 。1 ) 式直接给出，如果采用二维程序计算 三维轴对称模型，需要将( 3 2 1 ) 式计算出的二维曲率转换为旋转体的三维曲 i f ，( ，z ) = 0r = x 2 + _ y 2 ( 3 2 2 ) ；三r ，；上 。鼍1 7 x y 1 y ( 3 2 3 ) rr ，10 气、 x 2 2 ， o 2 了一7 一7 2 i 一7 妒，= ￥，y ，= 妒，0 ( 3 2 4 ) 三维非理想爆轰波传播过程的数值方法和数值模拟 “= y ，( ) 2 + _ ；f ，k 。= 少。h ) 2 + _ 】f ， 将( 3 2 3 ) ( 3 2 5 ) 代入三维曲率计算公式，得 芷= 业盟巫型罐筹学业型业趔 2 i v ：+ j + v ：r 1 r ：生嵩掣 b z 匀 r = 可一 l o o ， 2 + ；尸2 在以y 坐标轴为旋转轴的旋转体模型中，考虑x 为正的区域，于是在( 3 2 5 ) 中， 用g ，_ y ) 替换( 3 2 5 ) 中r ，z ) ，结合二维的( 3 2 1 ) 式，可以得到下面的旋转体曲率 加 2 ( 3 2 6 ) 式中，位于旋转轴上的点有z = 0 ，对于这样的点 等同于三维曲率。 3 3 h a m i l t o n j a c o b i 偏微分方程的差分方法 ( 3 2 6 ) 不难证明其二维曲率 方程( 3 1 8 ) 的左边为一个h a m i l t o n - j a c o b i 偏微分方程，可以采用h a m i l t o n 。 j a c o b i 方程的通用差分方法，考虑方程 f ，+ 日( 恢，) = 0 ( 3 3 t 1 ) 差分计算需将方程( 3 3 1 ) 数值离散，o s h e r 和f e d k i w u 5 1 给出了离散方法，对 方程( 3 3 1 ) ，时间采用e u l e r 向前差分的数值离散后可以表示为 竺半+ 存”，p ；崂，y ；) = o ( 3 蚴 ( 3 3 2 ) 式即( 1 3 5 ) 式，h ”p ；，；崂，；) 需满足相容性条件 h ”f ，y ；y 。，) = 日妙f ，y ，) ( 3 - 3 3 ) 其中妒一等项可以为一阶精度的向前或向后差分，也可以采用二阶精度的 e n o 格式和更高阶的格式，时间离散还可以采用高阶的r u n g ek u t t a 格式， 三维非理想爆轰波传播过程的数值方法和数值模拟 y ”1 还需满足单调性条件，贴体坐标系中( 3 3 2 ) 式的妙”1 满足单调性条件由 ( 善，7 7 ) 坐标系下的