一种针对任意构型爆源地运动稳态效应的计算方法与流程

文档序号:31053311发布日期:2022-08-06 09:04阅读:390来源:国知局
一种针对任意构型爆源地运动稳态效应的计算方法与流程

1.本发明涉及工程爆破效应数值计算领域,具体涉及一种针对任意构型爆源地运动稳态效应的计算方法。


背景技术:

2.工程爆破在基建、采矿等国民经济活动中应用广泛,各种类型和构型的爆源爆炸广泛应用于经济建设、国防和科研领域。准确计算和评估各种类型和构型的爆源在地介质中爆炸产生的地运动效应,对爆炸效应评估、爆破方案设计以及安全防护有重要意义。
3.以往的工程爆破地运动稳态效应计算方法只针对球形爆源地运动稳态效应,给出了一维的计算方法,但是其计算范围大、计算时间长,同时并不适用于非球型爆源爆炸产生的地运动稳态效应。因此,研究设计一种针对任意构型爆源爆炸产生的地运动稳态效应,且高效简单的计算方法是非常必要的。


技术实现要素:

4.本发明的目的是解决现有的对称爆源爆炸产生的地运动稳态效应的计算,采用的一维数值计算方法、特征取值点的选取以及地运动稳态效应计算,计算数据冗余、耗时长;且不适用于其他非球对称爆源爆炸产生的地运动稳态效应的问题,而提供了一种针对任意构型爆源地运动稳态效应的计算方法。
5.为达到上述目的,本发明采用的技术方案为:
6.一种针对任意构型爆源地运动稳态效应的计算方法,其特殊之处在于,包括以下步骤:
7.s1、依据实验条件或理论设计,获取爆源的特征以及周围地介质的类型;
8.所述爆源的特征包括爆源的形状、特征尺寸、类型、爆炸当量、整体位置和起爆点;
9.s2、根据爆源的特征和周围地介质类型,建立物理模型,并确定物理模型的最小计算范围和最小计算时间;
10.s3、简化物理模型,建立有限元计算模型,进行有限元网格划分,并选择合理的材料模型、材料参数及流固耦合算法;
11.s4、根据物理模型的最大塑性区半径,确定系列取值点的爆心距位置,确定取值点的位置及取值点个数;s5、根据s2中确定的最小计算时间,采用动力有限元程序进行计算,输出各个取值点的物理量的时间历程曲线;
12.所述物理量包括取值点的等效应力、等效应变或位移;
13.s6、对s5所输出系列取值点的物理量的时间历程曲线进行数据处理,得到爆源爆炸产生的地运动稳态效应值。
14.进一步地,步骤s2中:
15.所述物理模型的最小计算范围为:
[0016][0017]
式中,l为地介质的计算范围爆心距;w为爆源的爆炸当量;c
p
为地介质的纵波波速;
[0018]
所述物理模型的最小计算时间为:
[0019][0020]
式中,t为模型的计算时间。
[0021]
进一步地,步骤s3中,所述简化物理模型是根据爆源的几何特征,具体为:
[0022]
若爆源为轴对称,则简化为二维轴对称模型;
[0023]
若爆源有3个对称面,则简化为1/8对称模型;
[0024]
若爆源有2个对称面,则简化为1/4对称模型;
[0025]
若爆源仅有1个对称面,则简化为1/2对称模型;
[0026]
若爆源无对称性,则建立完整的三维模型。
[0027]
进一步地,步骤s4具体为:
[0028]
s41、利用动力有限元程序进行试算,根据计算结果确定物理模型的最大塑性区半径r
plastic

[0029]
s42:确定系列取值点的爆心距位置:
[0030]rpoints
=1.5r
plastic
ꢀꢀ
(3)
[0031]
式中,r
points
是系列取值点的爆心距位置,r
plastic
是s41中试算得到的最大塑性区半径;
[0032]
s43:根据s3中建立的物理模型特征以及s42中确定的系列取值点的爆心距大小,确定取值点的位置及取值点个数n。
[0033]
进一步地,若s3中简化的物理模型为二维轴对称,则取值点个数n=7,取值点位于爆心距为r
points
的环线上且从0
°
到90
°
每隔15
°
的等间隔角位置处;
[0034]
若s3中简化的物理模型为1/8对称,则取值点个数n=18,取值点分别位于计算模型的3个对称面上、爆心距为r
points
的环线上且每隔15
°
的等间隔角位置处;
[0035]
若s3中简化的物理模型模型为1/4对称,则取值点个数n=29,取值点分别位于计算模型的2个对称面和另外一个坐标面上、爆心距为r
points
的环线上且每隔15
°
的等间隔角位置处;
[0036]
若s3中简化的物理模型为1/2对称,则取值点个数n=45,取值点分别位于计算模型的1个对称面和另外两个坐标面上、爆心距为r
points
的环线上且每隔15
°
的等间隔角位置处;
[0037]
若s3中爆源无对称性,则建立完整的三维模型,则整体三维模型取值点个数n=66,取值点分别位于计算模型的3个对称面上、爆心距为r
points
的环线上且每隔15
°
的等间隔角位置处。
[0038]
进一步地,步骤s6具体为:
[0039]
s61、对每个取值点的物理量时间历程曲线进行规整化处理,对数据进行插值,按
等间距时间间隔δt输出物理量的时间历程曲线;
[0040]
s62:准确获得每个取值点的物理量的稳态值:
[0041]
读取各个取值点的物理量的时间历程曲线,在曲线大致稳定后在至少3个振动的时间周期范围内对物理量取算术平均值,将该值作为该取值点上的稳态值:
[0042][0043]
式中,pj为j取值点上的物理量的稳态值;pi为i时刻j取值点上物理量值;p
i+1
为i+1时刻j取值点上物理量值;δt为时间历程曲线数据的输出时间间隔;ts为取值的起始时刻;te为取值的终止时刻;
[0044]
s63、将所有取值点的物理量稳态值进行算术平均,得到该构型爆源爆炸产生的地运动稳态效应值
[0045][0046]
式中,n为s43中确定的取值点个数。
[0047]
进一步地,所述流固耦合算法采用带侵蚀算法的罚函数耦合算法,耦合单元内的积分点个数为3,并在所有方向上都进行耦合;
[0048]
所述流固耦合算法中,炸药和空气采用欧拉网格,使用流体输运算法,具体采用van leer算法和半指数漂移二阶精度算法;地介质采用拉格朗日算法。
[0049]
进一步地,若爆源的类型为tnt炸药,对tnt采用高能炸药燃烧材料模型和jwl状态方程描述炸药爆轰过程,其中压力p和内能e及其相对体积v满足:
[0050][0051]
其中,a、b、r1、r2和ω为jwl状态方程参数,e为炸药的体积比内能。
[0052]
进一步地,对空气采用空气材料模型,采用多项式状态方程描述:
[0053]
p=c0+c1μ+c2μ2+c3μ3+(c4+c5μ+c6μ2)e
ꢀꢀ
(7)
[0054]
其中,ρ0为空气的初始密度,ρ为空气的当前密度。
[0055]
进一步地,若地介质为花岗岩,采用塑性硬化材料模型描述,塑性硬化材料模型的屈服应力满足:
[0056][0057]
其中,表示强度的应变率效应;
[0058]
σ0为初始屈服强度;
[0059]
β为控制运动硬化和各向同性硬化的参数;
[0060]ep
为塑性硬化模量,
[0061]
为有效塑性应变。
[0062]
采用上述方法,可有效表征出爆源的任意构型产生的地运动稳态效应的差异,且计算范围小、计算时间短,成量级地提高了计算效率,可直接应用于实际工程问题中。
[0063]
与现有技术相比,本发明具有的有益技术效果如下:
[0064]
1、本发明提出的针对任意构型爆源地运动稳态效应的计算方法,不局限与目前工程爆破地运动稳态效应的计算方法仅适用于球形爆源的一维计算方法,本发明的计算方法适用于任意构型的爆源爆炸产生的地运动稳态效应,并且计算方法拓展到了二维或三维,计算模型可适用于地下爆炸、地面爆炸和空中爆炸,该方法使用范围更广,计算更为精确。
[0065]
2、本发明提出的针对任意构型爆源地运动稳态效应的计算方法,选取弹性区一定范围内的物理量的稳态值作为爆炸产生的地运动效应评估的参考量,有效降低了数值计算带来的多种计算误差。
[0066]
3、本发明提出的针对任意构型爆源地运动稳态效应的计算方法,其中系列取值点的确定,合理、有效地表征了爆源的任意构型的结构特征。
[0067]
4、本发明提出的针对任意构型爆源地运动稳态效应的计算方法,计算模型的计算范围小、计算时间短,成量级地提高了计算效率,更适用于实际工程问题的求解。
附图说明
[0068]
图1为本发明针对任意构型爆源地运动稳态效应的计算方法实施例设计流程图;
[0069]
图2为本实施例中的一种爆源的构型和尺寸示意图;
[0070]
图3为本实施例中的三维1/8对称计算模型示意图;
[0071]
图4为本实施例根据三维1/8对称爆源的塑性区最小计算范围确定取值点爆心距及系列取值点位置示意图;
[0072]
图5为本实施例中采用的对物理量的时间历程曲线取稳态值的数据取值区间示意图。
具体实施方式
[0073]
为使本发明的目的、优点和特征更加清楚,以下结合附图和具体实施例对本发明提出的一种针对任意构型爆源地运动稳态效应的计算方法作进一步详细说明。本领域技术人员应当理解的是,这些实施方式仅仅用来解释本发明的技术原理,目的并不是用来限制本发明的保护范围。
[0074]
如图1所示,本实施例提供了1kg tnt非球型装药在花岗岩地介质中爆炸产生地运动稳态效应的计算方法,具体包括以下步骤:
[0075]
s1、获取实验或理论设计中爆源形状及特征尺寸、爆源的类型(诸如tnt炸药或各种类型的岩石乳化炸药等)、爆炸当量、爆源的整体位置以及起爆点的位置等条件,以及获
取爆源周围的地介质类型(诸如土壤、砂石或不同类型的岩石介质等)及特征尺寸等参数。
[0076]
本实施例中,爆源的形状及特征尺寸,如图2所示,爆源的类型为tnt炸药,当量为1kg tnt,爆源起爆点位于炸药的几何中心;
[0077]
爆源周围的地介质为花岗岩c
p
=4250m/s;c
p
为地介质的纵波波速。
[0078]
s2、根据爆源的特征及周围地介质类型,确定物理模型的最小计算范围和最小计算时间;
[0079]
s2.1、根据爆源的特征以及爆源周围地介质类型,确定物理模型的最小计算范围,如式(1)所示,
[0080][0081]
式中,l为地介质的计算范围爆心距(单位为m),w为装药的tnt当量(单位为kg),c
p
为地介质的纵波波速(单位为m/s);
[0082]
本实施例中,当量为1kg tnt,花岗岩,根据式(1),确定模型的最小计算爆炸范围为8m;
[0083]
s22、根据爆源及药室的特征以及爆源周围地介质类型,确定物理模型的最小计算时间,如式(2)所示,
[0084][0085]
式中,t为模型的计算时间(单位为s);
[0086]
本实施例中,确定模型的计算时间为3ms。
[0087]
s3、合理简化物理模型并建立有限元计算模型;
[0088]
s31、根据爆源的几何特征,合理简化物理模型;
[0089]
简化原则为:若爆源的几何特征为轴对称,则简化为二维轴对称模型;若爆源有3个对称面,则简化为1/8对称模型;若爆源有2个对称面,则简化为1/4对称模型;若爆源仅有1个对称面,则简化为1/2对称模型;若爆源无对称性,则建立完整的三维模型;
[0090]
本实施例中,根据图2所示爆源的几何特征,模型在三维直角坐标系中有3个对称面,将计算模型简化为1/8对称模型,如图3所示;
[0091]
s32、对s31中建立的物理模型,采用有限元前处理软件建立几何模型;
[0092]
s33、对s32中建立的几何模型,进行合理的有限元网格划分;
[0093]
s34、采用合适材料模型、材料参数及流固耦合算法。
[0094]
本实施例中,对tnt采用高能炸药燃烧材料模型(*mat_high_explosive_burn)和jwl状态方程描述炸药爆轰过程,其中压力p和内能e及其相对体积v的关系,如式(6)所示,
[0095][0096]
其中,a、b、r1、r2和ω为jwl状态方程参数,e为炸药的体积比内能。采用的tnt的模型参数,如表1所示。
[0097]
表1 tnt炸药材料参数(单位:cm-g-microsecond)
[0098][0099]
本实施例中,对空气采用空气材料模型(*mat_null),采用多项式状态方程描述,如式(7)所示,
[0100]
p=c0+c1μ+c2μ2+c3μ3+(c4+c5μ+c6μ2)e
ꢀꢀ
(7)
[0101]
其中,ρ0为空气的初始密度,ρ为空气的当前密度。采用的空气模型参数,如表2所示。
[0102]
表2空气参数(单位:cm-g-microsecond)
[0103][0104]
本实施例中,对花岗岩地介质采用塑性硬化材料模型描述。塑性硬化材料模型的屈服应力表述如式(8),
[0105][0106]
其中,表示强度的应变率效应,σ0为初始屈服强度,β为控制运动硬化和各向同性硬化的参数,e
p
为塑性硬化模量,其表达式为为塑性硬化模量,其表达式为为有效塑性应变。岩石的模型参数,如表3所示。
[0107]
表3岩石参数(单位:cm-g-microsecond)
[0108][0109]
采用流固耦合算法进行计算,合理设置流体输运算法及流固耦合计算算法参数。tnt和空气设置为流体输运算法,具体采用van leer算法和半指数漂移二阶精度算法;地介质采用拉格朗日算法;欧拉网格和拉格朗日网格采用流固耦合算法;流固耦合计算的算法采用耦合单元内的积分点个数为3,并在所有方向上都进行耦合;采用带侵蚀算法的罚函数耦合算法。
[0110]
s4、选取反映爆源特征的系列取值点;
[0111]
s41、利用动力有限元程序进行试算,根据计算结果确定模型的最大塑性区半径rplastic

[0112]
本实施例中模型的最大塑性区半径r
plastic
为87.35cm;
[0113]
s42、确定系列取值点的爆心距位置,包括:根据s41中的计算结果,确定系列取值点的爆心距位置r
points
,如式(3)所示,
[0114]rpoints
=1.5r
plastic
ꢀꢀ
(3)
[0115]
确定系列取值点的爆心距位置r
points
=131cm。
[0116]
s43、根据s31中建立的物理模型特征以及s42中确定的系列取值点的爆心距大小,确定取值点的位置及取值点个数n;
[0117]
具体为:若s31中的模型为二维轴对称,则取值点个数n=7,取值点位于爆心距为r
points
的环线上且从0
°
到90
°
每隔15
°
的等间隔角位置处;若s31中的模型为1/8对称,则取值点个数n=18,取值点分别位于计算模型的3个对称面上、爆心距为r
points
的环线上且每隔15
°
的等间隔角位置处;若s31中的模型为1/4对称,则取值点个数n=29,取值点分别位于计算模型的2个对称面和另外一个坐标面上、爆心距为r
points
的环线上且每隔15
°
的等间隔角位置处;若s31中的模型为1/2对称,则取值点个数n=45,取值点分别位于计算模型的1个对称面和另外两个坐标面上、爆心距为r
points
的环线上且每隔15
°
的等间隔角位置处;若s31中建立的模型为无对称性,则建立整体模型,取值点个数n=66,取值点分别位于计算模型的3个对称面上、爆心距为r
points
的环线上且每隔15
°
的等间隔角位置处。不同模型的取值点个数及位置如表4中所示。
[0118]
表4不同模型的取值点个数及位置
[0119][0120][0121]
本实施例中,s31中建立的简化模型为1/8对称模型,分别在计算模型的3个对称面上在等取值点爆心距131cm的环向位置处每隔15
°
等间隔角共选取18个取值点(n=18),取值点位置,如图4所示。
[0122]
s5、根据s22中确定计算时间,采用动力学有限元程序进行计算,并输出每个取值点的物理量(如等效应力、等效应变或总位移等)的时间历程曲线。
[0123]
本实施例中,采用动力学有限元程序计算,计算总时间为3ms。本实施例中选取等效应力作为参考物理量,在计算过程中输出每个取值点上的等效应力时间历程曲线。
[0124]
s6、对所输出系列取值点的相关物理量的时间历程进行数据后处理,得到爆源爆炸产生的地运动稳态效应结果。
[0125]
s61、对每个取值点的物理量的时间历程曲线进行规整化处理,对数据进行插值,按等间距时间间隔δt输出物理量的时间历程;
[0126]
本实施例中,按等间距时间间隔δt=5
×
10-3
ms输出每个取值点的等效应力的时间历程;
[0127]
s62、准确获得每个取值点的物理量的稳态值。读取各个取值点的物理量的时间历程曲线,在曲线大致稳定后在至少3个振动的时间周期范围内对物理量取算术平均值,将该值作为该取值点上的稳态值,计算方法如式(4)所示,
[0128][0129]
式中,pj为j取值点上的物理量的稳态值,pi为i时刻j取值点上物理量值,p
i+1
为i+1时刻j取值点上物理量值,δt为时间历程曲线数据的输出时间间隔,ts为取值的起始时刻,te为取值的终止时刻。
[0130]
本实施例中,以1#取值点(element号:38473)为例,其等效应力的总时间历程曲线如图5中黑色曲线所示,从图中看出等效应力波形大致稳定后,其波形仍存在周期振荡,如图5中的红色曲线所示。在该取值点的27个振动的时间周期范围内对等效取算术平均值,取值点的起始时刻ts=1250
×
10-3
ms,终止时刻te=2710
×
10-3
ms,δt=5
×
10-3
ms,根据式(4),计算得到该点的稳态等效应力值为p1=10.4614mpa;采用相同的方法,得到其它17个取值点的稳态等效应力值分别为:p2=10.4130mpa、p3=10.7230mpa、p4=10.7656mpa、p5=11.1173mpa、p6=11.2284mpa、p7=10.9077mpa、p8=10.4523mpa、p9=10.7391mpa、p
10
=10.9369mpa、p
11
=10.5063mpa、p
12
=10.5752mpa、p
13
=10.4360mpa、p
14
=10.7552mpa、p
15
=11.0275mpa、p
16
=10.9144mpa、p
17
=11.3697mpa和p
18
=11.6663mpa。
[0131]
s63、将所有取值点的物理量稳态值进行算术平均,得到该构型爆源爆炸产生的地运动稳态效应值如式(5)所示,
[0132][0133]
式中,n为s43中确定的取值点个数。
[0134]
本实施例中,n=18,根据式(5)和s62中的计算结果,将所有取值点的物理量的稳态值进行算术平均,得到该构型爆源爆炸产生的地运动等效应力的稳态效应值为
[0135]
最后应说明的是:以上实施例仅用以说明本发明的技术方案,而非对其限制;尽管
参照前述实施例对本发明进行了详细的说明,本领域的普通技术人员应当理解:其依然可以对前述实施例所记载的技术方案进行修改,或者对其中部分或者全部技术特征进行等同替换;而这些修改或者替换,并不使相应技术方案的本质脱离本发明技术方案的范围。
当前第1页1 2 
网友询问留言 已有0条留言
  • 还没有人留言评论。精彩留言会获得点赞!
1