技术领域
[0001] 本
发明涉及一种射孔弹性能参数的欧拉获取方法，属于石油
测井技术领域。
背景技术
[0002] 随着石油勘探的不断深入，
地层及井况条件越来越复杂，对射孔弹的性能要求越来越高，因此获取射孔弹一些重要参数显得十分重要。
[0003] 石油射孔弹属于聚能装药的一种，是靠形成的高速射流来完成射孔。射孔弹的结构包括壳体、药型罩及炸药，作用过程一般经历炸药爆轰、药型罩压垮、射流形成及压垮等过程。由于射流形成过程在高温、高压及高速的情况下几十微秒内完成，所能获得的实验观测数据极为有限，并且实验周期长，成本高。随着计算机技术的发展，可以用计算机对射孔弹的作用过程进行数值模拟，动态地再现其形成过程，以获取射孔弹的一些重要参数。
[0004] 射孔弹作用过程的数学模型，包括基本的守恒方程组、应
力应变关系及
状态方程。
[0005] ①守恒方程组：
[0006] 忽略外力、外源及热传导，聚能射流形成问题的
流体弹塑性动力学模型，可用如下偏微分方程组表示：
[0008]
[0009] 动量守恒：
[0010]
[0012]
[0013] 其中，ρ为
密度，u为速度矢量，e为比内能，σ为
应力张量。
[0014] ②应力应变关系：
[0015] 在一定的加载速率条件下，材料的响应被分为体积
变形部分和应力偏量部分。一般地，应力分为体积应力和偏应力两部分：
[0016] σ＝-PI+S (4)
[0017] 式中，右边第一项-PI反映体积压缩和膨胀，从状态方程中求得；第二项S为应力偏量反映畸变，应力偏量变化率从应变率计算求得。
[0018] 在弹性范围内，用Cauchy应力率描述的广义Hooke定律为：
[0019]
[0020] 由于Cauchy应力率不满足客观性原理，不能用于本构方程中，因此，以Jaumann应力率代替Cauchy应力率，得到如下形式：
[0021]
[0022] 应变率张量：
[0023]
[0024] 旋转率张量：
[0025]
[0026] ③状态方程
[0027] 爆轰产物采用JWL状态方程：
[0028]
[0029] 式中，P为压力，E为爆轰产物的内能，V为爆轰产物的相对体积(即爆轰产物体积与初始体积之比)，A、B、R1、R2和ω为待定常数。
[0030] 空气采用理想气体状态方程：
[0031] p＝(ka-1)ρ·e (10)
[0032] 式中，ρ为密度，e为比内能，kα是空气的等熵指数，在计算中取ka＝1.4。
[0033] 对金属考虑在高温、高压、高应变率下表现的材料动态行为，采用状态方程描述，其一般形式如下：
[0034] P＝PH+ρe(Γ+Aμ) (11)
[0035] 其中：
[0036]
[0037]
[0038] 式中，Γ为 系数，A、S1、S2、S3为材料的特性常数，C为材料中的声速，ρ为当地密度，ρ0为初始密度，e为当地单位质量的比内能。
[0039] 射孔弹作用过程的数值模拟就是结合应力应变关系及状态方程来求解偏微分守恒方程组。然而由于三维数值模拟计算规模大，受制于单机内存及处理器主频限制，当前的单处理器性能无法满足计算需求，因此有必要在高性能集群上开展聚能射流形成大规模并行计算模拟研究。聚能射流形成属于大变形问题，对此问题，Euler方法有很大的优势。对于Euler方法的数值模拟，网格固定，物质在网格间流进流出，当前网格的输运会造成周边网格物理量的改变，隐含在输运
算法中的数据相关性及子区域间的关联性不易发现。如何确定网格内物质间的分界面以及如何确定当前网格向周
围网格的各种输运量，是Euler方法并行实现的关键和困难之处。
发明内容
[0040] 为了解决上述问题，本发明的目的在于提供一种射孔弹性能参数的获取方法，针对射孔弹作用过程进行三维数值仿真，开展大规模并行计算模拟研究，进而提高计算速度，增大计算规模，从而解决目前单处理器性能无法满足射孔弹作用过程中数值模拟的问题。
[0041] 本发明提供一种射孔弹性能参数的获取方法，该方法包括以下步骤：
[0043] 步骤2：收集确定步骤1所建立的模型中的相关参数；
[0044] 步骤3：根据步骤1建立的射孔弹模型确定计算域的大小，其大小为射孔弹模型尺寸的4倍以上，，并采用
正交的六面体网格来离散计算域，网格步长由药型罩壁厚确定，最大为药型罩壁厚的1/5；
[0045] 步骤4：将步骤3所得到的计算域沿三个坐标方向进行分区，得到的子区域数目与参与并行计算的计算机集群
节点数相对应(相等)；
[0046] 步骤5：获取步骤3所确定的各网格的物理量；
[0047] 步骤6：进行时间步内的计算，连续输出聚能装药性能参数，例如射流头部速度、射流长度及直径等，采用
可视化技术，再现射流的形成过程，以精确得到射孔弹侵彻深度及孔径大小。
[0048] 根据本发明提供的射孔弹性能参数的获取方法利用飞速发展的计算机技术成果，采用并行方法对射孔弹的作用过程进行数值仿真，极大提高了计算规模。采用此方法可以精确再现射孔弹的作用过程，获取射孔弹的一些重要参数。
附图说明
[0049] 图1示出
实施例1中数值模拟的计算几何模型示意图。
[0050] 图2(a)示出实施例1中所得闪光
X射线高速摄影图像，图2(b)示出实施例1中所得数值模拟结果图，其中图2(a)和图2(b)右侧示出t1＝25.5μs时刻的射流图像，左侧示出t2＝41.4μs时刻的射流图像。
具体实施方式
[0051] 下面结合附图，通过具体实施例对本发明进行详细说明。本发明的特点和优点将随着这些说明变得更为清楚、明确。
[0052] 在根据本方面的射孔弹性能参数的获取方法中，在步骤1中，根据工程中地层及井况条件，建立射孔弹计算几何模型，从而得到壳体及炸药的尺寸、以及药型罩的锥
角和厚度，供后续使用。
[0053] 在根据本方面的射孔弹性能参数的获取方法中，在步骤2中，收集确定步骤1所建立的模型中的相关参数，例如包括炸药、药型罩及空气的密度、初始比内能、状态方程参数，壳体及药型罩的屈服极限参数等，供后续使用。
[0054] 在根据本方面的射孔弹性能参数的获取方法中，在步骤3中，根据步骤1建立的射孔弹模型确定计算域的大小，其大小为射孔弹模型尺寸的4倍以上，优选4～5倍，并采用正交的六面体网格来离散计算域，网格步长由药型罩壁厚确定，最大为药型罩壁厚的1/5，优选为药型罩壁厚的1/10～1/5，进一步优选1/6～1/8。
[0055] 在根据本方面的射孔弹性能参数的获取方法中，在步骤4中，将步骤3所得到的计算域沿三个坐标方向进行分区，得到的子区域数目与参与并行计算的计算机集群节点数相对应(即子区域数目与计算机集群节点数相等)。
[0056] 在根据本方面的射孔弹性能参数的获取方法中，在步骤5中，获取步骤3所确定的各网格的物理量。
[0057] 具体而言，获取各网格的物理量通过采用算子分裂法求解偏微分守恒方程组来实现：
[0058] 将质量、动量、能量守恒方程统一写成如下形式：
[0059]
[0060] 式中，φ代表物理量密度(ρ)、能量(e)或速度(u)， 表示物理量φ在网格上的变化率； 为
对流项，即物理量φ在单位空间网格边界上的通量；等式右边H为源项。
[0061] 物理效应意义上的算子分裂法将方程(12)分裂成如下两个方程，同时在数值计算时也分成相应的两步完成：
[0062]
[0063]
[0064] 在优选的实施方式中，步骤5包括以下子步骤：
[0065] 子步骤5.1：并行求解方程(15)，即Lagrange步或压力效应步，不考虑对流项的影响，只考虑源项H(压力和偏应力的梯度效应)的作用，得到网格各物理量的中间值；
[0066] 子步骤5.2：在子步骤5.1更新中间物理量后，考虑方程(16)中对流项 的影响，即Euler步或输运步，通过计算质量、动量、能量在网格间的输运，对物理量在网格上进行重新分配；
[0067] 子步骤5.3：在子步骤5.2确定数据相关性及子区域关联性后，并行实现网格的输运。
[0068] 其中，在子步骤5.1中，具体求解方法为：对方程(15)离散，时间采用向前差分，空间采用一阶中心差分，此过程不产生数据相关性，各个网格并行独立求解，得到其中质量、动量、能量、密度、速度等物理量的更新值。
[0069] 在子步骤5.2中，由于该子步骤中网格的输运造成周边网格物理量的改变，会产生数据相关性及子区域间的关联性。在本发明中，确定数据相关性及子区域的关联性采用如下方法：
[0070] ①关于子区域间的关联性
[0071] 子区域的关联性来源于两个方面：网格间输运及数值格式带来的关联性。体现为子区域与邻近子区域的重叠程度，为了计算需要，往往要在子区域的边界增加一些网格来储存临近子区域的信息。
[0072] 网格物理量的改变量等于该网格在三个方向流进量与流出量的差，因此在Euler输运步里一个网格正确的更新包含该网格与邻近网格的输运。以单方向输运为例，假设流场速度沿着正方向，则单方向k网格一次正确的更新包含k-1和k网格的输运，即k网格物理量的改变等于k-1网格流进量(速度为负、方向值为负)与k网格流出量(速度为负、方向值为负)的差。因此对于一个三维的子区域，需要在子区域的左边、前边及下边各增加一层网格，且这些网格参与输运。
[0073] 由于在实际的更新运算中，采用的数值格式往往会需要周边的网格物理量。因此，考虑一个简单的数值格式：
[0074]
[0075] uk、uk+1、 分别表示第k、k+1个网格以及两个网格边界处的速度。
[0076] 则k网格更新需要k、k+1网格的物理量。因此对于一个三维的子区域，需要在子区域的右边、后边及上边各增加一层网格，这些网格不参与输运。
[0077] ②关于数据相关性分析
[0078] k-1网格的输运会造成k-1和k网格物理量的改变，因此当k网格输运时，k-1和k网格物理量已经改变。并且由于受制于CFL条件，当前网格的输运只会影响到周围26个网格，即当前网格的某种物质的输运量可写作如下函数关系式：
[0079] ΔVijk ＝ f(Δt，Δx，Δy，Δz，Vijk，jflagijk，uijk，Vi±1j±1k±1，iflagi±1j±1k±1，ui±1j±1k±1) (18)
[0080] 其中Δt为时间步长，Δx，Δy，Δz为网格步长，V为网格体积，iflag为网格介质标志。由于当前网格输运时，相邻及自身网格物理量已更新，若输运算法中V，iflag采用输运更新后的值，则算法中存在数据相关性。为消除数据相关性，采用V，iflag未更新的值，同时为了消除因此而带来的过量输运问题，将原串行算法中在一个三重空间循环下完成三个方向的输运，改为在一个三重空间循环下只进行一个方向的输运。
[0081] 在子步骤5.3中，在确定数据相关性及子区域关联性后，网格的输运可以并行实现。
[0082] 发明采用改进的“模糊界面法”实现重新分配过程。其基本思想是：对于含有多种介质、不确定各介质间界面的混合网格，把各介质占有网格的体积比作为模糊权重系数；根据模糊权重计算输运量；对介质进行模糊排序，决定输运优先级，按输运优先级进行输运。在上述建模和计算中应用模糊数学方法。
[0083] 具体实现过程包括：
[0084] ①网格分类
[0085] 将计算模型中含有的N种介质，编号为1，2，---，N，按照网格中含有介质种类个数不同将网格进行分类，得到 个类型的网格，分类如下：
[0086] 第零类，不含有介质的空网格：0(0)；
[0087] 第一类，含有一种介质的网格：1(1)、1(2)、---、1(N-1)及1(N)；
[0088] 第二类，含有两种介质的网格：2(1，2)、2(1，3)、---、2(1，N)，2(2，3)、2(2，4)、---、2(2，N)、---及2((N-1)，N)；
[0089] 第三类，含有三种介质的网格：3(1，2，3)、3(1，2，4)、---、3(1，2，N)、3(2，3，4)、3(2，3，5)、---、3(2，3，N)、---、3((N-2)，(N-1)，N)；
[0090] ------
[0091] 第N类，含有N种介质的网格：N(1，2，3，---，N)。
[0092] 当混合网格中含有三种以上的介质时，首先按照体积份额的大小对混合网格中各介质进行排序，将最大体积份额的介质序号规定为介质1，依次类推。然后删除第四种及其以后的介质，同部分的体积、质量、动量、能量等填充为介质1，这样最多余下三种介质。
[0093] ②制定输运方案
[0094] 首先根据模糊综合评判方法确定贡献网格中各介质的输运构型，以此确定各介质的输运优先级，优先级的评判是一个模糊综合评判，优先级高的介质优先输运，优先级相同时平均输运。
[0095] 当沿x轴正向输运时，设Vl1、Vr1、Vl2、Vr2、Vl3、Vr3分别为左边相临网格1介质、右边相临网格1介质、左边相临网格2介质、右边相临网格2介质、左边相临网格3介质、右边相临网格3介质的体积比(某介质的体积比是指网格中某介质的体积除以网格总体积所得到的百分比)。
[0096] 如果Vl1＞0则Il1＝1，否则，Il1＝0；如果Vr1＞0则Ir1＝1，否则，Ir1＝0；如果Vl2＞0则Il2＝1，否则，Il2＝0；如果Vr2＞0则Ir2＝1，否则，Ir2＝0；如果Vl3＞0则Il3＝1，否则，Il3＝0；如果Vr3＞0则Ir3＝1，否则，Ir3＝0；
[0097] 定义如下变量(sgn为符号函数)：
[0098] L1＝Il1·sgn(Vr1-Vl1)
[0099] R1＝Ir1·sgn(Vr1-Vl1)
[0100] L2＝Il2·sgn(Vr2-Vl2)
[0101] R2＝Ir2·sgn(Vr2-Vl2)
[0102] L3＝Il3·sgn(Vr3-Vl3)
[0103] R2＝Ir3·sgn(Vr3-Vl3)
[0104] 根据L1、R1、L2、R2、L3、R3的值，作出如表1所示的五种输运构型，这五种构型的输运优先级顺序为从高到低，即：①＞②＞③＞④＞⑤。例如，贡献网格中含有1和2两种介质时，若1介质属于构型③，2介质输运构型⑤，则1介质的输运优先级大于2介质，因此在输运时，优先输运1介质的输运量，不足部分以2介质进行补充。
[0105] 表1 介质输运构型表(k＝1，2，3)
[0106]
[0107]
[0108] 当贡献网格为纯网格时，输运因子β的值等于1；当贡献网格为混合网格时，网格内每种介质的输运因子β按下式计算：
[0109]
[0110] 式中，Vid、Vic分别为贡献网格和接受网格中同种介质的体积比。输运因子β表示从贡献网格向接受网格输运的该介质的体积占该介质总体积的百分比。
[0111] 具体实施时，模糊输运比的制定可以通过设定1、2及3三种介质具体所代表的物质来实现。
[0112] 在根据本方面的射孔弹性能参数的获取方法中，在步骤6中的优选实施方式中，子步骤5.1、子步骤5.2和子步骤5.3完成一个时间步内的计算，对于射孔弹问题一般要计算上万个时间步，子步骤5.1、子步骤5.2和子步骤5.3不断
迭代循环，设定输出聚能装药性能参数的步骤数和输出可视化图像的步骤数，从而连续输出聚能装药性能参数，例如射流头部参数如射流头部速度、射流长度及直径等，并连续输出可视化图像。
[0113] 以下通过范例性实例对进一步描述本发明。
[0114] 实施例1
[0115] 射孔弹锥角为60°。
[0116] (1)建立计算几何模型。
[0117] 数值模拟的计算几何模型如图1所示，药型罩锥角为60°，壁厚(壁厚均匀)为2.4mm，装药高度为21.9mm，装药直径为60mm。
[0118] (2)确定材料参数。
[0119] 炸药采用B炸药，金属罩采用45#
钢。B炸药、空气及45#钢分别采用公式(9)、(10)、(11)所示的状态方程，介质参数见下表2-表5。
[0120] 表2 B炸药爆炸性能参数
[0121]
[0122] 表3 B炸药炸药JWL状态方程参数
[0123]#
[0124] 表4 45 钢材料性能参数
[0125]
[0126] ρ0为初始密度，G为
剪切模量，C为材料声速。
[0127] 表5 45#钢材状态方程参数
[0128]
[0129] (3)划分有限差分网格。
[0130] 计算域大小为100mm×100mm×256mm，由于药型罩壁厚为2.4mm，三个方向的网格步长确定为0.4mm，采用等步长正交六面体网格离散计算域，共离散251×251×640个网格。
[0131] (4)选择节点及分配数据。
[0132] 网格数为4032万，选择四个节点共32个
进程参与计算，分区方式为2×2×8；主进程读取介质物理参数及网格剖分参数，并广播到其它进程。
[0133] (5)模糊界面应用于射孔弹聚能装药问题。
[0134] 根据模糊输运规则以及聚能射流形成的特点，设定45#钢为m1，B炸药为m2，空气为m3。
[0135] (6)启动计算。
[0136] 编程实现各网格的Lagrange步和Euler步的求解。输入上述三种介质的几何参数和材料参数，启动计算。共计算24000步，每100步输出一个可视化图片，及每15步输出射流头部信息。
[0137] (7)分析计算结果。
[0138] 图2(a)为闪光X射线高速摄影图像，图2(b)为数值模拟结果图，两图中右侧为t1＝25.5μs时刻的射流图像，左侧为t2＝41.4μs时刻的射流图像。
[0139] 下表6中示出数值模拟值及与试验值的比较结果。
[0140] 表6数值模拟结果及与试验值的比较
[0141]
[0142] 表6中的头部速度和射流直径取25.5μs和41.4μs时平均值。
[0143] 由图2和表6可以看出，数值模拟结果与实验结果总体上吻合的比较好。这也表明采用根据本发明的方法可以再现聚能射流形成过程，从而获取聚能装药的重要参数。
[0144] 以上通过具体实施方式和范例性实例已对本发明进行详细说明，不过这些实施方式和实例仅是说明性的，并不对本发明的保护范围构成任何限制，在不偏离本发明精神和范围的情况下，本领域技术人员能对本发明及其实施方式进行多种改进、等价替换或
修改，这些均落入本发明的保护范围内。本发明的保护范围以所附
权利要求书为准。
