基于坐标变换的强间断问题伪弧长算法_陈泽平.pdf

第 44 卷第 4 期兵 器 装 备 工 程 学 报2023 年 4月收稿日期:2022 05 23;修回日期:2022 06 13基金项目:国家自然科学基金项目(11822203)作者简介:陈泽平(1995)，男，硕士研究生，E-mail:chernzp_9637 outlook com。通信作者:马天宝(1981)，男，博士，教授，E-mail:madabal bit edu cn。doi:1011809/bqzbgcxb202304011基于坐标变换的强间断问题伪弧长算法陈泽平，王晨涛，李坤，马天宝(北京理工大学 爆炸科学与技术国家重点实验室，北京100081)摘要:针对爆炸冲击波数值模拟中经常要面临的强间断奇异性问题，在一维或多维空间中引入弧长参数，通过添加一个约束方程，致使原网格自适应移动变形而能够在强间断区域分布较多的网格点，以此提高大梯度物理量捕捉的分辨率。针对非均匀变形物理空间中格式构造困难的问题，提出了坐标变换，将物理空间中的控制方程映射至均匀正交的弧长计算空间。基于 MUSCL 格式编写了伪弧长算法相应的一维、二维程序，并对其开展人为解方法验证。将伪弧长算法(PALM)的数值结果同有限体积法进行比较，经过误差分析，验证了伪弧长算法在间断附近具有较高的分辨率。关键词:爆炸冲击波;数值模拟;伪弧长算法;坐标变换;精度验证本文引用格式:陈泽平，王晨涛，李坤，等 基于坐标变换的强间断问题伪弧长算法 J 兵器装备工程学报，2023，44(4):68 76Citation format:CHEN Zeping，WANG Chentao，LI Kun，et al esearch on pseudo arc-length algorithm for strong discon-tinuity based on coordinate transformation J Journal of Ordnance Equipment Engineering，2023，44(4):68 76中图分类号:O381文献标识码:A文章编号:2096 2304(2023)04 0068 09esearch on pseudo arc-length algorithm for strong discontinuitybased on coordinate transformationCHEN Zeping，WANG Chentao，LI Kun，MA Tianbao(State Key Laboratory of Explosive Science and Technology，Beijing Institute of Technology，Beijing 100081，China)Abstract:Aiming at strong discontinuity singularity often occurring in the numerical simulation ofexplosion shock wave，this paper introduces arc-length parameters in one-or multi-dimensional space Byadding a constraint equation，the original grid can move and deform adaptively to achieve distribution of agreat number of grid points in strong discontinuous areas so as to improve the resolution of large-gradientphysical quantity capture Concerning the difficulty in constructing a scheme in non-uniform deformationphysical space，a coordinate transformation is proposed to map the control equation in physical space intothe uniform orthogonal arc-length computational space In addition，the corresponding one-and two-dimensional programs of the pseudo arc-length algorithm are compiled based on MUSCL，and the method ofmanufactured solutions is developed to verify this algorithm The numerical results of the pseudo arc-lengthalgorithm are compared with the finite volume method Through error analysis，it is verified that the pseudoarc-length algorithm has high resolution near the discontinuityKey words:explosion shock wave;numerical simulation;pseudo arc-length algorithm;coordinatetransformation;accuracy verification0引言爆炸与冲击是在高温高压和相变等极端条件下，气液固多介质间强耦合作用的瞬态动力学问题，数值求解该问题模型对航空航天等国防工业及武器装备的研制开发均具有重要的基础应用价值1 2。计算爆炸力学很重要的一个问题在于精确捕捉爆轰波的波阵面及波阵面前后物理量的变化，但对于爆炸冲击波的双曲型守恒律方程，其解随着时间的演化往往会出现奇异性即存在强间断，例如激波、接触间断，或者稀疏波之类的弱间断。为了降低这种奇异性，提高间断分辨率，近年来发展了许多的理论和计算方法，像谱方法、奇异摄动理论、小波分析法3 以及某些高分辨率的数值格式等。在早期，数值计算捕捉激波通常采取一阶精度的差分格式，后来，Van Leer4 将一阶 Godunov 格式进行推广，率先提出了二阶精度 MUSCL(monotone upstream-centred schemes for conservation laws)，随后近 40 年，高精度、高分辨率数值格式蓬勃发展，出现诸如TVD、PPM、ENO、KDG、CE/SE、WENO、WENO-Z 及各类杂交格式等。高精度数值格式一般色散较强，容易产生非物理性振荡，通常需要额外的限制振荡的方法，比如人工粘性法，可人工粘性的添加有时会使数值解的耗散比较严重，间断的分辨率不够清晰，其也非从根本上解决振荡。此外，提升 Eulerian 法数值求解精度的另一个角度则是网格加密。若采用均匀网格计算求解，便需对整个计算域进行加密，使得计算资源极大地浪费，基于这个矛盾，网格自适应技术应运而生。网格自适应技术大致分三类，h-型(额外增加节点)、p-型(增加逼近多项式的阶数)和 r-型(移动网格节点来达到网格重新分布)，此外还有正在发展的结合性方法 h-p 方法及 h-r 方法等，这里着重介绍 r-型，也即移动网格方法。移动网格法发展至今，其所凭借的网格移动策略及网格泛函逐渐多样化，比如基于变量扩展的等势方法、调和映射、坐标变换的雅可比矩阵思想、基于所谓等分布和对齐条件的泛函等。最近几年，Luo 等5 又将 DG(discontinuousgalerkin)方法同移动网格偏微分方程法相组合，提出了一种针对双曲守恒律方程的拟拉格朗日移动网格间断 Galerkin法，从旧网格到新网格的物理变量并不需要插值;Lopez 等6 还提出一种基于 MMPDE 法的并行变分网格改进方案。本文中的伪弧长算法可归结为 r-型方法。由于伪弧长算法涉及网格的自适应移动，进而导致原始均匀正交的物理空间发生扭曲变形，这给格式的重构与插值带来了困难。为了避免直接在变形的非结构网格中重构数值格式(需要大量的模板，计算效率较低)，根据弧长映射关系，借助坐标变换将物理空间映射至均匀正交的弧长计算空间，然后在弧长计算空间中，基于维数分裂思想可以用较少的模板完成高精格式的重构，从而保证了伪弧长算法的计算效率。伪弧长算法能够将物理空间中的强间断问题转变成均分弧长空间中正常的弱奇异流体问题，而计算空间的转换更保证了原本数值格式的运用，在此过程中还能巧妙地避开虚假振荡。算法程序的可靠性一直是数值模拟领域亟待解决的重难点，诸如数理模型的简化、边界条件的设置近似以及迭代格式的选取都有可能对计算结果造成偏差7。人为解方法(method of manufactured solution，MMS)早期经典工作可现于Oberkampf、Trucano 及 oy 等。Blais 等8 提出了一种针对VANS 方程(the volume-averaged Navier-stokes equations)的人为解方法框架以验证其算法程序，并能同任何 CFD 技术相结合。Choudhary 等9 介绍了基于旋度而满足无散度约束的人为解方法，以验证两相不可压缩流控制方程。刘学哲等10 利用其构造的一类二维人为解模型针对性地验证辐射流体力学程序，该办法中动量、能量方程包含源项而质量方程无源项。本文中研究了二维情况下基于 MUSCL 的伪弧长算法模型建立过程，针对网格自适应移动造成的物理空间扭曲变形而不易构造高精度格式的难题，基于坐标变换的思想，将变形的物理空间映射至一套正交的弧长计算空间，从而在弧长计算空间中实现控制方程的求解与新网格守恒变量的插值过程，提高了计算效率。编写了该算法相应的一维、二维程序，在人为解方法基础上验证其精度，借助某些经典算例，将PALM 数值结果对比分析于有限体积法。1伪弧长算法本节伪弧长算法包含了 2 个部分。第一部分，借助有限体积法给出物理空间中控制方程的时空离散格式。如前描述，在考虑网格尺度影响的非均匀网格下，格式重构不易，尤其在二维及其以上的情况，因此，第二部分在网格自适应移动之后，采取坐标变换的策略将物理量映射至弧长计算空间中，并在该空间进行格式重构，求解完之后再逆变换映射回原物理空间。1 1离散格式考虑如下双曲守恒系统:wt+xF(w)+yG(w)=S(w)(1)式中:w 为质量、动量、能量组成的守恒变量;F(w)，G(w)，S(w)为 w 的函数。式(1)在网格单元 Ki，j上进行积分，得到:Ki，jwtd+Ki，jF(w)x+G(w)()yd=Ki，jS(w)d(2)式中:d 为面积微元。96陈泽平，等:基于坐标变换的强间断问题伪弧长算法设?wi，j为网格单元的物理量均值，利用 Green-Gauss 定理对式(2)进行变换，有:Ki，jt?wi，j+Ki，j(w)ni，jds=Ki，jS(?wi，j)(3)式中:Ki，j为网格 Ki，j的面积;Ki，j为网格的边界;ds 为网格边界长度微元;ni，j为边界Ki，j的单位外法向量;(w)=(F(w)，G(w)。数值通量借助局部 Lax-Friedrichs 格式，定义为:h(u，v，n)=12(u)n+(v)n c(v u)(4)式中:u 和 v 为守恒变量的内外重构值;n 为网格每一条边的单位外法向量;c=maxu，v(u)n。式(4)应符合守恒性跟相容性:h(u，v，n)=h(u，v，n)，h(u，u，n)=(u)n(5)式(3)写成半离散格式，有:Ki，jt?wi，j=4k=1kKi，jh(wint(k)i，j，wext(k)i，j，nki，j)ds+Ki，jS(?wi，j)=4k=1h(wint(k)i，j，wext(k)i，j，nki，j)