-
摘要: 为解决松软煤体爆破孔成形困难导致深孔预裂爆破技术无法应用的问题,选择在运输巷底板距煤 层较近的岩石中进行爆破,利用数值模拟方法进行了相关理论探讨。分析了煤岩介质和单煤体介质中布孔的 差异,建立了单煤体和煤岩体深孔预裂爆破的5个数值计算模型。研究了单煤体和煤岩介质爆破孔与抽放孔 连心线上有效应力随距离的变化,从岩孔爆破传播到煤层其应力衰减的程度较单煤层中大得多,但在靠近抽 放孔附近煤层,二者的差距变小。抽放孔轴线方向所受有效应力的大小是决定爆破效果的重要因素,煤岩介 质中爆破孔与抽放孔间距为2.0m 时,抽放孔轴线方向上的平均有效应力与单煤层中爆破孔与抽放孔间距 为3.0m 时的相当,可作为实现岩孔爆破效果的布孔参数。Abstract: Toovercomethedifficultyinusingthedeep-holepre-splittingblastingtechnologyinthesoft coalseamspronetoresultingintheboreholecollapse,anewmethodwasappliedtoincreasingthegas drainagerate,inwhichtherocknearbythecoalseamwasexplodedatthebottomofthelaneway.The theoriesforthetechniquewereexploredbynumericalsimulation.Layoutsofblastholesanddrain holeswerecomparedthesinglecoalseam andthecoal-rock media.Fivenumericalcomputational modelswereconstructedbyDYNA3D,includingasinglecoalmodelandfourcoal-rockmodelswith thedifferentdistancebetweentheblastholeandthedrainhole.Thecorrespondingeffectivestress wascomputedalongthelinkinglineoftheblastholeandthedrainhole.Thecomputationalresults showthattheattenuationoftheeffectivestressfromtheblastholetothedrainholeismoreinthe coal-rockmediathaninthesinglecoal,whereastheeffectivestressdifferencebetweenthetwomedia decreasesinthecoalseamsadjacenttothedrainholes.Thedistributionoftheeffectivestressalong theaxisofthedrainholeissignificanttothedraineffect.Themeaneffectivestressalongtheaxisof thedrainholeincoal-rock media,whenthedistancebetweentheblastholeandthedrainholeis 2.0m,isappropriatelythesameasthatinthesinglecoalseam,whenthedistancebetweentheblast holeandthedrainholeis3.0m.Anditcanbeusedasthelayoutparametertokeepthesameblast effectinthecoal-rockmediaasinthesinglecoalseam.
-
Key words:
- mechanicsofexplosion /
- dynamicstress /
- computationalmodel /
- coalseam /
- rock /
- blasting
-
在包含运动边界的非定常流动模拟中,流场物理空间随时间变化,边界运动导致网格也随时间变化,典型应用包括多体分离、机翼抖振颤振、舱门启闭等[1]。CFD计算一般采用运动嵌套网格、笛卡尔自适应网格和变形动网格等技术处理网格变化。在采用上述技术处理网格变化时,均涉及流场信息的插值问题[1-2],即利用旧网格上的流动信息插值获得重构区域流动参数。此外,在爆轰流场模拟中,爆轰计算网格与冲击波传播计算网格差两个量级[3-4],为提高计算效率,往往需要将物理量由爆轰计算时使用的密网格插值传递到冲击波传播模拟中的疏网格。
插值过程一般应具备守恒性和单调性等基本特性。对爆轰计算等守恒性要求高的问题,一般的插值不能保证计算区域内物理量的守恒,会降低计算精度,且无法捕捉激波的正确位置[5],需要发展守恒插值方法。守恒插值方法主要可分为两类,基于面通量(SFB/DC,simplified face-based donor-cell)的方法和基于单元相交(CIB/DC, cell-intersection-based donor-cell)的方法[6]。
SFB/DC方法[6]效率较高,但要求新、旧网格具有相同的拓扑结构,不适用于网格重构等拓扑发生变化的情况。CIB/DC方法思想简单直接,对计算网格的类型、结构不进行任何假设,适用范围很广,并且具有保号性,插值过程中不会产生新的极值[6]。但该方法推广至三维情况时,由于三维网格单元重叠区域切割计算非常复杂,给这种方法的应用带来了很大困难。为了实现严格意义下的三维CIB/DC算法,P.E.Farrell等[7]进行了一系列卓有成效的工作,在2011年进一步发展为“局部超网格”(Local Supermesh)方法,成功实现了非结构网格下的三维CIB/DC算法[8]。
本文中在P.E.Farrell等和S.Menon等[9]工作的基础上,构造了应用于二维/三维混合网格的CIB/DC算法,实现了基于精确网格切割的积分守恒插值。验证结果表明,该方法能够保证插值过程中计算域内物理量的严格守恒,具有比常规二阶插值更高的精度。
1. 守恒插值基本思想
守恒插值要求在插值过程中保证守恒量的积分值保持不变,假设计算域为Ω,网格Ta是计算域的一种划分,则在网格Ta中,Ω被划分为Na个网格单元Ta=(T1a, T1a, …, TNaa),并满足如下条件:
Ω=Na⋃iTai,Tai∩Tbj=∅i≠j 在网格Tia中,物理量ϕ的分布函数为:
ϕa(x)=ϕai(x)ifx∈Tai 式中:ϕia(x)可以指定多种形式,通常为多项式分布。
网格Tb是计算域Ω的另一种形式的划分,划分形式与要求同Ta,则同样有Tb上的物理量分布ϕb(x)。如果对于Ω的任意子区域D,满足如下条件:
∫Dϕa(x)dV=∫Dϕb(x)dV∀D⊂Ω (1) 则称ϕb(x)是ϕa(x)的守恒插值结果。
对于一般的情形,要使式(1)满足,必须有Tb≡Ta,ϕb(x)≡ϕa(x),但这就失去了在两套网格间进行插值的意义。为此,对D进行限制,不能是Ω的任意子区域,而必须是网格Tb中的某个网格:
∫Dϕa(x)dV=∫Dϕb(x)dV∀D⊂{Tbi|i=1,2,⋯,Nb} (2) 上式就是一般意义上的守恒插值关系式。
对式(2),按照基于单元相交的思想,对D进行分解:
D=Tbi=Na∑j=1(Tbi∩Taj)=Na∑j=1Vij 式中:Vij=Tib∩Tja。将上式代入式(2),则有:
∫Dϕb(x)dV=∫Dϕa(x)dV=Na∑j=1∫Vijϕa(x)dV (3) 式(3)就是基于单元相交守恒插值方法(CIB/DC)的基本依据。从式(2)到式(3)未引入任何近似。因此,在给定分布ϕ(x)的情况下,CIB/DC方法是严格守恒的。
从式(3)还可以看出,如果已知网格Ta上的物理量分布ϕa(x),如何计算网格单元间的相交区域Vij便成为CIB/DC方法的关键。
2. “超网格”基本思想
图 1为文献[8]给出的超网格示意图,超网格是由Ta和Tb所有结点以及它们边界交点构成的网格,具有如下特性:
(1) 每一个超网格不再和Ta或Tb中任何网格单元相交,即它总可以在Ta中找到一个网格包含或等于它,也总可以在Tb中找到一个网格包含或等于它;
(2) Ta的每一个网格总可以明确表示为一个或多个超网格的集合,Tb的每一个网格总可以明确表示为一个或多个超网格的集合。
这样一来,对于Tb任意网格D可以写出它所包含的超网格,每一个超网格就是Vij=Tib∩Tja,可以在Ta上找到具体位置和对应的物理量,采用ϕb(x)≡ϕa(x),按照D所包含的超网格进行求和,就得到D的物理量,从而实现了从Ta到Tb的信息传递。由于Ta和Tb在超网格的交集完全等价,使传递的物理量严格守恒,实现精确守恒插值。
3. 单元相交算法
由上节可知,获得超网格是实现精确守恒插值的关键。下面给出计算平面凸多边形和三维凸多面体交集的具体算法。
3.1 平面凸多边形相交算法
如图 2所示,将目标三角形A1A2A3视为切割多边形,使用切割三角形B1B2B3对目标三角形进行切割,每次切割时只保留目标多边形中与切割多边形在切割线同侧的部分。图中依次使用△B1B2B3的3条边B1B2、B2B3、B3B1切割△A1A2A3,切割后保留的多边形依次为CA3A2G、CDEA2G、CDEFB1,多边形CDEFB1即为△A1A2A3和△B1B2B3的交集。
在使用直线去切割平面凸多边形时,以图 3为例,切割直线经过点O,其法向单位矢量为n,切割直线将平面划分为两部分,定义法向n所指的半平面为正半平面,反方向的半平面为负半平面。如果切割之后目标多边形仅保留正半平面的部分,那么切割操作等效于求目标多边形和正半平面的交集。定义平面上任意一点P到切割直线的距离为:
d=xOP⋅n=(xP−xO)⋅n 式中:xP和xP分别为点P和点O坐标。
根据上述定义,正半平面的点到切割直线的距离均为正数,负半平面的点到切割直线的距离均为负数,切割直线上的点到切割直线的距离为零。如果凸多边形各顶点到切割直线的距离均为非负数(图 3(a)),则切割结果即为原多边形;如果凸多边形各顶点到切割直线的距离均为非正数,则切割结果为空集;除上述两种情况外,凸多边形与切割直线必然有两个交点(图 3(b))。具体如下:
(1) 设目标多边形为N,其顶点依次为N1, N2, …, NN;切割直线为l,切割结果为R;
(2) 计算N的各个顶点到l的距离,记顶点Ni到l的距离为di;
(3) 初始化一个空的有序点集P;
(4) 从顶点N1起,依次对多边形N的边N1N2, N2N3, …, NNN1进行判定。
对线段NiNj的判定过程为:如果didj<0,则线段与切割直线有新的交点,记交点为Nnew;依次将Ni、Nnew、Nj中距离为非负数的顶点加入点集P;如果didj≥0,则依次将Ni、Nj中距离为非负数的顶点加入点集P;
(5) 完成所有边的判定后,如果P为空集,则R为空集;若P为非空集,则依次连接P中的点形成的多边形就是所求的切割结果R。
3.2 三维凸多面体相交算法
在三维流体计算中,常用的网格包括四面体、四棱锥、三棱柱、六面体等多种类型。其中,四棱锥、三棱柱、六面体这3种网格都具有四边形表面,网格生成过程中不能严格保证四边形表面的四个顶点位于同一平面内,会给多面体表面的定义带来歧义。
为了简化算法的复杂性,同时消除歧义,在计算两个任意凸多面体的交集时,首先将凸多面体分解为多个小四面体,如图 4所示,其中点N6为点N2、N3、N4、N5的重心;然后通过计算这些小四面体的交集,来获得多面体的交集。计算多面体Pa和Pb的交集的具体步骤为:
(1) 将Pa分解为一组四面体{Tai},将Pb分解为一组四面体{Tbj},并且满足:
Na⋃i=1Tia=Pa,Tia∩Tja=∅∀1≤i,j≤Na且i≠jNb⋃j=1Tjb=Pb,Tib∩Tjb=∅∀1≤i,j≤Nb且i≠j (2) 计算{Tai}和{Tbi}中任意两个四面体的交集:
Tijc=Tia∩Tjb1≤i≤Na,1≤j≤Nb 并记录交集Tcij的体积和重心分别为Vcij,xcij。
(3) 多面体Pa和Pb交集可表示为:
Pa∩Pb=(Na⋃i=1Tia)∩(Nb⋃j=1Tjb)=Na⋃i=1Nb⋃j=1(Tia∩Tjb)=Na⋃i=1Nb⋃j=1Tjb∩Tijc 交集的体积和重心分别是:
Vc=∑i,jVijc,xc=1Vc∑i,jVijcxijc 该方法的核心在于如何计算两个四面体的交集Tcij。相对于平面多边形的情形,即使四面体是最简单的多面体,其交集的计算也是十分复杂的:两个四面体相交,最多可能产生一个12顶点的8面体。本文中使用切割法计算两个四面体的交集:将一个四面体看做目标四面体,另一个四面体为切割四面体,依次使用切割四面体的4个表面去切割目标四面体,最终留下的多面体就是这两个四面体的交集。
4. 相交单元的搜索算法
相对于通常的流场插值方法,守恒插值的计算效率较低,主要有两个原因:(1)计算两个单元交集的算法较为复杂;(2)对同一个新网格单元,普通插值方法只需用到一个旧网格单元(称为贡献单元)的信息,而守恒插值需要用到多个贡献单元的信息。前文已对单元相交算法进行了详细的描述,本节将介绍守恒插值中多个贡献单元的搜索算法。
在贡献单元的搜索算法中,首先需要快速定位出新网格单元在旧网格中的位置。对此,可采用四叉树[11]或八叉树[12](三维)方法加速搜寻,如图 5(a)所示,完全填充网格为旧网格,其上的非填充网格为新网格中的某个单元,通过四叉树方法,可快速定位出新网格单元的格心所在的旧网格单元(图中填充单元),这是新网格单元的第1个贡献单元;然后,逐个搜索各个贡献单元的直接相邻单元,如果相邻单元与新网格单元相交,则将该单元加入贡献单元的列表,如图 5(b)中浅色填充单元所示;当所有贡献单元的相邻单元(图 5(c)与(b)中不同的填充单元)都不与新网格单元相交时,搜索结束。
5. 守恒插值验证算例
设计了3个算例来验证本文发展的守恒插值方法。第1个算例主要考核相同拓扑结构下网格点位置变化的影响;第2个算例主要考核网格形状的影响;第3个算例主要考核三维情况下的应用情况。
算例 1
考虑如下的测试函数:
\phi (x, y) = 1 + \sin (2{\rm{ \mathsf{ π} }}x)\sin (2{\rm{ \mathsf{ π} }}y) 计算区域为[0, 1]×[0, 1],计算网格为一组拓扑结构相同的非结构三角形网格,分别表示为M0、M1、…、MN,其中M0如图 6所示。网格MN中的第(i, j)个网格点的坐标由下式给出:
\left\{ {\begin{array}{*{20}{l}} {x(\xi , \eta , t) = (1 - \alpha (t))\xi + \alpha (t){\xi ^2}}\\ {y(\xi , \eta , t) = (1 - \alpha (t))\eta + \alpha (t){\eta ^2}} \end{array}} \right. 式中:α(t)=0.5sin(4πt), ξ=(i-1)/(imax-1), η=(j-1)/(jmax-1)t=n/T, 其中i=1, 2, …, imax; j=1, 2, …, jmax; n=0, 1, …, N。
测试过程中,在网格M0上给定测试函数的准确分布,然后依次插值到网格M1、M2、…、MN上,共进行N次插值,最后比较网格MN上的函数分布与准确值之间的误差。在本文计算中,取imax=jmax=65,N=200,T=20。
图 7给出了计算误差随插值次数的变化情况,其中纵坐标代表计算区域内各单元函数分布插值后计算值与理论值之差的1范数。两种插值方法中,插值误差都随插值次数的增加而增大,但守恒型插值方法的误差始终较小,不到常规二阶插值误差的一半。
图 8为两种方法的物理量积分值比较,其中纵坐标代表每个单元计算值总和与理论值的比,从图中明显可以看出,守恒插值保证了计算域内物理量的守恒。
算例2
考虑如下的测试函数:
\phi (r) = 2 + \cos ({\rm{ \mathsf{ π} }}r/L) 计算区域是边长为1的正方形,L为正方形的对角线长度,r是与正方形中心的距离。测试中使用两种不同的网格,网格1是均匀分布的结构网格,网格2是非结构网格,并与网格1具有相同的边界网格点分布,如图 9所示。测试中,首先在网格1上给出测试函数的准确分布,然后将其插值到网格2上,再插值回网格1,如此反复,进行200次插值后,将物理量分布与初始的准确分布进行比较,统计插值过程中的误差。分别采用常规二阶插值和守恒型二阶插值进行测试,并使用不同疏密的网格以分析插值方法的精度。测试中使用的4组网格的网格量如表 1所示。
表 1 计算网格参数Table 1. Parameters of the computational grid序号 结构网格 非结构网格 节点数 单元数 1 33×33 695 1 260 2 65×65 2 593 4 928 3 129×129 10 044 19 574 4 257×257 39 687 78 348 图 10是计算误差与网格步长之间的关系,其中横坐标是表 1中4个不同结构网格计算单元的边长,纵坐标代表计算区域内各单元函数分布插值后计算值与理论值之差的1范数。可以看出在相同的网格步长下,守恒型插值的误差更小,而且误差与网格步长之间的变化斜率更大,说明其具有更高阶的精度。计算结果表明,常规二阶插值精度约为1.4阶,而守恒型二阶插值精度约为1.87阶。表 2是不同算例的守恒误差比较,守恒插值方法的守恒误差接近机器零,明显优于普通二阶插值。
表 2 不同算例的守恒误差Table 2. Conservation errors indifferent numerical cases序号 二阶插值 守恒型二阶插值 1 3.724×10-3 0 2 -1.447×10-3 0 3 -3.985×10-4 1×10-14 4 5.636×10-5 -2×10-14 算例3
考虑如下的测试函数:
\varphi (x, y, z) = 1 + \sin (2{\rm{ \mathsf{ π} }}x)\sin (2{\rm{ \mathsf{ π} }}y)\sin (2{\rm{ \mathsf{ π} }}z) 计算区域是边长为1的正方体。每次测试使用两套网格,首先在网格1上给出测试函数的准确分布,然后将其插值到网格2上,再插值回网格1,如此反复,进行10次插值后,将物理量分布与初始的准确分布进行比较,统计插值过程中的误差。分别采用普通二阶插值和守恒型二阶插值进行测试,并使用不同疏密的网格以分析插值方法的精度。测试中使用的4组网格的网格量如表 3所示,其中网格步长由下式得到:
{\rm{d}}x = {\left( {\sqrt {{N_{E1}}{N_{E2}}} } \right)^{1/{N_D}}} 表 3 计算网格参数Table 3. Parameters of the computational grid序号 网格1 网格2 dx 节点数 单元数 节点数 单元数 1 914 4 033 676 2 711 0.067 2 2 222 10 632 1 424 6 075 0.050 3 6 527 32 913 3 733 16 681 0.035 4 26 099 138 542 12 600 59 347 0.022 式中:NE1和NE2分别为两套网格的单元个数,ND=3为网格的维数。
图 11是测试函数在三维空间中的分布云图,图 12是计算误差与网格步长之间的关系,其中横坐标是表 3中dx的值,纵坐标代表计算区域内各单元函数分布插值后计算值与理论值之差的1范数。在三维条件下,守恒型插值依然具有更小的计算误差和更高的计算精度。
6. 结论
采用在新旧网格相交单元剖分构建出来的超级网格作为网格之间信息传递的基础,采用四叉树(二维)或八叉树(三维)方法搜寻算法提高计算效率,基于以上方法提出了网格之间流场物理量的守恒插值算法。该算法理论上可以实现网格之间物理量的精确守恒,二维和三维算例也表明产生的误差接近机器零,明显优于常用的二阶插值。本文中提出的方法可以解决爆轰流场模拟中网格尺度和物理量相差较大情况下的插值难题。
期刊类型引用(7)
1. 余昕,吴乘胜,金奕星. 面向交界面网格的超网格算法研发及CFD应用. 中国造船. 2023(03): 233-246 . 百度学术
2. 王子维,刘东健,陈逖. 单级跨声速压气机旋转失速模拟方法. 气体物理. 2021(01): 30-37 . 百度学术
3. 常思源,白晓征,崔小强,刘君. 一种改进的非定常激波装配算法. 航空学报. 2020(02): 163-176 . 百度学术
4. 刘耘臻,万志强,杨超. 飞行载荷外部气动力的二次规划等效映射方法. 北京航空航天大学学报. 2020(03): 541-547 . 百度学术
5. 孙伟,陈泽栋. 风载荷对飞行器尾罩分离影响的数值模拟分析. 弹箭与制导学报. 2019(06): 102-106 . 百度学术
6. 崔鹏程,唐静,李彬,马明生,邓有奇. 基于超网格的重叠网格守恒插值方法. 航空学报. 2018(03): 23-35 . 百度学术
7. 靳晨晖,王刚,王泽汉. 子母弹多体分离过程的非定常CFD/RBD数值仿真. 气体物理. 2018(04): 47-63 . 百度学术
其他类型引用(3)
-
计量
- 文章访问数: 2386
- HTML全文浏览量: 80
- PDF下载量: 225
- 被引次数: 10