模型的建立基于以下假设：①不考虑喷嘴壁面与超临界CO₂之间的换热；②磨料颗粒为规则的刚性球体，无质量损失且与射流流体无质量交换；③颗粒与壁面的碰撞视为弹性碰撞；④本研究主要探究颗粒轨迹和旋转射流流场特征的变化，不考虑因射流冲击造成的几何形变。

超临界CO₂与水不同，其密度、黏度等物理性质随温度和压力的变化而变化，进而影响射流流场分布规律及颗粒运动特征。数值模拟计算需采用针对可压缩气体与超临界态CO₂流体的状态方程，考虑温度、压力对流体密度、黏度的影响。本研究采用美国国家标准与技术研究所 (NIST) 开发的隐式方程来计算超临界态CO₂在喷射过程中的物理性质变化。与Peng-Robinson方程以及维里系列状态方程相比，Span-Wagner模型通过引入无因次亥姆霍兹自由能来计算超临界CO₂状态参数，方程适用温度与压力范围更广、计算精度更高，亥姆霍兹自由能无因次表达式为式(1)^[22]。

基于无因次表达式，可推导出真实气体的压力与定压热容变化方程如式(2)~式(3)所示^[22]。

超临界CO₂流体黏度与导热系数分别采用Fenghour方法与Vesovic方法计算，黏度和导热系数表达式分别如式(3)~式(4)所示^[23-24]。

式中：ϕ为无量纲亥姆霍兹能，ϕ₀是理想气体条件下亥姆霍兹能，ϕ_r是残余亥姆霍兹能，δ是流体密度与临界条件下密度比值，τ是流体温度与临界温度比值；p是流体压力，MPa；ρ是流体密度，kg·m⁻³；R是理想气体常数, J·(mol·K)⁻¹；T是流体温度，K；M是相对分子量，kg·mol⁻¹；C_p是定压热容，J·(kg·K)⁻¹；μ是流体黏度，Pa·s；λ是导热系数，W·(m·K)⁻¹；μ₀和λ₀为零密度条件下流体黏度与导热系数；Δμ_c和Δλ_c流体接近临界点时黏度和导热系数的增量。

本研究对象为高雷诺数下的冲击射流流场，并且由于旋流特性，在外流场中存在剧烈的涡旋流动，因此直射流流场条件下常用的standard k-ε两方程湍流模型适用条件受限。本研究采用Realizable k-ε湍流模型计算旋流射流的流场分布特征，Realizable k-ε湍流模型已被证实在多种复杂流动条件下具有较高的准确性，包括喷射流动、旋转剪切流、自由流和分离流等，模型湍动能k与湍流耗散ε率的计算公式如式(6)~式(8)所示^[25]。

式中：k是湍流动能，J；u_j是速度分量，m·s⁻¹；μ_t是湍流粘度系数，Pa·s；σ_k是普朗特数在k上的分量；G_k是湍流动能k因平均速度梯度的引出项，kg·J·(m³·s)⁻¹；G_b是湍流动能k因浮力的引出项，kg·J·(m³·s)⁻¹；Y_M是可压缩湍流中波动膨胀对总耗散率的贡献，kg·J·(m³·s)⁻¹；S_k和S_ε是自定义的源项，kg·J·(m³·s)⁻¹；ε是湍流耗散率；σ_ε是普朗特数在ε上的分量；C₁，C₂，C_1ε，C_3ε为常数。

关于固-液两项流动的模拟，常见的方法有2种：一种是将固体颗粒视为连续介质或拟流体，另一种是将固体颗粒视为离散相介质。在本研究模型中，固体相颗粒体积分数较小，为了捕捉旋转射流颗粒的运动轨迹和速度，采用离散相模型 (discrete phase model，DPM) 进行模拟。通过计算每个颗粒的受力平衡方程，推导颗粒的速度和轨迹分布，单个磨粒的受力平衡方程如式(9)所示^[26]。

式中：F_D (u_f −u_p)为颗粒的拖拽力项；$ \dfrac{{g\left( {{\rho _p} - {\rho _f}} \right)}}{{{\rho _p}}} $为重力项，F_x为颗粒附加力项，F_D表达式如式(10)所示。

为了定量对比不同参数条件下旋转磨料射流冲击下岩体的损伤程度，引入了颗粒侵蚀模型，该模型常用于评估管道运输过程中的颗粒对管道壁面的磨损情况^[27]。在本研究中的射流冲击模型中，虽然侵蚀模型不能直接反映实际岩体的破碎效率，但它仍然可以为我们评估颗粒撞击效果提供参考，磨料侵蚀方程可以表示为式(11)^[28]。

式中：u_p和u_f分别为颗粒和流体速度，m·s⁻¹；ρ_p和ρ_f分别为颗粒和流体密度，kg·m⁻³；F_D (u_f −u_p)为拖曳力项，m·s⁻²；$ \dfrac{{g\left( {{\rho _p} - {\rho _{f} }} \right)}}{{{\rho _{p} }}} $为重力项，m·s⁻²；F_x为附加力项，模型中考虑了虚质量力和压力梯度力，m·s⁻²；d_p为颗粒直径，m；C_D为两相间拖曳力系数；Re是雷诺数；N_p为颗粒数量；m_p为质量流量，kg·s⁻¹；C(d_p)是粒度分布函数；f(α)是颗粒冲击角度函数；v_p是颗粒速度，m·s⁻¹；b(v)是颗粒速度函数；A_face是冲击面积，m²。

基于实际喷嘴几何结构，提取流道特征建立了前混式旋转磨料射流和锥形磨料射流流场几何模型。如图1所示，射流冲击流场可分为喷嘴内流场和外流场2部分组成。其中旋转射流喷嘴内流场较为复杂，主要包含入口段、叶轮段和收缩发展段，其中流体首先经叶轮段导流后产生周向速度，再经历收缩发展段加速后呈螺旋放射状进入外流场；锥形喷嘴由入口段、收缩段和直管段3部分组成，高压流体从喷嘴入口流入，经直管段加速后形成高速直射流进入外流场，喷嘴内收缩夹角为30°，喷嘴出口直径5 mm。外流场区域设定为长200 mm、宽200 mm、高50 mm的长方体，上部与喷嘴相连，下部设定为冲击壁面，四周为流场压力出口，高速射流冲击至底部平面后产生漫流，并沿径向方向流出外流场。

本研究主要探究超临界CO₂旋转磨料射流在外流场中的颗粒携带特征，为提高外流场网格质量，采用结构化网格分块构建拓扑结构，将整体网格划分为喷嘴内流场网格和外流场网格2部分，两者之间采用“Interface”界面连接。旋转喷嘴内流场的叶轮段由多组弯曲流道构成，结构较为复杂，将其进一步分段切割，采用旋转对称的方式对模型进行网格离散。在外流场中，旋转射流的径向扩散特征会形成复杂的涡旋结构，为了捕捉旋转射流流场结构变化，对外流场射流剪切层内的网格进行了局部加密如图1(a)所示，在此分块网格划分方案下，外流场网格质量可达到1.0，从而精细表征射流流场的结构特征。依据此网格划分方案，外流场网格尺寸是影响模型精度的关键因素，因此针对外流场进行了网格无关性验证，以流体速度作为标定参数，得到了如图2所示的速度随网格尺寸变化规律，综合考虑计算效率与计算精度，外流场加密区网格尺寸设定为1 mm。此外，考虑边界层影响，对喷嘴壁面网格和冲击面底部网格也进行了加密处理，加密完成后旋转射流模型整体流场网格数量为88×10⁴，锥形射流模型网格数量为59×10⁴。

在本研究中射流模型采用压力入口和压力出口边界条件，压力入口即泵注压力设定为50 MPa，压力出口即环境围压，对于页岩气和深部煤层气等非常规油气藏^[29]，井底温度和压力都显著高于CO₂的临界条件，井底围压设定为30 MPa，环境温度为340 K，超临界CO₂物性参数变化采用NIST真实气体模型。喷嘴内部壁面为无滑移绝热固壁条件，颗粒撞击后发生弹性碰撞，不产生能量损失。在进行流场模拟时，为保证模型的稳定性和准确性，首先计算不添加磨料颗粒的超临界CO₂射流流场特征，待流场稳定后将此状态作为初始条件导入超临界CO₂磨料射流模型，而后采用离散相 (DPM) 模型定义磨料颗粒属性，假设磨料颗粒为圆球形，本研究采用油田上常用的支撑剂类型和粒径作为磨料颗粒，包括石英砂、陶粒、自悬浮支撑剂和超低密度支撑剂，研究不同颗粒的跟随特性，详细参数见表1。在离散相中设置磨料颗粒的材料、粒径以及进入流场的位置、方向和速度等参数，其中磨料进入流场时的速度大小与完全发展的高速射流相比可忽略不计，掺混后的颗粒速度完全由高速射流决定，因此在本模型中，忽略了磨料颗粒入口速度的影响。

在基准算例条件下 (石英砂，40/70目，0.10 kg·s⁻¹) ，对超临界CO₂旋转磨料射流流场进行了模拟，得到了射流速度场、颗粒轨迹以及射流冲击压力和颗粒冲蚀分布云图 (图3) 。如图3(a)所示，由于旋转喷嘴内部特殊的叶轮结构，流体速度在外流场呈现旋转扩散状，这种射流结构可有效扩大冲击面积 (图3(c)) ，产生平行于岩石表面的切向应力，从而提高破岩效率。受旋转射流流场影响，颗粒轨迹同样呈扩散状进入外流场，这进一步增大了射流有效冲蚀面积 (图3(b)) 。在磨料射流中，颗粒侵蚀是岩石破碎的主要机制，本研究通过引入冲蚀模型 (公式(10)) 对磨料颗粒的撞击效果进行评估。如图3(d)所示，超临界CO₂旋转磨料射流的冲蚀云图呈现环状特征，在射流中心可以观察到一个弱侵蚀带，环状图案主要是由于粒子轨迹的扩散，导致较少的粒子撞击到中心区域，这一现象与旋转水射流的破岩实验结果一致^[31]。在相同的边界条件下，与旋转磨料水射流相比，超临界CO₂得益于其较低的黏度和密度，射流过程中的压力损失较小。因此，超临界CO₂射流速度、冲击压力和最大冲蚀速率均高于水射流。但是超临界CO₂的特殊物性也导致了磨料颗粒冲蚀范围的降低，由公式(8)和公式(9)可知，流体对颗粒的拖曳力受流体密度、流体黏度和流体速度综合影响。超临界CO₂较低的黏度和密度使得其对颗粒的携带能力降低，虽然颗粒的扩散范围受限，但同时这也增大了颗粒在同一位置的冲击频率，提高破岩效率。作为对照，锥形射流流场形态及冲蚀特征如图4所示，与旋转射流不同，锥形射流的速度和颗粒轨迹更加集中，虽然冲蚀面积显著缩小 (减小17倍) ，但喷射过程中的能量耗散更低，在相同压差条件下，锥形射流冲击压力较旋转射流提高29.34%。此外，与旋转射流变化规律相同，锥形水射流的在颗粒速度、冲击压力以及冲蚀速率等方面均低于超临界CO₂射流，再一次凸显了超临界CO₂流体的独特优势。

颗粒粒径是影响射流冲蚀效果的关键因素，本节对比了不同磨料粒径下旋转射流和锥形射流的速度分布和冲蚀特征。如图5所示，绘制了不同粒径下最大射流速度和最大颗粒速度的分布曲线。由图可知，随磨料粒径增加，旋转磨料射流的颗粒速度略有提升，这是由于磨料粒径增大，单个颗粒的受力面积增加，流体对颗粒的拖曳力增大。但在本研究模型中由于设定颗粒质量流量为定值，粒径变化对流体整体做功影响较小，因此颗粒速度变化较小，特别是锥形喷嘴的颗粒速度几乎保持不变，这与张世昆等^[32]关于液氮磨料射流的研究结果一致，而在旋转喷嘴中，由于其较长的流道结构，为颗粒加速提供了更长的距离，从而放大了粒径变化对速度的影响。

虽然粒径变化对喷嘴内的颗粒加速过程影响较小，但当颗粒进入外流场后受超临界CO₂携岩能力限制，其对底部的冲蚀特征和强度变化作用显著。如图6(a)所示，在超临界CO₂旋转磨料射流冲击作用下，随磨料粒径增加，射流冲蚀强度逐渐增加，但冲蚀面积逐渐缩小，这是由于随着磨料粒径增加，单颗粒质量增加，导致超临界CO₂对颗粒的携带能力逐渐下降，因此颗粒轨迹的扩散性减弱，而在同一位置处的颗粒冲蚀频率逐渐增加，冲蚀强度提升。旋转磨料水射流的变化趋势与超临界CO₂基本相似，但是由于水对磨料颗粒的携带能力相对较强，随粒径增加依然能保持较大的冲蚀面积 (图6(b)) 。在本模型设置条件下，优选30/50目作为旋转水射流的磨料粒径，可在提高冲蚀速率的基础上保持较大的冲击面积。锥形磨料射流则与旋转磨料射流表现出相反的规律，如图7所示，随磨料粒径增加，锥形磨料射流在冲蚀强度和冲蚀范围均呈现减小的趋势。这是由于锥形射流和旋转射流不同的流场结构特征，锥形磨料射流的颗粒轨迹更加集中，在相同质量流量条件下，大目数 (小粒径) 磨料的颗粒数量更多，即在喷射过程中颗粒对岩体的冲击频率增加，而磨料速度几乎不变，依据本研究引入的冲蚀模型，冲蚀强度受冲击频率和冲击速度影响显著，因此在大目数条件下锥形磨料射流表现出较高的冲蚀速率。与之相对应的旋转射流由于颗粒轨迹的扩散范围较广，虽然增加了冲击频率，但在较大的冲击范围上依然表现出较低的冲蚀速率。因此，在本模型设置条件下，超临界CO₂锥形磨料射流优选70/140目、锥形水射流40/70目。

为优选超临界CO₂磨料射流颗粒类型，以密度为指标，对比了现场常用的4类磨料颗粒，探究了不同磨料类型对超临界CO₂磨料射流速度分布和冲蚀速率特征的影响。如图8(a)所示，颗粒密度变化对超临界CO₂射流影响显著，在相同射流压差和颗粒质量流量条件下，随颗粒密度增加，流体速度和颗粒速度均逐渐降低。这是由于颗粒密度增大，质量增加，但是流体物性参数和颗粒直径保持不变，即流体对颗粒的拖曳力不变，因此加速颗粒至相同速度所需的时间增加。此外，为定量表征流体对颗粒的携带能力，定义速度比值 (r_v) =最大颗粒速度/最大流体速度，4种情况下的速度比值分布为99.9%、93.3%、93.0%和92.8%。超低密度支撑剂的速度变化与其他3种磨料颗粒显著不同，其颗粒速度与流体速度更为接近。这是由于超临界CO₂的低黏特性使其颗粒携带能力较水等流体介质更弱，因此在对高密度颗粒加速时会存在一定的滞后，加之旋转喷嘴较长的叶轮流道结构进一步增强了这一特性，但是由于超低密度支撑剂的密度与流体更为接近，因此其在喷嘴内的运移过程中能够保持较好的跟随特性。与旋转喷嘴不同，锥形喷嘴流道结构较为简单，由图8(b)速度分布曲线可知，随颗粒密度增加，流体速度和颗粒速度未发生明显变化，即在锥形喷嘴条件下，流体对颗粒的携带能力受颗粒密度影响较小。

图9为不同磨料类型下的颗粒冲蚀云图，由公式(10)可知，颗粒侵蚀强度主要与颗粒数量、质量流量以及颗粒冲击速度等参数有关。虽然锥形射流的颗粒出口速度小于旋转射流，但是由于颗粒轨迹更加集中，在同一位置的颗粒冲击频率更高，因此表现出较高的冲蚀强度。在旋转射流中，随颗粒密度的增加，颗粒的最大冲蚀速率不断增大，冲蚀面积逐渐降低。这是由于颗粒密度较低时，流体对颗粒的携带能力更强，随旋转射流的扩散范围更大，但同时由于单颗粒质量较小，虽然颗粒出口速度较大，然而单颗粒的冲击强度却不高，因此随颗粒密度增加，冲蚀速率整体表现出逐渐增大的趋势。因此，实际现场应用时，磨料类型的选择需综合考虑冲蚀速率与冲蚀范围之间的平衡。锥形磨料射流的冲蚀面积同样随颗粒密度的增加而降低，但是最大冲蚀速率与旋转射流表现出不同的变化趋势，密度为2 485 kg·m⁻³时冲蚀速率最大，而密度为3 300 kg·m⁻³时冲蚀速率最低，这主要是由于锥形射流冲蚀面积更集中，在低磨料密度条件下，虽然颗粒质量较小，但是由于更大的冲击频率和冲击速度，整体显现出更高的冲击强度，而高密度颗粒则相反。因此，当锥形喷嘴用于超临界CO₂喷射压裂时，选择低密度或超低密度磨料较为合理。

如前文所述，超临界CO₂兼有气体和液体的双重特点，其密度接近液体，而粘度近似于气体，扩散系数是液体的近百倍，这些特性为超临界CO₂的应用提供了显著优势。然而，这也意味着超临界CO₂在颗粒携带能力上相对较弱。在现场施工过程中，磨料浓度即颗粒质量流量也是影响射流冲蚀强度的关键因素，本节对比了不同颗粒质量流量对超临界CO₂射流流场以及冲蚀强度的影响。如图10和图11所示，随质量流量增加，射流速度和颗粒速度总体上呈线性减小。这是由于颗粒质量流量增加，单位时间内通过喷嘴的颗粒数量增多，喷射过程中单位体积流体的负荷量提高，导致流体动量损失增大，进而导致射流速度减小，颗粒速度相应降低。此外，由速度变化幅度可以看出，由于旋转射流流道结构更为复杂，其随颗粒质量流量变化较锥形射流更为敏感，在现场应用时更宜选择低磨料浓度场景。

虽然质量流量较低时可以获得更高的喷嘴出口速度，增大单个颗粒对井底岩石的冲击作用，但是较低的磨料浓度，同样会降低颗粒对岩石的冲击频率，进而影响磨料射流冲击破岩的效果。依据本研究引入的冲蚀模型 (公式(10)) ，计算了不同颗粒质量流量下的冲蚀云图 (图12) 。由图可知，在本研究模型设定条件下，旋转射流和锥形射流随颗粒质量流量增加，射流冲蚀范围和冲蚀强度均呈现增大趋势，在较高的冲击频率下弥补了颗粒速度损失带来的影响。但是不能为了提高破岩效率而过量增加磨料颗粒浓度，由图11(a)可知，当颗粒质量流量超过0.20 kg·s⁻¹时，旋转喷嘴内易出现颗粒堵塞，井下发生砂堵的风险增大。同时，旋转射流的流场形态出现收缩，颗粒轨迹的扩展特性减弱，颗粒侵蚀面积减小。因此，在现场应用中应当合理选择射流颗粒质量流量，在获得较好破岩效果的同时降低砂堵和喷嘴磨蚀的风险。

通过上文对比超临界CO₂旋转磨料射流和锥形磨料射流的冲蚀性能，可以发现旋转磨料射流在扩大冲蚀范围方面具有优势，而锥形磨料射流表现为深穿透。因此，超临界CO₂旋转磨料射流可主要应用于径向水平井套管开窗、水平井造穴 (造腔) 应力释放、碳酸盐岩储层近井地带酸化解堵等，采用清洁流体实现控近扩远。超临界CO₂锥形磨料射流可主要应用于超临界CO₂喷砂射孔-压裂联作，采用超临界CO₂作为射孔液和压裂液 (或前置液) ，免去了炮弹射孔和化学添加剂对于储层的伤害和污染，为环境友好的新型流体介质。

针对于磨料粒径的优选，超临界CO₂旋转磨料射流需要兼顾冲蚀速率与冲蚀范围之间的平衡关系，根据不同的工程应用场景和工程目标，以及喷枪上喷嘴的排布和水力参数进行优化；旋转水射流的磨料粒径优选30/50目，可在提高冲蚀速率的基础上保持较大的冲击面积。超临界CO₂锥形磨料射流对小粒径的磨料颗粒携带性能较好，建议选择70/140目石英砂，可兼顾较强的冲蚀速率与较大的冲蚀范围；锥形磨料水射流应用于套管射孔时，建议选择40/70目，具有较集中的冲击能力和较强的冲蚀强度。对于磨料类型的选择，超临界CO₂旋转磨料射流同样需要兼顾冲蚀速率与冲蚀范围之间的平衡关系，针对不同地质储层，选择与其相适应的磨料类型；而超临界CO₂锥形磨料射流建议采用低密度或超低密度磨料进行套管射孔。针对于质量流量，无论是超临界CO₂旋转磨料射流还是超临界CO₂锥形磨料射流，质量流量越大越有利于提高冲蚀速率和冲蚀范围。因此，在井底砂浓度确定的条件下，在限压范围内尽可能提高排量，以提高冲蚀效果。

1) 超临界CO₂旋转磨料射流在射流速度、冲击压力和磨料冲蚀速率等方面均优于磨料水射流。与常规锥形射流相比，旋转磨料射流能显著扩大冲蚀面积，适用于径向水平井套管开窗、水平井造穴应力释放、碳酸盐岩储层近井地带酸化解堵等工程场景。超临界CO₂锥形磨料射流能量更为集中，其深穿透的性能更适用于喷砂射孔-压裂联作，有效减少储层伤害和污染。

2) 针对于磨料参数的优选，超临界CO₂旋转磨料射流需要兼顾冲蚀速率与冲蚀范围之间的平衡关系，根据不同的工程应用场景和工程目标，选择合适的磨料类型和粒径。喷射过程中通过增加磨料粒径和密度，可以提升超临界CO₂旋转磨料射流的冲蚀强度，但同时也会受到超临界CO₂携岩能力的限制，导致射流冲蚀范围逐渐降低。提高磨料浓度，单位体积流体负荷增大，虽然流体速度和颗粒速度会降低，但较高的颗粒冲击频率能够有效弥补速度损失，从而增强射流的冲蚀范围和强度。然而质量流量过高时，喷嘴内易出现颗粒堵塞，井下发生砂堵的风险增大。因此，在现场应用中应当合理选择射流颗粒质量流量，在获得较好破岩效果的同时降低砂堵和喷嘴磨蚀的风险，同时，在井口限压范围内尽可能提高排量，以提高冲蚀效果。