基于坐标变换的强间断问题伪弧长算法_陈泽平.pdf
《基于坐标变换的强间断问题伪弧长算法_陈泽平.pdf》由会员分享,可在线阅读,更多相关《基于坐标变换的强间断问题伪弧长算法_陈泽平.pdf(9页珍藏版)》请在咨信网上搜索。
1、第 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)摘要:针对爆炸冲击波数值模
2、拟中经常要面临的强间断奇异性问题,在一维或多维空间中引入弧长参数,通过添加一个约束方程,致使原网格自适应移动变形而能够在强间断区域分布较多的网格点,以此提高大梯度物理量捕捉的分辨率。针对非均匀变形物理空间中格式构造困难的问题,提出了坐标变换,将物理空间中的控制方程映射至均匀正交的弧长计算空间。基于 MUSCL 格式编写了伪弧长算法相应的一维、二维程序,并对其开展人为解方法验证。将伪弧长算法(PALM)的数值结果同有限体积法进行比较,经过误差分析,验证了伪弧长算法在间断附近具有较高的分辨率。关键词:爆炸冲击波;数值模拟;伪弧长算法;坐标变换;精度验证本文引用格式:陈泽平,王晨涛,李坤,等 基于坐
3、标变换的强间断问题伪弧长算法 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
4、 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 singula
5、rity 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 disc
6、ontinuous 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-
7、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
8、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引言爆炸与冲击是在高温高压和相变等极端条件下,气液固多
9、介质间强耦合作用的瞬态动力学问题,数值求解该问题模型对航空航天等国防工业及武器装备的研制开发均具有重要的基础应用价值1 2。计算爆炸力学很重要的一个问题在于精确捕捉爆轰波的波阵面及波阵面前后物理量的变化,但对于爆炸冲击波的双曲型守恒律方程,其解随着时间的演化往往会出现奇异性即存在强间断,例如激波、接触间断,或者稀疏波之类的弱间断。为了降低这种奇异性,提高间断分辨率,近年来发展了许多的理论和计算方法,像谱方法、奇异摄动理论、小波分析法3 以及某些高分辨率的数值格式等。在早期,数值计算捕捉激波通常采取一阶精度的差分格式,后来,Van Leer4 将一阶 Godunov 格式进行推广,率先提出了二阶
10、精度 MUSCL(monotone upstream-centred schemes for conservation laws),随后近 40 年,高精度、高分辨率数值格式蓬勃发展,出现诸如TVD、PPM、ENO、KDG、CE/SE、WENO、WENO-Z 及各类杂交格式等。高精度数值格式一般色散较强,容易产生非物理性振荡,通常需要额外的限制振荡的方法,比如人工粘性法,可人工粘性的添加有时会使数值解的耗散比较严重,间断的分辨率不够清晰,其也非从根本上解决振荡。此外,提升 Eulerian 法数值求解精度的另一个角度则是网格加密。若采用均匀网格计算求解,便需对整个计算域进行加密,使得计算资源极
11、大地浪费,基于这个矛盾,网格自适应技术应运而生。网格自适应技术大致分三类,h-型(额外增加节点)、p-型(增加逼近多项式的阶数)和 r-型(移动网格节点来达到网格重新分布),此外还有正在发展的结合性方法 h-p 方法及 h-r 方法等,这里着重介绍 r-型,也即移动网格方法。移动网格法发展至今,其所凭借的网格移动策略及网格泛函逐渐多样化,比如基于变量扩展的等势方法、调和映射、坐标变换的雅可比矩阵思想、基于所谓等分布和对齐条件的泛函等。最近几年,Luo 等5 又将 DG(discontinuousgalerkin)方法同移动网格偏微分方程法相组合,提出了一种针对双曲守恒律方程的拟拉格朗日移动网格
12、间断 Galerkin法,从旧网格到新网格的物理变量并不需要插值;Lopez 等6 还提出一种基于 MMPDE 法的并行变分网格改进方案。本文中的伪弧长算法可归结为 r-型方法。由于伪弧长算法涉及网格的自适应移动,进而导致原始均匀正交的物理空间发生扭曲变形,这给格式的重构与插值带来了困难。为了避免直接在变形的非结构网格中重构数值格式(需要大量的模板,计算效率较低),根据弧长映射关系,借助坐标变换将物理空间映射至均匀正交的弧长计算空间,然后在弧长计算空间中,基于维数分裂思想可以用较少的模板完成高精格式的重构,从而保证了伪弧长算法的计算效率。伪弧长算法能够将物理空间中的强间断问题转变成均分弧长空间
13、中正常的弱奇异流体问题,而计算空间的转换更保证了原本数值格式的运用,在此过程中还能巧妙地避开虚假振荡。算法程序的可靠性一直是数值模拟领域亟待解决的重难点,诸如数理模型的简化、边界条件的设置近似以及迭代格式的选取都有可能对计算结果造成偏差7。人为解方法(method of manufactured solution,MMS)早期经典工作可现于Oberkampf、Trucano 及 oy 等。Blais 等8 提出了一种针对VANS 方程(the volume-averaged Navier-stokes equations)的人为解方法框架以验证其算法程序,并能同任何 CFD 技术相结合。Cho
14、udhary 等9 介绍了基于旋度而满足无散度约束的人为解方法,以验证两相不可压缩流控制方程。刘学哲等10 利用其构造的一类二维人为解模型针对性地验证辐射流体力学程序,该办法中动量、能量方程包含源项而质量方程无源项。本文中研究了二维情况下基于 MUSCL 的伪弧长算法模型建立过程,针对网格自适应移动造成的物理空间扭曲变形而不易构造高精度格式的难题,基于坐标变换的思想,将变形的物理空间映射至一套正交的弧长计算空间,从而在弧长计算空间中实现控制方程的求解与新网格守恒变量的插值过程,提高了计算效率。编写了该算法相应的一维、二维程序,在人为解方法基础上验证其精度,借助某些经典算例,将PALM 数值结果
15、对比分析于有限体积法。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(
16、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 为守恒变量
17、的内外重构值;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)kKi,j+Ki,jS(?wi,j)(6)式中:kKi,j为网格 Ki,j的第 k 条边,k=1,2,3,4;kKi,j为第 k 条边的长度;nki,j为边界kKi,j的单位外法向量;wint(k)i
18、,j,wext(k)i,j为网格在边界的内外逼近值,如图 1 所示,且当wint(k)i,j=?wi,j,wext(k)i,j=?wi,j时,上式(6)退化为一阶格式。现对网格单元边界进行二阶重构,具体 MUSCL 形式11 如下:wint(1)i,j=?wi,j+12(?wi,j?wi,j+1,?wi,j1?wi,j)wext(1)i,j=?wi,j1+12(?wi,j1?wi,j,?wi,j2?wi,j1)wint(2)i,j=?wi,j+12(?wi,j?wi1,j,?wi+1,j?wi,j)wext(2)i,j=?wi+1,j+12(?wi+1,j wi,j,?wi+2,j?wi+1,
19、j)wint(3)i,j=?wi,j+12(?wi,j?wi,j1,?wi,j+1?wi,j)wext(3)i,j=?wi,j+1+12(?wi,j+1?wi,j,?wi,j+2?wi,j+1)wint(4)i,j=?wi,j+12(?wi,j?wi+1,j,?wi1,j?wi,j)wext(4)i,j=?wi1,j+12(?wi1,j?wi,j,?wi2,j?wi1,j)(7)式中:为非线性限制器函数。在计算中 取4:(a,b)=(sign(a)+sign(b)|ab|a|+|b|+(8)式中:为小量,可取 =109。对于时间的离散,忽略掉网格索引 i,j,采取如下的三阶TVD-K 格式:?
20、w(1)=?w(z)+tL(?w(z)?w(2)=34?w(z)+14?w(1)+14tL(?w(1)?w(z+1)=13?w(z)+23?w(2)+23tL(?w(2)L(w)=(4k=1kh(wint(k),wext(k),nk)ds+|K|S(w)|K|(9)式中:t 为时间步长;?w(z)为上步物理量,?w(z+1)为更新的物理量;L(?w(k)为微分算子,?w(1)、?w(2)为子步插值。图 1网格单元计算示意图Fig 1 Schematic diagram of grid cell calculation1 2网格自适应移动与计算空间的转换以 x=(x,y)和 =(,)分别代表物理
21、空间坐标和弧长空间坐标,(xi 1,j 1,xi 1,j,xi,j,xi,j 1)为网格单元 Ki,j的4 个顶点,从计算空间 c到物理空间 p的一一对应坐标映射为(x,y)=(x(,),y(,)。多维空间要同时照顾到网格的移动速度、方向和尺寸问题,且期望网格朝物理量梯度大之处进行移动,借助变分原理可知,网格的自适应移动将受泛函 E(x,y)最小值影响12。E(x,y)=12c(TxG1x+TyG2y)dd(10)式中:G1,G2为给定的正定对称矩阵;=(,)T,且有=/,=/,则式(10)对 应 的 Euler-Lagrange 方程13 为:(G1x)+(G1x)=0(G2y)+(G2y)
22、=0(11)它对应于泛函临界点。定义伪弧长控制函数 =1+a1|w|2+a2|w|2,这里采用单一物理量(其可为,p,化学反应度等)及其梯度相结合的模型,其中 a1,a2为可调参数。借助低通滤波器进行光滑打磨,二维情况表达式为:i,j14i,j+18(i+1,j+i,j+1+i1,j+i,j1)+116(i+1,j+1+i+1,j1+i1,j+1+i1,j1)(12)07兵 器 装 备 工 程 学 报http:/bzxb cqut edu cn/具体光滑次数视情况而定,这里不再赘述。令 G=wI,如此,进一步简化式(11),有:(x)=0(y)=0(13)对式(13)采用 Gauss-Seid
23、el 迭代法进行计算,具体步骤如下:xz+1i,j=i+1/2,jxzi+1,j+i1/2,jxz+1i1,j+i,j+1/2xzi,j+1+i,j1/2xz+1i,j1i+1/2,j+i1/2,j+i,j+1/2+i,j1/2i1/2,j=12(w(z)i,j)+(w(z)i,j1)i+1/2,j=12(w(z)i+1,j)+(w(z)i+1,j1)i,j1/2=12(w(z)i,j)+(w(z)i1,j)i,j+1/2=12(w(z)i,j+1)+(w(z)i1,j)(14)式中:xz+1i,j为下一时刻网格的坐标,其中,N,N代表弧长空间,方向的网格总数,1iN,1jN,且对于边界(i=
24、0,j=0,i=N,j=M),网格保持不动。由式(14)知网格移动速度表达式:x=x(z+1)i,j x(z)i,j=x(z+1)i,j x(z)i,j t(15)网格移动完之后,需要得到物理量在新网格上的值,这一步称作物理量重映,本文中采用基于通量的重映方法。如图 2 二维新旧网格转化图所示。图 2二维新旧网格转换示意图Fig 2 Schematic diagram of old grid transformingto new grid in 2D设 Dk为经过一次迭代物理空间中网格变化导致边kKi,j所扫过的面积,则有:Dk=124n=1(xnyn+1 xn+1yn)(16)基于前人的工作
25、14,考虑到新网格 xk+1 与旧网格 xk 之间通量守恒,得到如下守恒插值格式:K(z+1)i,j?w(z+1)i,j=K(z)i,j?w(z)i,j+4k=1Fk(wint(k)i,j,wext(k)i,j)(17)式中:K(z+1)i,j代表新网格的面积;w(z+1)i,j代表新网格的物理量;wint(k)i,j,wext(k)i,j的计算方法见式(7),守恒插值算法的核心在于 wint(k)i,j,wext(k)i,j的精确构造,考虑到变形物理空间中不易构造高精度格式,因此将重构过程转换至均匀正交的弧长计算空间(详见图 2(b)中。为书写方便而不失一般性,此处略掉上下角标,Fk(,),
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 基于 坐标 变换 间断 问题 伪弧长 算法 陈泽平
1、咨信平台为文档C2C交易模式,即用户上传的文档直接被用户下载,收益归上传人(含作者)所有;本站仅是提供信息存储空间和展示预览,仅对用户上传内容的表现方式做保护处理,对上载内容不做任何修改或编辑。所展示的作品文档包括内容和图片全部来源于网络用户和作者上传投稿,我们不确定上传用户享有完全著作权,根据《信息网络传播权保护条例》,如果侵犯了您的版权、权益或隐私,请联系我们,核实后会尽快下架及时删除,并可随时和客服了解处理情况,尊重保护知识产权我们共同努力。
2、文档的总页数、文档格式和文档大小以系统显示为准(内容中显示的页数不一定正确),网站客服只以系统显示的页数、文件格式、文档大小作为仲裁依据,平台无法对文档的真实性、完整性、权威性、准确性、专业性及其观点立场做任何保证或承诺,下载前须认真查看,确认无误后再购买,务必慎重购买;若有违法违纪将进行移交司法处理,若涉侵权平台将进行基本处罚并下架。
3、本站所有内容均由用户上传,付费前请自行鉴别,如您付费,意味着您已接受本站规则且自行承担风险,本站不进行额外附加服务,虚拟产品一经售出概不退款(未进行购买下载可退充值款),文档一经付费(服务费)、不意味着购买了该文档的版权,仅供个人/单位学习、研究之用,不得用于商业用途,未经授权,严禁复制、发行、汇编、翻译或者网络传播等,侵权必究。
4、如你看到网页展示的文档有www.zixin.com.cn水印,是因预览和防盗链等技术需要对页面进行转换压缩成图而已,我们并不对上传的文档进行任何编辑或修改,文档下载后都不会有水印标识(原文档上传前个别存留的除外),下载后原文更清晰;试题试卷类文档,如果标题没有明确说明有答案则都视为没有答案,请知晓;PPT和DOC文档可被视为“模板”,允许上传人保留章节、目录结构的情况下删减部份的内容;PDF文档不管是原文档转换或图片扫描而得,本站不作要求视为允许,下载前自行私信或留言给上传者【自信****多点】。
5、本文档所展示的图片、画像、字体、音乐的版权可能需版权方额外授权,请谨慎使用;网站提供的党政主题相关内容(国旗、国徽、党徽--等)目的在于配合国家政策宣传,仅限个人学习分享使用,禁止用于任何广告和商用目的。
6、文档遇到问题,请及时私信或留言给本站上传会员【自信****多点】,需本站解决可联系【 微信客服】、【 QQ客服】,若有其他问题请点击或扫码反馈【 服务填表】;文档侵犯商业秘密、侵犯著作权、侵犯人身权等,请点击“【 版权申诉】”(推荐),意见反馈和侵权处理邮箱:1219186828@qq.com;也可以拔打客服电话:4008-655-100;投诉/维权电话:4009-655-100。