柯年峰，孙伟，刘小龙，张华，李辉

（730070 甘肃省 兰州市 甘肃农业大学 机电工程学院）

0 引言

随着农业种植技术的发展，地膜覆土这一项技术得到了广泛应用，而膜上覆土技术会造成覆盖的泥土发生板结，引起缺苗或者烧苗的情况。使用打孔地膜会有效改善膜上土壤板结，针对地膜打孔问题设计了一种热熔式成孔器。热熔式成孔器是基于工具与工件之间传热现象，热量从高温热熔式成孔器传递到地膜上，达到对地膜形状和尺寸的加工。有热量传递就会引起温度的变化。地膜结构由于温度变化导致其产生热变形并形成热应力。同时热变形也会改变其边界，又导致地膜发生破裂[1]。处于热环境下的加工，热应力主要受两方面影响：一是其内部温度变化，二是其本身的材料参数的变化。通过ANSYS Workbench 进行瞬态热结构耦合不仅能够直观地得到热熔式成孔器对地膜的力学特性，还能观测地膜热应力的变化情况。

学者对薄膜热应力进行了大量研究。何志童等[2]研究了热应力对FeCo 薄膜磁晶各向异性的影响，分析得出较大热应力的引入对于提高磁晶各向异性能有显著的效果；何涛[3]利用有限元方法研究了不锈钢316L/Ti/ta-C 膜系从沉积温度冷却到室温时的热应力分布规律，分析得出不同沉积温度、钛过渡层厚度、不锈钢基底厚度和ta-C 薄膜厚度情况下，ta-C 薄膜和钛过渡层所受热应力均为压应力；吴靓臻等[4]采用有限元模型对SiO2薄膜热应力进行模拟薄膜热应力的大小和分布，分析了不同镀膜温度、不同膜厚和不同基底厚度生长环境下热应力的大小。

从现有研究看，对不同薄膜的热应力研究都比较广泛，但是对地膜的热应力研究还未开始，因此研究地膜热应力在不同情况下的变化，包括不同温度热熔式成孔器、不同高度热熔式成孔器，对于解析地膜的物理性质有重要意义。

本文采集热熔式成孔器以及地膜特性的相关数据，建立热熔式成孔器-地膜的三维模型，设置热熔式成孔器与地膜热结构耦合模型，基于ANSYS Workbench 有限元分析软件对将要打孔的0.01 mm 黑膜进行数值模拟，获得不同情况下的热熔式成孔器对地膜作用时地膜所产生的热应力。为研究地膜热应力提供了一定的参考依据。

1 地膜热应力问题分析

1.1 瞬态温度场有限元求解方程

求解瞬态温度场问题的核心是利用相应的数值方法求解线性常微分方程组，需要通过求解线性常微分方程组，接着通过推导三维温度场得到一个以时间t 为独立变量的线性常微分方程组[5]：

式中：C——热容矩阵；φ——节点温度对时间的导数列阵；K——热传导矩阵；——节点温度列阵；P——温度载荷列阵矩阵。C，K，P 的元素由单元相应的矩阵元素集成[6]，为

通过计算，将时间域和空间域的偏微分方程问题在空间域内转变成N 个节点温度常微分方程的初值问题。

1.2 热-结构耦合分析理论

热应力是由温度场和应力场相互作用的结果，温度场和应力场互相作用为热固耦合问题。针对这些耦合问题可以从物理方程出发，热-结构分析计算的有限元方程为[7]

式中：M——质量矩阵；u——位移；T——温度载荷；C——结构阻尼矩阵；Ci——比热矩阵；K——结构刚度矩阵；Kt——热传导矩阵；F——总等效结点力列阵；Q——总等效结点热流率向量。

1.3 耦合热应力有限元分析

仿真采用顺序耦合的方法进行分析。首先建立模型的结构，通过ANSYS Workbench 自带的模块建立模型，但是由于其操作复杂，可以尝试通过SolidWorks 建立模型导入到ANSYS Workbench；其次定义结构的材料属性，使模型属性符合实际；接着对建立的模型划分合适的网格，网格不合适会影响计算精度[8]。在进行瞬态热分析时，需要添加一系列条件，对2 个结构的初始温度、换热系数以及换热方式进行设置。在得到瞬态温度场分布以后，将得到的温度载荷输入到结构力学分析，利用给定的温度载荷，添加固定支撑条件以及z 轴方向的位移条件，从而完成耦合应力的计算。最后对得到的结果进行分析和检验。

2 热熔式打孔机构对地膜热应力数值计算

2.1 模型建立及网格的划分

在SolidWorks 中建立模型导入到ANSYS Workbench。热熔式成孔器长、宽、高分别为10，30，30 mm，其下半部分的圆柱半径为20 mm。地膜的长、宽、厚分别为20，10，0.01 mm。如图1 所示。

图1 热熔式成孔器-地膜的三维模型Fig.1 Three-dimensional model of hot-melt hole former-mulching film

模型建立以后对模型进行网格划分。网格划分的好坏直接影响到计算的精度和速度。热熔式成孔器与地膜属于不均匀接触过程，因此对地膜处的网格需要设置较密。由于热熔式成孔器是一个不规则形状的模型[9]，需要对其网格进行特殊处理，同时为了提高计算精度且节约计算资源，需要分别划分网格。地膜处的网格采用Face Meshing 的方法，将选中面生成映射网格。Face Meshing 的优点是能够快速生成高质量的六面体网格。由于热熔式成孔器下方呈现圆环状，系统默认采用扫掠网格划分法，但是无法生成映射网格，需要在设置此部分网格时点击mesh—insert—Face Meshing，在Face Meshing中设置Internal Number of Divisions 为6，选中圆环的环面生成网格。最后得到有限元模型中包含节点个数为30 048，单元个数为6 999。网格生成图如图2 所示。

图2 网格生成图Fig.2 Grid generation diagram

2.2 材料热物性参数

温度场是在x，y，z 直角坐标系中连续介质在各个地方同一时刻的温度分布。温度场模拟所需要的材料热物性参数有导热系数λ（W/m·℃）、比热容c（J/kg·℃）、密度ρ（kg/m3）、热膨胀系数α（1/℃）等。热熔成孔器主要由结构钢构成，参考现有的材料参数[10]，其导热系数为44.19 kW/(m·c)，比热容约为600 J/kg·℃，密度7 850 kg/m3，热膨胀系数为1.1×10-5/℃。地膜参考聚乙烯材料参数[11]，导热系数为0.42 kW/(m·c)，比热容约为0.002 J/kg·℃，密度为950 kg/m3，热膨胀系数为0.23×10-5/℃，熔点为105～115℃。将其材料参数在瞬态热分析的Engineering Data 内选定后，进行下一步骤。

2.3 热-结构耦合模块设置

在进行耦合之前，需要对两物体接触进行分析。由于热熔式成孔器与地膜表面相切，并且会沿z 轴方向相对移动，因此两物体属于非线性接触。将模型导入Mechanical 后，系统自动将二者设置为Frictionless。Frictionless 接触类型能够迭代多次且允许法向分离。在Frictionless 设置中还需要选择接触面与目标面，接触面选择热熔式成孔器的圆弧面，目标面选择地膜的上表面。其中接触面显示红色，目标面显示蓝色。为了防止在接触过程中热熔式成孔器穿透地膜表面，需要选择非对称接触选项。生成接触对以后，还需要对接触行为、接触算法等方式进行设置。

完成接触设置后，进行热结构耦合相关设置。仿真分为2 部分：第1 部分为分析热熔成孔器与地膜换热时的温度场；第2 部分为将所求的温度场作为载荷输入到结构力学进行热应力分析。

第1 部分分析结构随时间的变化而变化的瞬态温度场。在瞬态热结构内设置热熔式成孔器模型为恒定热源，表面温度一直维持在设定值，同时设置地膜温度为22 ℃。计算时，对于复杂的分析需要设置一些控制器选项。在Transient Thermal 中的Anslysis Settings 进行热分析设置。对瞬态热的载荷步控制采用时间作为跟踪参数，表示真实时间历程的变量在变化，计算得到的结果是时间的相关函数。设置结束的时间为10 s，初始时步为0.1 s，最大时步为1 s，最小时步为0.001 s，同时需要开启时间积分设置。在求解器的设置上，系统有直接求解器Direct 和迭代求解器Iterative 供选择。此处选择程序控制，即在计算过程中由系统自动选择两者中的一种求解方式，有利于得到较为精确的结果。

第2 部分是将所得的温度载荷输入到瞬态结构力学模块计算得到应力结果。为了约束实体在x，y，z 方向的不规则移动，需要添加固定支撑条件。在对模型进行固定支撑约束时，选取地膜的4 个面作为固定约束。由于热熔成孔器是自上而下移动与地膜接触传热，因此需要添加位移约束，即约束实体在x，y，z 方向的移动，在z 轴方向输入一个具体的值，表示在该方向强制移动此距离，其他2 个方向的位移约束均设置为0，便能模拟出热熔成孔器的运动过程。在Transient 中的Anslysis Settings 进行结构分析设置，载荷步控制同样采用时间作为跟踪参数，设置结束的时间为10 s，初始时步为0.1 s，最大时步为1 s，最小时步为0.001 s，与Transient Thermal 的设置保持一致。在求解器的设置上，与Transient Thermal 设置有所不同，需要将求解器控制中的Weak springs 打开。因为在有限元软件中，经过网格离散化，几何体可能处于微小的不平衡状态，打开后便能使得模型在一个平衡状态下计算。由于地膜的材料属性受热易变形，还需要将Large Deflection 打开。

Workbench 结果后处理功能很强大，能得到用户所需的多种结果。如变形、应力应变、接触工具、疲劳工具等。由于本仿真探求的是关于地膜的应力问题，需要在Solution 工具条Strain 中的Stress 查看应力命令。等效应力/应变又称冯米斯(von-Mises)等效应力/应变，是塑形材料计算时最常用的结果，计算结果需满足第四强度理论σVon≤σs[12]，所以只需查看等效应力的命令便能得到地膜应力的结果。

由于热熔成孔器与地膜接触时具有速度快、接触时间短、加热面积小等特点，为了更直观地观测地膜温度变化，可从地膜表面温度场进行分析。

在热熔成孔器的表面温度分别设置为200，210，220℃，地膜表面温度22℃，位移条件为0.008 m。取2.5，5，7.5，10 s 四个时间研究温度场的分布。通过ANSYS Workbench 后处理分别得到4 个时间段的温度场，如图3 所示。

图3 不同时刻地膜表面温度分布Fig.3 Surface temperature distribution of mulch film at different times

刚开始加热时，热熔成孔器温度还不能及时传递到地膜上，2.5 s 时中心处的温度分布不均匀，直到10 s 时中心处的温度才逐渐呈均匀分布状态，随着时间变化，热熔成孔器与地膜接触时间长，与地膜交互的温度增大。由图3 可知，地膜表面温度随热熔成孔器表面温度同步增大，当热熔成孔器表面温度达到220 ℃时，地膜表面温度达到其熔点。

由于热熔式成孔器与地膜接触传热，导致地膜的温度发生变化，这种变化使得地膜内部产生热应力。将上一步计算得到的温度载荷导入到瞬态结构模块，并结合所给定的固定支撑条件以及位移条件进行地膜的热应力分析。

在热熔成孔器的表面温度设置成为220 ℃，地膜表面温度22℃，位移条件分别为0.004，0.006，0.008 时的热应力分布情况。取2.5，5.0，7.5，10 s 四个时间研究热应力的分布。通过ANSYS Workbench 后处理分别得到4 个时间段的热应力图如图4 所示。

从图4 可以看出，考虑位移约束以后地膜的应力变化随高度增大而变大。采用约束条件以后，可以分析出不同情况下应力的变化情况。在位移初期，由于交互位移变化较小，应力的轮廓面比较不规则。随着时间的增大，位移向下移动与地膜接触面积增大，应力轮廓面呈出一个椭圆形的形状，且边缘较为平整。在移动过程中，由于热熔式成孔器下移高度的变化，其应力也是随之变化的。在下移8 mm时，呈现的应力集中在热熔式成孔器与地膜接触处，这可能造成地膜的破裂。轮廓面上的应力图尺寸符合打孔尺寸。

在热应力的作用下，地膜会发生变形或者破裂。孙伟[13]研究发现地膜在应力8.49e6 Pa 时会发生破裂。在对轮廓面的处理上，设置Capped Isosurface 的值为8.49e6 Pa，并且点击Top 选项便能从轮廓面上看到大于选定值的区域，这样能够直观地看到应力作用于地膜造成破裂的形状及尺寸大小，如图5 所示。

3 结论

通过ANSYS Workbench 将瞬态热分析与瞬态结构力学分析相结合的方法，开展针对地膜在200～220℃内的瞬态温度场仿真，以及针对热熔打孔器在-4，-6，-8 mm 的瞬态力学仿真。结论如下：

（1）热熔式打孔器在工作过程中与地膜接触，由于温度梯度而引起内能交换，会产生受热不均的现象，如果没达到合适温度会造成成孔不完整或未穿透现象。

（2）通过瞬态热分析计算，得到了热熔式成孔器对地膜传热的温度场变化，在其他条件相同的情况下热熔是成孔器达到220℃时，地膜与热熔式成孔器接触处表面最高温度达到其临界温度。温度的变化主要集中在最初的10 s 内，后期的变化较平缓。且在升温过程中，通过对地膜热结构耦合计算，得到了地膜的应力场分布，整个过程中地膜应力场分布均匀，只在两工件接触区域附近应力值较大，其余大部分区域无应力或者应力变化较小。当下降高度为8 mm 时所产生的应力作用得到的孔洞符合尺寸要求。

（3）对地膜受热过程中发生的变形进行分析，引起地膜变形的包括轴向变形和径向变形，由于设置固定支撑及位移约束，轴向变形随高度的变化而改变。通过观察仿真轮廓图地膜中心处的变形区域，对实际生产起到指导作用。