0 引言
岩芯是研究地下地质和矿产情况重要的实物地质资料。通过对岩芯的观察、描述和研究，可以直接地了解地下岩层的岩性、岩相、物性和油气水的分布特征[1]。由于受天然岩芯来源的限制，在石油科学技术研究中，绝大多数试验都采用人造岩芯，人造岩芯的制造通常有石英烧结、磷酸铝石英烧结和环氧树脂压制胶结法3种[2]。笔者提出了一种新的岩芯制作方法——3D打印岩芯，即以模型的三维数据为基础，通过打印机喷嘴挤出材料，逐层打印增加材料，以实现人造岩芯的3D打印。在3D打印人造岩芯试验中，有很多因素制约着打印精度。在诸多因素中，打印装置的喷嘴部分设计至关重要。传统的喷嘴设计优化取得了一定的成果。张振等[3]建立了轴对称的喷嘴几何模型，利用Fluent软件对不同收缩角的喷嘴进行仿真，得出收缩角不同对喷嘴流场的轴线速度影响较大的结论；魏德军等[4]研究了喷嘴形状对喷射能力的影响，得出喷嘴入口形状为流线性和等变速的喷射效果为最佳的结论。笔者仅对岩芯3D打印装置中的一个重要组成部分——喷嘴进行设计和优化，对不同结构和几何参数的射流喷嘴进行仿真，研究各个参数对喷嘴射流速度的影响。所得结论可为3D岩芯打印装置的进一步改进提供参考。
1 3D打印装置简介岩芯3D打印装置分为3个部分：3D打印机、出料装置和供料装置。岩芯3D打印装置整体如图 1所示，其中供料装置连接到图中① 处。
|
| 图 1 岩芯3D打印装置整体图 Fig.1 Overall view of 3D core printing device |
1.1 3D打印机
岩芯3D打印的基础装置是一款Kossel_Mini型3D打印机(见图 2)。这款3D打印机结构简单，容易改进。工作过程如下：用CAD软件建立需要打印的模型，并生成STL格式文件，将生成的STL文件用切片软件(如Slic3r)转换为能够控制打印机各硬件的GCODE文件之后，通过上位机将GCODE文件传入控制器中，通过控制器驱动步进电机等部件，使打印机喷头按照指定的路线运动。
|
| 图 2 Kossel_Mini型3D打印机 Fig.2 Kossel_Mini 3D printer 1—电机；2—联轴器；3—轴承座；4—套管；5—喷嘴。 |
1.2 出料装置
出料装置由电机、轴承座、套管、螺杆和喷嘴等组成，采用单螺杆挤出原理，电机带动螺杆，实现材料的连续均匀挤出。把出料装置安装到3D打印机上，打印机带动出料装置按照指定的轨迹运动，同时喷嘴均匀出料，就可以完成材料的层层堆积，打印所要求的样本。图 3为出料装置整体图。
|
| 图 3 出料装置整体图 Fig.3 Discharge device |
1.3 供料装置
供料装置是为出料装置提供材料的推动装置，此处采用的是电动推杆，用电动推杆为出料装置提供均匀稳定的材料。
2 喷嘴的设计 2.1 喷嘴参数的确定喷嘴的作用是流体加速和改变所出料的孔径，从而改变打印精度。喷嘴设计的好坏在很大程度上影响打印精度，因为喷嘴为最终出料口，出料的粗细程度制约着打印样品的精度。喷嘴按结构的不同可以分为圆柱形喷嘴、锥形喷嘴、锥柱形喷嘴和锥柱扩散形喷嘴等。通过流体动力学分析软件对各种喷嘴进行模拟，发现在相同入口和出口条件下，锥形喷嘴出口速度分布比较均匀，出口紊动强度与壁面紊动强度分布比较小[5]，所以圆锥形喷嘴是目前应用较广泛的一种喷嘴。岩芯3D打印装置采用圆锥形喷嘴，决定其特性的参数有入口直径D、出口直径d、收缩角α和筒壁长L，如图 4所示。
|
| 图 4 圆锥形喷嘴参数图 Fig.4 Dimension of the conical nozzle |
实验室已加工的喷嘴具体参数为：入口直径D=0.010 m，出口直径d=0.002 m，收缩角α=45°，筒壁长L=0.010 m。以现有参数组作为对比，进行喷嘴的优化设计，从而改进喷嘴结构。
由于连接出料装置的套管直径已预先设计好，所以入口直径设为定值，于是研究另外3个参数对喷嘴性能的影响。又因为打印精度和打印机尺寸等限制，出口直径d取值范围为0.002~0.005 m，收缩角α取45°~60°，筒壁长L取0.01~0.02 m。
2.2 正交试验设计用响应面分析法对喷嘴进行方案设计和优化。响应面分析法采用多元二次回归方法作为函数估计的工具，将多因子试验中因子与指标的相互关系用多项式近似, 依此可对函数的响应面和等值线进行分析, 研究因子与响应面之间、因子与因子之间的相互关系[6]，具有试验周期短、求得的回归方程精度高及能研究几种因素间交互作用等优点。
试验根据中心组合设计原理，选取出口直径d、收缩角α和筒壁长L对出口速度v影响显著的3个因素，采取三因素三水平响应面分析法进行分析。当入口速度一定时，研究各个自变量(d、α、L)对出口速度的影响。表 1为喷嘴参数试验影响因素及范围。
| 影响因素 | 水平 | ||
| -1 | 0 | +1 | |
| L/m | 0.010 | 0.015 | 0.020 |
| α/(°) | 30 | 45 | 60 |
| d/m | 0.002 0 | 0.003 5 | 0.005 0 |
正交试验安排结果如表 2所示。
| 试验组 | L/m | α/(°) | d/m |
| 1 | 0.010 0 | 45 | 0.005 0 |
| 2 | 0.010 0 | 60 | 0.003 5 |
| 3 | 0.010 0 | 30 | 0.003 5 |
| 4 | 0.010 0 | 45 | 0.002 0 |
| 5 | 0.015 0 | 30 | 0.005 0 |
| 6 | 0.015 0 | 45 | 0.003 5 |
| 7 | 0.015 0 | 30 | 0.002 0 |
| 8 | 0.015 0 | 60 | 0.002 0 |
| 9 | 0.015 0 | 60 | 0.005 0 |
| 10 | 0.020 0 | 45 | 0.005 0 |
| 11 | 0.020 0 | 60 | 0.003 5 |
| 12 | 0.020 0 | 30 | 0.003 5 |
| 13 | 0.020 0 | 45 | 0.002 0 |
2.3 基于Fluent的喷嘴流体模型分析 2.3.1 Fluent分析
用Solidworks软件分别绘制不同参数的喷嘴模型，导入到ANSYS软件中划分网格，并定义入口、壁面及出口，导入到Fluent软件中进行流体分析。
将网格导入Fluent软件中，设置为3D，用多相流模型来模拟水泥流体[7]，定义液体水和水泥材料，水泥材料密度为3 100 kg/m3，黏度为10-3 kg/ (m·s)。边界条件给定入口为速度入口，出口为outflow，壁面无滑移，设置迭代次数为100，设置好相关参数后进行分析，最终得出不同正交试验组喷嘴出口流体速度，填入Design-Expert软件中进行结果分析，模拟后参数如表 3所示。由表可以看出现有参数组的出口速度情况(表 3中第4组数据为实验室已有喷嘴尺寸)，得出现有参数组出口速度v=0.172 0 m/s，以便与后面的优化结果进行对比。
| 试验组 | L/ m | α/ (°) | d/ m | v/ (m·s-1) |
| 1 | 0.010 0 | 45 | 0.005 0 | 0.031 6 |
| 2 | 0.010 0 | 60 | 0.003 5 | 0.044 6 |
| 3 | 0.010 0 | 30 | 0.003 5 | 0.054 0 |
| 4 | 0.010 0 | 45 | 0.002 0 | 0.172 0 |
| 5 | 0.015 0 | 30 | 0.005 0 | 0.030 8 |
| 6 | 0.015 0 | 45 | 0.003 5 | 0.057 0 |
| 7 | 0.015 0 | 30 | 0.002 0 | 0.160 0 |
| 8 | 0.015 0 | 60 | 0.002 0 | 0.139 5 |
| 9 | 0.015 0 | 60 | 0.005 0 | 0.024 4 |
| 10 | 0.020 0 | 45 | 0.005 0 | 0.031 3 |
| 11 | 0.020 0 | 60 | 0.003 5 | 0.044 4 |
| 12 | 0.020 0 | 30 | 0.003 5 | 0.054 3 |
| 13 | 0.020 0 | 45 | 0.002 0 | 0.171 2 |
2.3.2 网格无关性
J.H.FERZIGER等[8]指出，随着网格质量的提高，性能预测的误差将会逐渐降低，为确保数值计算的正确性，需对网格无关性进行检查。
现以其中一组数据为代表，即表 3中第6组数据，分别选择4种数量的网格导入到Fluent软件中进行计算，其他参数设置保持一致，单纯看网格数量对结果的影响。网格情况及模拟结果见表 4。
| 组别 | 1 | 2 | 3 | 4 |
| 网格数 | 62 926 | 555 795 | 893 805 | 2 415 575 |
| v/(m·s-1) | 0.055 6 | 0.057 3 | 0.057 1 | 0.057 0 |
由表 4可以看出，随着网格数量的增加，出口速度发生了一定变化，当网格数量达到55万以后，随着网格数的增加，出口速度整体呈现下降趋势；当网格数量选定为较大的数量，模拟的出口速度大小并无明显的变化，且下降幅度越来越小，在网格数约为240万时能够得到出口速度与网格无关性的结果。
2.4 试验模型的建立采用响应面法对试验数据进行统计分析，用多元二次回归方程来拟合因子和响应值之间的函数关系，得出结果。软件依据拟合误差最小原则选取理想模型对表 3中的试验数据进行拟合。图 5为拟合模型的学生化残差(学生化残差是评价各观测数据相对于回归拟合曲线是否为异常点的依据，残差值越大，图中的残差点偏离直线越远，表示拟合效果越差)分布图。从图可以看出，其残差各点的分布几乎在一条直线上，模型拟合较好。
|
| 图 5 模型的学生化残差分布图 Fig.5 Internally studentized residuals of the model |
依据拟合误差最小的原则，选取理想模型对表 3中的试验数据进行拟合，最终得到用二次方程拟合效果最好，以v代表喷嘴出口流体速度，建立的拟合方程为：
|
(1) |
将出口速度的模拟值和方程拟合值进行比较，结果见图 6。2条曲线虽然存在差别，但大体比较接近，说明方程的拟合效果较好，拟合方程合理。
|
| 图 6 模拟值和拟合值对比图 Fig.6 Comparison of simulation and fitting values |
拟合方程中出口速度的方差分析结果如表 5所示。由表 5可看出各因素对出口速度的影响程度：出口直径d＞收缩角α＞筒壁长L；也可以根据数学模型中各影响因素的系数的绝对值大小来研究各影响因素对喷嘴出料速度的影响，可以得到相同的结果。可见，出料口直径对出料口速度影响最大。
| 变差来源 | 平方和 | 自由度 | 均方 | 虚拟项F0 |
| 回归 | 0.042 | 9 | 0.004 7 | 195.030 |
| 因素L | 1.250×10-7 | 1 | 1.250×10-7 | 5.165×10-3 |
| 因素α | 2.668×10-4 | 1 | 2.668×104 | 11.030 |
| 因素d | 0.034 | 1 | 0.034 | 1 421.510 |
| 因素Lα | 6.250×10-8 | 1 | 6.250×10-8 | 2.583×10-3 |
| 因素Ld | 6.250×10-8 | 1 | 6.250×10-8 | 2.583×10-3 |
| 因素αd | 4.970×10-5 | 1 | 4.970×10-5 | 2.050 |
| 因素L2 | 2.819×10-5 | 1 | 2.819×10-5 | 0.316 |
| 因素α2 | 4.434×104 | 1 | 4.434×10-4 | 18.320 |
| 因素d2 | 7.405×10-3 | 1 | 7.405×10-3 | 306.000 |
2.5 参数优化
岩芯3D打印装置的喷嘴出口速度优化问题属于多变量、有约束及非线性的最小化问题，其目标函数为：
|
(2) |
式(2) 中，筒壁长L、收缩角α和出口直径d的取值范围如下：0.01 m≤L≤0.02 m, 30°≤α≤60°, 0.002 m≤d≤0.005 m。
非线性不等式约束为：
|
(3) |
针对试验的特点，在优化过程中选取序列二次规划法，即SQP算法。该算法不仅具有全局收敛性, 而且具有超线性收敛的速度。SQP法实现主要分3个阶段：① 拉格朗日函数的Hessian矩阵(Hk)的更新；② 二次规划QR问题的求解[9]; ③ 通过求解一系列的二次规划问题以逼近原方程，形成QP子问题，即：
|
(4) |
|
(5) |
|
(6) |
其中：xk是当前问题的迭代点，得到一个搜索方向dk，然后经过线搜索求得步长ak，于是下一个迭代点xk+1=xk+akdk。该优化问题给定初始点为x0=[0.01，30，0.002]，得到搜索方向d0，从而进行下一次迭代，直到满足迭代终止条件。用SQP方法求解问题的计算流程如图 7所示。
|
| 图 7 SQP方法求解问题的计算流程 Fig.7 The solution process of SQP method |
应用MATLAB优化工具箱的SQP法求解该优化问题时，可直接调用fmincon函数(过程略)。由于试验装置的3D打印机打印速度的限制，喷嘴的出口速度要配合打印机的打印速度，同时考虑打印精度等因素的影响。装置要求喷嘴出口速度在一定范围内趋于最小，即0.01 m/s≤v≤0.04 m/s，得到的优化结果为：筒壁长L=0.015 m，收缩角α=56.045°，出口直径d=0.005 m，出口速度v=0.01 m/s。
将优化后的结果与现有的参数组(表 3第4组)进行对比，发现优化后出口速度v=0.01 m/s，现有参数组出口速度为0.172 0 m/s，优化后的结果明显接近要求的出口速度，优于现有的参数组。最后得出最优解为：L=0.015 m, 收缩角α=56°，出口直径d=0.005 m，出口速度v=0.01 m/s。
3 结论(1) 岩芯是研究地下地质和矿产情况重要的实物地质资料。人造岩芯是试验研究必不可少的物质。3D打印技术正在发展起来，不同的打印材料正在开发中。
(2) 用响应面分析法对喷嘴进行设计及优化，可以极大程度地减少工作量，提高效率。
(3) 在喷嘴的几个参数中，对喷嘴出口速度影响最大的是出口直径，是今后设计喷嘴时主要考虑的参数。
(4) 通过SQP算法的优化求解，并与现有参数组进行对比，最终得到的喷嘴最优参数组合为L=0.015 m，收缩角α=56°，出口直径d=0.005 m，出口速度v=0.01 m/s。
| [1] | 梁万林. 人造岩心制备技术研究[J]. 石油仪器, 2008, 22(2): 72–76. |
| [2] | 皮彦夫. 特殊人造岩心的类型及应用[J]. 科学技术与工程, 2010, 10(31): 7724–7727. DOI: 10.3969/j.issn.1671-1815.2010.31.024 |
| [3] | 张振, 章巧芳. 高压喷嘴的射流仿真研究[J]. 机电工程, 2013, 30(2): 185–187. |
| [4] | 魏德军, 曹仲文, 李萍. 喷嘴形状对喷射能力的影响[J]. 辽宁化工, 2006, 35(7): 416–427. |
| [5] | 施春燕, 袁家虎, 伍凡, 等. 射流抛光喷嘴的设计[J]. 光电工程, 2008, 35(12): 131–135. DOI: 10.3969/j.issn.1003-501X.2008.12.027 |
| [6] | 陶有俊, DANIELT, 赵跃民, 等. 采用Design-Expert设计进行优化Falcon分选试验[J]. 中国矿业大学学报, 2005, 34(3): 343–348. |
| [7] | 陈卓, 曾良才, 湛从昌, 等. 基于Fluent的钻井纠斜液压系统导向套叶片形状优化选择[J]. 武汉科技大学学报(自然科学版), 2015, 38(4): 279–281. |
| [8] | FERZIGER J H, MILOVAN P. Computational methods for fluid dynamic[M]. Berlin: Springer Berlin Heidelberg, 2002. |
| [9] | 叶秉良, 郭绍义, 戚金明. 基于SQP法的斜齿圆柱齿轮减速器可靠性优化设计[J]. 浙江理工大学学报, 2008, 25(2): 187–190. |