具有损伤演化特性的荸荠粘弹性本构模型的建立与验证

《Postharvest Biology and Technology》：Development and validation of a viscoelastic constitutive model with damage evolution for water chestnut

　　•为莲子建立了一种考虑损伤演变的改进型伯格斯模型。•将损伤因子设定为应力与加载时间的函数。•通过稳定性分析确定冲击能量作为内在损伤指标。•开发了用于粘弹性-损伤耦合仿真的VUMAT子程序。•有限元仿真验证了该模型在预测蠕变损伤和碰撞损伤方面的有效性。

1. 引言
莲子（Eleocharis dulcis）属于莎草科的块茎作物，因其可食用及药用价值而备受重视。它在中国南方以及印度、泰国等亚洲热带和亚热带国家广泛种植。作为重要的特色经济作物，2021年中国莲子出口额达1000万美元（Li等人，2023年）。然而，在机械化采收过程中，块茎需要从土层中挖掘出来并分离，这使其容易受到碰撞和摩擦造成的表面损伤。在后续的运输、分选和储存过程中，莲子还会遭受堆放挤压和振动冲击等各种机械应力，从而导致组织结构受损，增加病原体侵入、腐烂和交叉污染的风险，最终造成批次损失（Zhang，2023年）。因此，研究莲子在机械载荷作用下的变形行为和损伤机制，对于优化采收设备设计、减少产后损失具有重要意义。

本构模型作为描述材料应力-应变关系的数学表达式，为分析农业材料的力学行为提供了理论基础。像莲子这样的农业材料通常具有粘弹性特征，其力学响应既包含弹性固体的可逆变形，又包含粘性流体的时间依赖性行为（Xie等人，2024年）。目前，包括麦克斯韦模型、开尔文模型和伯格斯模型在内的经典流变学模型已被广泛应用于描述苹果（Ji等人，2019年）、番茄（Zhang等人，2018年）和梨（Celik，2017年）等农产品的粘弹性特性。其中，伯格斯模型由串联的麦克斯韦单元和开尔文-沃伊特单元组成，能够同时表征材料的瞬时弹性变形、蠕变变形和永久变形（Sakamoto等人，2007年；Sakamoto等人，2009年；Higashimori等人，2010年）。不过，传统流变学模型假设材料参数恒定，没有考虑载荷作用下材料内部微观结构中的损伤演变对粘弹性行为的影响。Xu等人（2025年）采用了经过指数修正和幂函数修正的伯格斯模型来描述苹果果肉在受损前后的粘弹性特性，结果表明这些修正模型在描述应变饱和现象方面更具优势。

在从生产到加工的整个过程中，莲子等农业材料常常因坠落、碰撞、振动和挤压等机械作用而出现组织损伤或变形，从而降低其品质和商业价值。根据所受外力的特点，农业材料的损伤可分为三类：静态载荷损伤、振动损伤和冲击损伤。静态载荷损伤是指材料在静态条件下长期受压时出现的组织变形或失效，主要发生在堆放储存和包装运输过程中。研究人员通常通过静态压缩试验来模拟加载过程，以变形量、应力-应变关系和组织变化作为评估指标（Jahanbakhshi等人，2020年）。振动损伤则是由于长期处于周期性动态载荷作用下，材料在交变应力作用下逐渐产生疲劳并积累损伤，这一过程具有明显的时间依赖性（Mir和Shahbazi，2022年；Keyhan等人，2024年）。Zhao等人（2024年）基于疲劳累积损伤理论分析了梨的振动损伤特性，确定了其振动疲劳常数和损伤加速阈值。冲击损伤则源于瞬态载荷的短暂高强度作用，作用时间极短，可能导致材料瞬间发生塑性变形或脆性断裂，多发生在采收、分级以及装卸作业过程中（Pathare和Al-Dairi，2021年；Pathare和Al-Dairi，2022年；Hou等人，2023年；Pathare等人，2023年）。在摆锤冲击研究方面，Ye?ilo?lu Cevher和?ztekin（2022年）通过摆锤冲击试验研究了桃子的碰伤形成机制，发现基于冲击能量的损伤体积预测模型精度最高。Li等人（2024年）对蜜桃进行了单摆碰撞试验，获得了包括冲击力和接触时间在内的关键力学参数，为选择损伤评估指标提供了参考。目前，大多数损伤表征仍停留在对比受损前后的水平，缺乏对农业材料损伤演化过程的动态描述。

将数值模拟技术与实验方法相结合，为农业材料损伤研究提供了一种有效途径（Rashvand等人，2022年）。有限元方法已成功应用于苹果（Diels等人，2016年；Ji等人，2019年；Zulkifli等人，2020年；Chen等人，2025年）、猕猴桃（Du等人，2019年）和蓝莓（Hou等人，2024年）等农产品的力学行为分析。Miraei Ashtiani等人（2019年）建立了柚子的双层有限元模型，用于研究其在外部压力作用下的机械损伤敏感性。Ping等人（2023年）使用ANSYS Workbench的显式动力学模块以及双线性各向同性硬化模型，研究了落高、碰撞方式和碰撞物体对马铃薯损伤的影响。Li等人（2024年）将蜜桃视为由果肉和核组成的双层模型，通过碰撞仿真从冲击面提取了最大力、损伤面积、损伤体积和吸收能量等力学参数。然而，现有的仿真研究大多采用有限元软件内置的标准粘弹性材料模型，无法准确描述莲子等农业材料在载荷作用下的损伤演化行为。

综上所述，目前针对农业材料的力学研究存在以下局限性：（1）传统的粘弹性本构模型没有考虑损伤演变对材料参数的影响，难以准确描述材料在载荷作用下的真实力学响应；（2）损伤研究大多停留在阈值标定和定性分析层面，缺乏基于材料本构特性的定量损伤预测方法；（3）有限元仿真使用的标准材料模型无法反映农业材料独特的粘弹性-损伤耦合行为。为解决这些问题，本研究以莲子为研究对象，旨在建立一种能够准确描述其粘弹性行为和损伤演化过程的本构模型，并配套开发相应的有限元仿真方法。主要工作包括：（1）通过不同应力水平下的蠕变和恢复试验分析莲子的粘弹性特性，提出一种考虑损伤演变的改进型伯格斯模型；（2）通过摆锤冲击试验筛选莲子的内在损伤评价指标，并确定其临界阈值；（3）将改进型伯格斯模型扩展到三维空间，用Fortran语言开发Abaqus VUMAT用户材料子程序，以实现莲子粘弹性-损伤耦合行为的有限元仿真；（4）通过蠕变仿真和碰撞仿真验证模型的有效性。研究结果可为莲子机械化加工设备的优化设计提供理论依据。

2. 材料与方法

2.1 试验材料
试验材料选用中国湖北省黄冈市团风县生产的莲子。这些莲子在收获时已完全成熟，购买后立即存放在低温冰箱（0±5°C）中。为保证样本的新鲜度及材料性质的稳定性，所有试验均在48小时内完成。试验前挑选形状大小均匀、表面无机械损伤的莲子样本，其平均横向直径为40–48毫米，厚度为20–27毫米，单个果实质量为20–30克。为尽量减少生物个体差异对试验结果的影响，每组试验均使用同一批次的莲子样本。

2.2 蠕变与恢复试验
蠕变试验采用图1(a)所示的蠕变试验装置进行，该装置主要由高精度数字表头（精度为0.001毫米）、数据采集系统、螺纹砝码和样品支撑平台组成。试验过程中，数据采集系统实时记录莲子在压缩载荷作用下的蠕变变形情况，并将数据传输至计算机进行存储和分析。螺纹砝码可确保静态载荷均匀且稳定地作用于样品表面，而支撑平台则为试验提供稳定的基础，减少环境振动对测量结果的干扰。

考虑到莲子在实际运输和加工过程中通常处于自然状态受力，本研究直接对完整莲子果实进行蠕变加载试验，而非先将其加工成标准几何形状再测试。这种方法更具工程实用价值，能更准确地反映完整材料的真实力学性能。在正式试验之前，为消除微小间隙对试验结果的影响，先让样品承受200克标准砝码的预加载30秒，然后再卸载。当样品不再出现蠕变恢复现象后，开始正式试验。正式试验中，使用图1(b–f)所示的不同质量的砝码，施加F=2.94牛顿、4.90牛顿、6.98牛顿、9.10牛顿和11.53牛顿五种压缩载荷。如图1(a)所示，平头圆柱形压头直径为d=3毫米。由于接触面为平面，且有预加载工序，加之变形量相对于压头直径而言较小，因此在整个加载过程中接触面积基本保持不变。压头下方的名义应力可通过公式σ=4F/πd2计算，得到的数值分别为0.416兆帕、0.693兆帕、0.987兆帕、1.287兆帕和1.631兆帕。蠕变试验的采样频率为5赫兹，总试验时间为180秒，其中包括90秒的加载时间和90秒的卸载恢复时间，以便全面掌握莲子的蠕变和恢复特性。90秒的加载和卸载时间是根据初步试验确定的，在大约90秒后应变率显著下降并进入准稳态阶段。选择这个时间长度是为了在保证实验效率的同时提升准确性。需要指出的是，90秒可能无法完全捕捉更长时间内的延迟弹性恢复现象，但90秒时的残余变形可以较为保守地估算不可逆变形量，这对于工程损伤评估来说是可行的。应变量是通过公式ε=D/l计算的，其中D为蠕变变形量，l为莲子的厚度。为减少个体差异对试验结果的影响，不同的载荷依次作用于同一莲子的不同位置。对于每个载荷等级，随机选取5个莲子样本进行重复试验，最终以平均值作为结果。

2.3 摆锤冲击试验
摆锤冲击损伤试验装置是在指针式摆锤冲击试验机（Yishite，ST-XBL-LID）的基础上改造而成的，以满足莲子碰撞试验的需求，如图2(a)所示。在旋转中心安装了一个角编码器（Witmotion，JY-ME02，角度精度为0.01°），其传动轴通过单盘联轴器与试验台旋转轴相连，以此提高轴的对齐度，减少冲击振动对编码器的影响。摆锤的初始释放角度通过指针和角度表进行校准。在冲击端安装了一个压力传感器（Ruilide，RDF-TM15B，测量范围为0–200牛顿，精度为0.1牛顿），用于监测冲击过程中的碰撞力变化。由于碰撞过程持续时间极短，因此使用了高速数据采集卡（Ruilide，RDC-DAQ，采样频率为5000赫兹）进行数据采集，以确保能够完整记录碰撞过程中的力学响应。为适应莲子的大小和形状，设计了兼具垂直和水平固定功能的专用夹具，用于固定莲子进行冲击试验。

图2. (a) 摆锤冲击试验装置；(b) 不锈钢摆锤；(c) 铝合金摆锤；(d) 铝合金+50克砝码；(e) 铝合金+100克砝码；(f) 铝合金+150克砝码；(g) 损伤直径W1和W2；(h) 损伤深度d。为筛选导致莲子损伤的内在因素并确定其临界阈值，试验采用了五种不同质量的摆锤进行碰撞试验，分别是不锈钢摆锤（425.65克）、铝合金摆锤（123.36克）、铝合金+50克砝码、铝合金+100克砝码以及铝合金+150克砝码，如图2(b–f)所示。铝合金摆锤的两侧对称地安装了配重，以确保摆锤的质心始终位于中心线上。使用SolidWorks制作了摆锤和配重的比例模型，并结合电子天平测量的实际质量（精度为0.01克），校准了每个摆锤的质心到旋转中心的距离：不锈钢摆锤为243.37毫米，铝合金摆锤为243.03毫米，铝合金加50克配重为268.85毫米，铝合金加100克配重为280.81毫米，铝合金加150克配重为288.33毫米，这些数据为计算方程式(1)中的冲击能量和吸收能量等损伤指标提供了依据。测试在多个初始释放角度下进行：5°、10°、15°、20°、25°、30°、35°和40°。每种测试条件重复六次，取平均值作为最终结果。通过最大变形角、反弹角、最大冲击力、冲击速度、冲击能量、吸收能量以及恢复系数等指标来评估水栗的损伤程度。最大变形角乘以摆锤旋转中心到冲击中心的距离，可换算出水栗在冲击过程中的最大变形量。反弹角是指摆锤反弹过程中达到的最大角度。最大冲击力则是冲击过程中压力传感器记录的最大值。为避免空气阻力和摩擦造成的能量损耗影响测试结果，采用角度编码器数据来计算冲击速度和能量指标：(1)v=?θ/Δt，Q1=1/2mv12，Q2=1/2mv12?1/2mv22，e=v2/v1，其中v为冲击点的最大速度，单位为米/秒；l为旋转中心到冲击点的距离，单位为米；v1和v2分别为摆锤质心的最大速度和最大反弹速度，单位均为米/秒；m为摆锤的质量，单位为千克；Q1为冲击能量，单位为焦耳；Q2为吸收能量，单位为焦耳；e为恢复系数。水栗的表皮呈深棕色，因此冲击后表面难以直接观察到淤伤损伤。为此采用了切片法对损伤体积进行定量分析（Opara和Pathare，2014年）。碰撞处理后，使用游标卡尺测量冲击点处的表面压痕尺寸W1和W2。随后将水栗置于室内保存24小时，再切开受损区域，以变白深度作为损伤深度d，如图2(g–h)所示。实验观察表明，表面淤伤大致呈椭圆形，内部的变白组织沿冲击方向向内延伸。基于这种形态特征，损伤体积Vd的计算公式为：(2)Vd=πd2/43W1W2+4d2，该公式常用于估算表面损伤呈椭圆形且向内发展的园艺作物的淤伤体积。2.4 粘弹性损伤本构模型Burgers模型由一个麦克斯韦单元（弹性单元E1与粘性单元η1串联）和一个开尔文-沃伊特单元（弹性单元E2与粘性单元η2并联）串联组成。它能描述三种典型的粘弹性行为：瞬时弹性变形、延迟弹性变形以及粘性流动，因此是表征材料粘弹性特性的经典模型。其加载阶段和卸载阶段的一维本构方程分别为：(3)ε(t)load=σE1+σE2/(1?e^?t/τ)+ση1t/εtunload=εt/(1?ση1t/(1?e^?t?t1/τ))+ση1t1，其中ε(t)为应变；σ为应力，单位为兆帕；E1为表征材料瞬时弹性响应的瞬时弹性模量，单位为兆帕；E2为表征材料延迟弹性响应的延迟弹性模量，单位为兆帕；η1为表征材料不可逆粘性流动的粘度系数，单位为兆帕·秒；η2为延迟粘度系数，单位为兆帕·秒；τ=η2/E2为迟滞时间，单位为秒；t和t1分别为蠕变时间和卸载时间，单位均为秒。在方程式(3)中，加载阶段本构方程的第一项代表瞬时弹性变形，这种变形在加载后立即发生，卸载时也会立即恢复。第二项代表延迟弹性变形，它会随时间逐渐增大，卸载后也会逐步恢复。第三项代表粘性流动，它随时间线性增加，卸载后无法恢复。本研究将基于蠕变试验数据对Burgers模型进行参数识别，并根据实验观察结果对模型进行必要修改（详细内容见3.1节）。2.5 有限元仿真方法使用手持式3D扫描仪对去掉了顶芽的整个水栗果实进行多方向扫描（图3(a)），以获取完整的点云数据。对扫描得到的模型进行表面轮廓去噪和平滑处理，同时减少三角形面的数量（图3(b)）。之后使用C3D8R八节点线性六面体简化积分元素对优化后的3D模型进行网格划分，整体网格尺寸为0.5毫米，共生成199388个元素（图3(c)）。为验证网格的合理性，采用了四种不同的全局网格尺寸（1.5毫米、1.0毫米、0.5毫米和0.3毫米）进行了网格收敛性分析。如图3(d)所示，0.5毫米和0.3毫米网格尺寸下的变形-时间曲线几乎重合，最大变形量的相对差异仅为0.7%，这表明所选的网格尺寸能够实现良好的收敛效果。下载：下载高分辨率图像（295KB）下载：下载全尺寸图像图3. 水栗几何模型的构建过程。(a) 3D扫描；(b) 模型优化；(c) 网格生成；(d) 网格收敛性分析；(e) 有限元模型。为在有限元软件中实现所建立的粘弹性损伤本构模型，用Fortran语言开发了一个Abaqus VUMAT用户材料子程序。该子程序在假设材料各向同性的基础上，将一维本构方程扩展为三维形式，通过应力与应变的体积偏微分分解，将弹性模量转换为剪切模量，将体模量转换为剪切模量，以此分别描述剪切变形和体积变形。三维本构方程的推导及子程序的实现细节见3.3节。蠕变模拟是使用Abaqus/Explicit求解器进行的，共包含四个分析步骤：初始步骤（1秒，模拟瞬时加载响应）、蠕变步骤（90秒，模拟压缩蠕变行为）、卸载步骤（1秒，模拟瞬时卸载响应）以及恢复步骤（90秒，模拟蠕变恢复行为）。边界条件的设置如下：底部表面的垂直位移被限制，而侧面则可以自由膨胀；顶部表面则施加与实验中压头相同接触面积的局部分布应力载荷。碰撞模拟也是通过Abaqus/Explicit求解器来进行的。水栗的材料参数与蠕变模拟中所使用的参数一致。摆锤则根据其实际的几何形状和质量进行建模，被定义为离散的刚体（R3D4元素），并置于相应的初始释放角度位置。同时施加重力场，使摆锤能够自由下垂并撞击水栗，从而再现实际的实验过程。接触属性设置为：法向为硬接触，切向则采用摩擦系数为0.6的惩罚摩擦算法。碰撞有限元仿真模型如图3(e)所示。模拟完成后，通过Abaqus的后处理功能，根据吸收能量密度分布来确定损伤区域。超过临界阈值的元素会被标记为受损元素，将其体积相加即可得到预测的损伤体积。确定损伤阈值以及进行损伤体积定量分析的方法详见3.2节“内在损伤指标的筛选与阈值确定”以及3.3节“有限元仿真结果与分析”。3. 结果与讨论3.1 粘弹性特性分析及本构模型建立3.1.1 蠕变试验数据分析在9.10牛顿的载荷作用下，水栗的蠕变与恢复曲线如图4(a)所示。在加载瞬间，会出现瞬时弹性变形；在保持载荷的阶段，变形会随时间逐渐增大，呈现出典型的蠕变行为；卸载时则会出现瞬时弹性恢复，随后还会继续恢复，逐渐趋于稳定状态。值得注意的是，试验结束时仍存在明显的不可逆变形。为进一步验证这一现象，进行了时长为30分钟加载和30分钟恢复的长期蠕变试验（图4(b)）。结果表明，经过长时间蠕变后，不可逆变形显著增加，这证实了水栗具有典型的粘弹性特征。下载：下载高分辨率图像（259KB）下载：下载全尺寸图像图4. 水栗蠕变特性的分析。(a) 在9.10牛顿载荷作用下的短期蠕变曲线（180秒）（对应应力为1.287兆帕）；(b) 长期蠕变曲线（60分钟）；(c) 标准Burgers模型及其理论曲线；(d) 水栗的实际蠕变曲线及修正后的Burgers模型。尽管水栗的蠕变行为总体上符合Burgers模型的预测趋势，但详细分析发现，加载时的瞬时弹性变形与卸载时的瞬时弹性恢复之间存在显著差异（图4(d)），而标准Burgers模型则认为这两者应该相等（图4(c)）。出现这种现象的原因是，水栗在静态压缩载荷作用下会遭受不可逆损伤，从而导致材料参数发生变化，最终使得卸载结束时的不可逆变形程度高于标准Burgers模型的预测值。因此，有必要对标准Burgers模型进行修正，以便更准确地描述水栗的粘弹性特性。3.1.2 损伤因子的确定通过对蠕变参数对曲线形态影响的分析发现，瞬时变形由瞬时弹性模量E1决定，而卸载时的塑性变形则由粘度系数η1决定。对不同载荷水平下的蠕变试验数据进行分析后发现，加载阶段和卸载阶段的延迟弹性变形大致相等，这说明在蠕变过程中，延迟弹性模量E2和延迟粘度系数η2保持稳定。因此，后续的损伤修正仅聚焦于瞬时弹性模量E1和粘度系数η1。根据Burgers模型的物理意义，从加载阶段和卸载阶段的蠕变曲线中分别提取出E1和η1。加载时的瞬时弹性模量是根据加载时的瞬时应变计算得出的：E1=σ/ε0。卸载时的瞬时弹性模量则是根据卸载时的弹性恢复情况计算得出的：E1′=σ/?εe，其中?εe为卸载时的弹性恢复应变。加载时的粘度系数是通过对方程式蠕变应变曲线进行线性拟合，根据稳态应变率计算得出的：η1=σ/k，其中k为稳态蠕变阶段的应变率。卸载时的粘度系数则是根据试验结束时的塑性应变计算得出的：η1′=σ·t1/εp，其中εp为最终的不可逆塑性应变，t1为蠕变时间。图5(a)和(b)展示了加载阶段和卸载阶段各自的瞬时弹性模量与粘度系数，以及它们随载荷水平的变化情况。结果表明，水栗在蠕变过程中会发生压缩损伤，这使得卸载时的E1大于加载时的值（瞬时弹性恢复程度降低），而卸载时的η1则小于加载时的值（塑性变形程度增加）。此外，加载阶段与卸载阶段之间的E1和η1差异都会随着外部载荷的增加而增大，这说明损伤程度与蠕变应力呈正相关关系。下载：下载高分辨率图像（250KB）下载：下载全尺寸图像图5. 蠕变参数的变化及损伤因子拟合曲线。(a) 瞬时弹性模量E1随应力的变化；(b) 粘度系数η1随应力的变化；(c) DE1损伤因子拟合曲线；(d) Dη1损伤因子拟合曲线。基于上述实验结果，引入了损伤因子来描述瞬时弹性模量E1和粘度系数η1随蠕变应力的变化情况。修正后的Burgers模型如图4(d)所示，与标准模型不同的是，弹性单元E1和粘性单元η1不再是常数，而是与蠕变应力和时间相关的变量。根据E1增大而η1减小的趋势，损伤因子的定义如下：(4)E1′=E1/(1?DE1)，η1′=η1/(1?Dη1)，其中E1′为受损后的瞬时弹性模量，单位为兆帕；η1′为受损后的粘度系数，单位为兆帕·秒；DE1和Dη1分别为瞬时弹性模量和粘度系数的损伤因子；E1和η1则为初始未受损时的数值。根据方程式(4)，计算了不同载荷水平下的损伤因子。对于每一种应力水平，都在五颗随机选取的水栗样本上重复进行蠕变试验，然后在提取参数之前先对所得的蠕变-恢复曲线取平均值。表1显示，对于E1损伤因子，幂函数拟合效果最佳（R2=1.000），而线性函数的R2值均高于0.99。对于η1损伤因子，线性函数拟合效果最好（R2=0.777）。综合考虑拟合精度和公式简洁性，两种参数均选择线性函数：DE1=(0.1776σ+0.0717)，Dη1=(0.3051σ?0.0497)。表2展示了粘度系数损伤因子Dη1的拟合结果。3.1.3 损伤演化方程的建立与验证水葫芦在蠕变过程中的损伤不仅与外部应力有关，还是随时间逐渐积累的结果。为使本构模型更符合实际的损伤演化过程，采用指数饱和函数来描述损伤因子随时间的变化规律，该函数能够体现初期快速增长、后期趋于稳定的特征。以不同载荷水平下测试结果的平均值作为拟合依据，得到了损伤因子演化方程：DE1(t)=(0.1776σ+0.0717)×(1?e?t/100)，Dη1(t)=(0.3051σ?0.0497)×(1?e?t/45)，其中DE1(t)和Dη1(t)为随时间变化的损伤因子，σ为蠕变应力，单位为MPa；t为加载时间，单位为秒。基于上述分析，将损伤演化方程代入标准Burgers模型，得到了水葫芦的粘弹性本构模型：εtload=σE101?(0.1776σ+0.0717)×(1?e?t/100)+σE21?e?t/τ+ση1(t)t，εtunload=εt1?ση101?0.3051σ?0.0497×(1?e?t/45)t1e?t?t1/τ+ση1(t1)，其中εt为应变，εt1为卸载时的应变；E10=113MPa，η10=35000MPa·s分别为初始未受损时的瞬时弹性模量和粘度系数；τ=η2/E2为迟滞时间；E2=400MPa，η2=4200MPa·s为延迟弹性模量和延迟粘度系数；t和t1分别为蠕变时间和卸载时间。基于该本构模型，计算了不同载荷水平下的理论蠕变曲线，并将结果与实验数据进行了对比，见图6(a)。在0.416、0.693、0.987、1.287和1.631 MPa的应力水平下，该模型的决定系数R2分别为0.876、0.986、0.964、0.976和0.990，表明该本构模型能够准确描述水葫芦在不同载荷水平下的蠕变行为。为进一步验证该本构模型的普适性，在5.54 N的载荷水平（对应应力为0.784 MPa）下进行了额外的蠕变试验，共重复5次并取平均值。实验结果与模型预测的对比见图6(b)。对于瞬时应变、最大应变和塑性应变这三个指标，模型预测值与实验值之间的相对误差分别为1.16%、1.87%和7.04%，这表明所建立的粘弹性本构模型能够准确描述水葫芦在不同载荷水平下的蠕变行为。下载：下载高分辨率图像（239KB）下载：下载全尺寸图像图6。(a) 不同应力水平下本构模型理论值与实验结果的对比；(b) 5.54 N载荷水平下蠕变实验结果与本构模型理论预测的对比。3.2. 内在损伤指标的筛选与阈值确定3.2.1 不同摆锤的碰撞损伤试验结果采用不锈钢摆锤和铝合金摆锤对水葫芦进行了碰撞损伤试验。使用不锈钢摆锤时，不同初始释放角度下水葫芦的碰撞试验结果见表3。试验中发现，当摆锤初始释放角度为15°时开始出现冲击压痕（见图7(a)），在20°时压痕变得极为明显（见图7(b)），而在30°时则出现直接破裂（见图7(c)）。经过24小时室内储存后再进行切割观察，发现15°时已开始出现白色损伤，这是由于水葫芦的组织结构（细胞骨架）遭到破坏导致局部水分迅速流失所致。因此，确定在水葫芦受到不锈钢摆锤碰撞时的损伤阈值为15°。在所有测试的冲击角度中，表面瘀伤大致呈椭圆形，内部变白组织则沿冲击方向向内延伸（见图7），这与用于计算实验中瘀伤体积的方程(2)的几何假设一致。表3展示了使用不锈钢摆锤的碰撞试验结果。释放角度最大变形角度（°）反弹角度（°）最大冲击力（N）冲击速度（m/s）冲击能量（J）吸收能量（J）恢复系数5°0.1221.8402.8700.1780.00650.00460.54010°0.2253.1574.3260.3390.02340.01870.44615°0.3304.2935.6000.5090.05280.04210.45020°0.6333.8137.7800.6600.08890.08310.29225°0.7185.7228.0830.8240.13880.12480.31830°0.9026.00310.1900.9850.19790.18030.298下载：下载高分辨率图像（289KB）下载：下载全尺寸图像图7. 不同释放角度下水葫芦受到不锈钢摆锤撞击后的表面及截面损伤情况：(a)15°，(b)20°，(c)30°。铝合金摆锤碰撞试验的结果见表4。分析表明，在25°的初始释放角度下，只有少数水葫芦出现轻微的白色损伤，而在30°时则有80%的水葫芦出现损伤。因此，确定在水葫芦受到铝合金摆锤碰撞时的损伤阈值为30°。表4展示了使用铝合金摆锤的碰撞试验结果。释放角度最大变形角度（°）反弹角度（°）最大冲击力（N）冲击速度（m/s）冲击能量（J）吸收能量（J）恢复系数15°0.2984.3022.0100.5000.01450.01200.41220°0.3624.9932.2190.6600.02530.02210.35425°0.5634.7272.5590.7790.03530.03060.36430°0.6175.5133.2950.9490.05230.04700.31835°0.7937.4053.4391.0990.07020.06160.35040°0.8607.7824.3241.2420.08960.08080.314对比两种摆锤的试验结果可知，尽管它们的损伤阈值角度不同（不锈钢为15°，铝合金为30°），但在各自阈值角度下的冲击能量却非常接近（分别为0.0528 J和0.0523 J），这说明冲击能量可能是引发水葫芦损伤的内在指标。3.2.2 内在损伤指标的确定为了筛选出能够引发水葫芦损伤的内在指标，又进行了补充碰撞试验，在铝合金摆锤的两侧对称地悬挂了50克、100克和150克的标准重物。试验结果表明，当悬挂50克重物时，所有水葫芦在23°的释放角度下都会出现白色损伤；悬挂100克重物时，大多数水葫芦在18°的释放角度下就会出现明显损伤；而悬挂150克重物时，有80%的水葫芦在16°的释放角度下出现白色损伤。考虑到水葫芦作为一种生物材料存在个体差异，本研究以大多数样本出现损伤的角度作为每种摆锤配置的损伤阈值。表5汇总了五种摆锤配置在各自损伤阈值条件下的碰撞参数。表5. 不同摆锤在临界损伤阈值下的碰撞参数。摆锤类型损伤指标损伤角度（°）变形角度（°）反弹角度（°）冲击力（N）冲击速度（m/s）冲击能量（J）吸收能量（J）恢复系数铝合金300.6175.5133.2950.9490.05230.04700.318不锈钢150.3304.2935.6000.5090.05280.04210.450铝合金+50克230.6026.9124.2780.7240.05640.04430.463铝合金+100克180.3905.1775.0550.5630.04980.03720.504铝合金+150克160.5124.0675.7080.5130.05450.0449.420表6. 损伤体积预测对吸收能量密度阈值的敏感性分析。释放角度吸收能量密度阈值1.0 mJ/mm31.2 mJ/mm31.386 mJ/mm31.5 mJ/mm31.8 mJ/mm330°218.37%99.86%19.59%?63.91%?100.00%40°5.09%?10.82%?26.00%?34.30%?53.03%50°?3.19%?15.73%?24.63%?31.32%?45.60%60°78.40%47.13%26.89%14.00%?14.64%70°78.07%48.70%28.27%16.57%?9.52%平均绝对误差76.62%44.45%25.08%32.02%44.56%为定量评估不同指标在筛选内在损伤因子时的稳定性，设计了一个指标稳定性评估公式：vardev=varmax?varminvaraverage，其中varmax、varmin和varaverage分别表示不同摆锤条件下同一评估指标的最大值、最小值和平均值。vardev值越小，说明该指标的变异幅度越小，稳定性越好。计算结果显示：变形角度的变异幅度为0.586，反弹角度为0.548，冲击力为0.504，冲击速度为0.675，冲击能量为0.124，吸收能量为0.227。其中冲击能量的相对变异幅度最小，因此确定冲击能量为引发水葫芦损伤的内在评价指标。由于水葫芦是一种生物材料，其损伤阈值会在一定范围内波动。综合五种摆锤类型的测试结果，确定水葫芦的冲击能量损伤阈值范围为0.0498–0.0564 J。3.3 有限元仿真结果与分析3.3.1 本构模型的三维扩展与VUMAT实现在开发VUMAT子程序时，需要将一维本构方程扩展为三维形式。修正后的Burgers模型由串联的Maxwell模型和Kelvin-Voigt模型组成，总应变可分解为三部分：ε=ε?+ε??+ε??，其中ε?为弹性应变（Maxwell弹簧E1），ε??为粘性流动应变（Maxwell阻尼器η1），ε??为Kelvin-Voigt组分应变。进行三维扩展时需要考虑以下关键问题：(1)根据泊松比将杨氏模量转换为剪切模量和体模量；(2)将应力和应变分解为体积分量和偏量分量；(3)假设材料具有各向同性。由于难以在整体水果压缩过程中捕捉完整水葫芦的尺寸变化，因此通过肉样测量泊松比。将边长为10毫米的立方体样品置于TMS-Pro纹理分析仪上进行单轴压缩试验。泊松比的计算公式为μ=?X/?L，其中?X和?L分别为横向和轴向的变形量。该试验重复进行了10次，平均泊松比为0.368，该数值随后被用于三维本构模型中。杨氏模量、剪切模量与体模量之间的关系为：G=E21+νK=E31?2ν。应力和应变的分解公式为：p=13σ11+σ22+σ33，Sij=σij?pδij，εvol=ε11+ε22+ε33，eij=εij?εvol3δij，其中p为静水压力，Sij为偏量应力张量，εvol为体积应变，eij为偏量应变张量，δij为克罗内克德尔塔函数。引入损伤演化方程后，即可得到包含损伤的线弹性应力-应变关系式：σij=λ1εkkδij+2G1εij，其中λ1=K1?2G13，G1=11321+ν1?DE1，K1=11331?2ν1?DE1。包含损伤的Maxwell粘性阻尼器的三维扩展方程为：ε˙iipl=σij?p2×35000×1?Dη1i=1,2,3，ε˙ijpl=σij2×35000×1?Dη1i≠j。Kelvin-Voigt模型的三维扩展方程为：ε˙iikv=σii2η2?G2η2εiikv(i=1,2,3)，ε˙ijkv=σij2η2?G2η2εijkvi≠j，其中G2=E221+ν=40021+ν为Kelvin-Voigt组分的剪切模量。VUMAT中损伤因子的数值更新公式为：DE,incr=0.1776σeq+0.0717(e?t?dt100?e?t100)，Dη,incr=0.3051σeq?0.0497(e?t?dt45?e?t45)，DE,new=minDmax,DE,old+max0,DE,incr，Dη,new=min?Dmax,Dη,old+max0,Dη,incr，其中σeq为等效应力，dt为时间步长，Dmax为损伤的上限值。VUMAT子程序的实现流程图如图8所示。下载：下载高分辨率图像（188KB）下载：下载全尺寸图像图8. VUMAT子程序实现流程图。3.3.2 蠕变有限元仿真结果为验证VUMAT子程序的正确性，首先构建了一个单元素立方体模型（C3D8R单元）进行数值验证。如图9(a)所示，对该模型的底部表面施加固定约束，顶部表面则施加均匀分布的应力载荷。在与蠕变试验对应的五个应力水平下进行了仿真，结果如图9(b)所示。仿真曲线与实验数据总体吻合良好，蠕变应变的最大相对误差为12.19%，出现在1.631 MPa的应力水平下，这证明了该本构模型和VUMAT程序的有效性。下载：下载高分辨率图像（420KB）下载：下载全尺寸图像图9. 蠕变有限元仿真验证。(a) 单元素验证模型；(b) 单元素模型在不同应力水平下仿真结果与实验结果的对比；(c) 1.287 MPa载荷作用下三维水葫芦模型的位移等值线图（上图：加载90秒后的最大变形；下图：卸载90秒后的残余变形）；(d) 三维水葫芦模型在不同应力水平下仿真结果与实验结果的对比。在此基础上，进一步验证了VUMAT子程序在三维几何模型中的应用效果。通过对三维水葫芦模型施加应力载荷，进行了蠕变仿真。图9(c)展示了在1.287 MPa载荷下水葫芦加载90秒及卸载90秒后的位移等值线情况。如图所示，水葫芦的受载区域出现了明显的蠕变变形，卸载后仍存在一定的残余变形。图9(d)展示了在五个应力水平（0.416、0.693、0.987、1.287和1.631 MPa）下，仿真结果与实验结果的变形-时间曲线对比，蠕变变形的决定系数R2分别为0.93、0.91、0.84、0.85和0.76。这些结果表明，所建立的粘弹性损伤本构模型不仅能揭示受载区域的变形集中现象和恢复行为，还能准确描述水葫芦在不同应力水平下的三维粘弹性响应。3.3.3 碰撞有限元仿真与损伤体积预测以铝合金摆锤以50°的初始释放角度撞击水葫芦为例进行了有限元仿真。水葫芦的吸收能量等值线图如图10(a)所示。仿真得到的最大变形量为3.909毫米，与实验测量的5.220毫米相比，相对误差为25.11%，这说明所建立的材质模型能够较好地描述水葫芦在碰撞过程中的变形行为。通过能量分析进一步验证了该模型的可靠性：沙漏形能量是由数值误差或约束不完善导致的虚拟能量，其与内能的比值通常应控制在5%–10%之间（Wallmeier等人，2015年）。在该模拟中，沙漏形能量与内能的比值仅为0.114%，远低于允许范围，从而证明了碰撞模型的数值稳定性。下载：下载高分辨率图像（308KB）下载：下载全尺寸图像图10. 莲子的碰撞有限元模拟及损伤体积计算。(a) 初始释放角度为50°时莲子吸收的能量等值线；(b) 不同初始释放角度下模拟得到的损伤体积与实验结果的对比。在确认碰撞模拟的可行性后，进一步对损伤体积进行了定量预测。根据3.2节中得出的内在损伤指数筛选结果，冲击能量具有最高的评估价值；不过由于吸收能量可直接通过有限元模拟计算，因此在冲击能量之后评估价值第二的吸收能量被选为定量损伤指标。吸收能量密度阈值是通过计算每个摆体在各自临界损伤状态下的吸收能量与损伤体积的比值，再对五个摆体的结果取平均值得出的，其阈值为1.386 mJ/mm3。VUMAT子程序根据变形能和耗散能计算出莲子的吸收能量，并将其输出到SENER变量中（即积分点应变能量密度）。根据SENER值，那些超过临界阈值（1.386 mJ/mm3）的单元会被标记为受损单元。图10(b)展示了不同初始释放角度下莲子受损单元的分布情况以及损伤体积的预测结果。与实验测量值相比，30°、40°、50°、60°和70°这五种条件下的绝对相对误差分别为19.59%、26.00%、24.63%、26.89%和28.27%。各冲击强度下的预测误差相对均匀，平均相对误差为25.08%。此外，模拟得到的损伤区域呈现出近似椭圆形的接触区，且受损单元会向内部扩展（见图10），这与图7中的实验损伤形态相似，进一步证明了使用公式(2)进行损伤体积比较的合理性。对损伤体积预测性能的分析表明，有限元模拟的整体预测精度是可靠的。预测误差主要源于两个方面：首先，作为生物材料，莲子在个体间的几何特征和材料属性存在显著差异；其次，本研究中所建立的粘弹性损伤模型是基于损伤演化方程的，它在线性损伤阶段能够实现高精度预测，但当出现严重的非线性损伤（如大变形或破裂）时，预测精度会下降。为了研究损伤体积预测对阈值选择的敏感性，研究者测试了五种不同的吸收能量密度阈值（1.0、1.2、1.386、1.5和1.8 mJ/mm3）。结果表明，所采用的1.386 mJ/mm3阈值在所有测试值中对应的平均绝对相对误差最低，为25.08%。当阈值偏离这一数值时，预测误差会双向上升：过低的阈值会导致损伤体积被高估，而过高的阈值则会导致其被低估。这一分析证明，所选的阈值对于损伤体积预测而言是经过合理校准的。总体而言，利用本研究建立的损伤演化方程和能量密度阈值来预测碰撞作用下的莲子损伤程度，能够准确量化不同碰撞条件下的损伤严重性。3.4 讨论本研究建立的粘弹性损伤本构模型与现有研究相比具有若干独特之处。传统的流变学模型，如标准Burgers模型，假设材料参数恒定，因此无法反映材料微观结构在持续受力作用下的损伤演化过程。Xu等人（2025）采用了经指数修正和幂函数修正的Burgers模型，用于描述苹果果肉在遭受碰伤前的粘弹性特性。不过他们的模型参数是为未受损状态和受损状态分别校准的，并未建立损伤因子随应力和时间持续演化的数学表达式。相比之下，本研究提出的修正版Burgers模型引入了两个损伤因子DE1和Dη1，能够基于应力水平和加载时间连续动态地描述材料参数的演化过程，从而更准确地呈现农业材料在实际加载过程中的渐进性损伤累积行为。在损伤评估方面，现有研究大多采用落差高度、碰撞速度或碰撞力作为损伤指标。然而这些指标依赖于具体测试装置，无法作为材料的固有属性。Ye?ilo?lu Cevher和?ztekin（2022）在五种不同的冲击能量水平（0.015–0.160 J）下对桃子进行了测试，并建立了多项线性回归模型用于预测碰伤体积。他们的研究结果表明，以冲击能量作为主要自变量的模型（R2 = 0.76–0.78）的预测效果优于以最大接触力作为自变量的模型（R2 = 0.63–0.65），这一结论证实了冲击能量是预测水果碰伤损伤更为可靠的指标。与这一发现一致，本研究通过使用五种不同质量的摆体进行系统性的指标稳定性分析，确定冲击能量是莲子最稳定的固有损伤指标，其临界阈值范围为0.0498–0.0564 J。这一阈值可以与其他农产品的已有数据相对照：Stropek和Go?acki（2020）报告称，梨子在出现碰伤时的内部损伤能量约为0.019–0.023 J，这一数值低于莲子的损伤能量。莲子更高的损伤阈值与其相比梨子和桃子这类软质水果更坚硬的质地以及更为致密的组织结构有关，这也凸显出损伤阈值具有物种差异性，以及制定针对特定作物的损伤评判标准的必要性。有限元模拟结果显示，在不同的加载条件下，该模型的预测精度相当不错。在蠕变模拟中，单元素模型与实验数据的吻合度较好（最大相对误差为12.19%），而三维模型的平均R2值为0.86。单元素模型与三维模型性能之间的差异，主要源于从均匀的单轴应力状态（与校准条件一致）到直径较小的平面压头作用于弯曲的莲子表面所产生的复杂三维应力场的转变。在碰撞模拟中，损伤体积预测的平均相对误差为25.08%。这些误差与现有关于农产品损伤的有限元研究中的报告值相当；例如，Li等人（2024）在使用弹塑性材料模型对蜜桃进行碰撞模拟时，报告的损伤体积最大误差为21.8%，最大力值的最大误差为25.1%。预测误差在所有冲击强度下都相对均匀，这说明该模型在极端条件下的性能依然稳定，没有出现明显下降。这些误差主要源于莲子作为生物材料本身存在的个体差异、将准静态本构参数外推到动态碰撞条件中，以及忽略了皮层与髓部结构差异而做出的各向同性材料假设。需要指出的是，模拟结果与实验结果在局部接触力响应方面仍存在差异：模拟得到的峰值力更高，接触持续时间更短，而冲量方面的吻合度则相对较好。这种差异主要源于将准静态本构参数外推到了动态碰撞条件中，因为当前模型无法完全捕捉莲子组织随应力变化而产生的硬化效应和耗散效应。值得注意的是，本研究中的碰伤预测是基于吸收能量密度而非峰值力来进行的。尽管局部力特性存在差异，但预测得到的损伤体积与实验结果的平均相对误差仍然为25.08%，这表明基于能量密度的损伤预测方法具有较高的可靠性。因此，未来的工作将致力于开发考虑应力变化率的本构模型，以提高动态力-时间响应的预测精度。本研究也存在一些局限性，需要予以承认。首先，测试样本仅来自湖北省团风地区的某个特定品种，没有考虑不同采收期、生长环境以及品种类型对莲子力学性能的影响。其次，目前的本构模型尚未将储存温度和储存时间等环境因素作为内在参数纳入考虑。初步实验表明，莲子的蠕变变形会随着储存温度的升高和储存时间的延长而显著增加；未来的研究可以将这些因素引入本构模型中，从而更全面地反映农产品在实际采后条件下的力学行为。第三，本研究将莲子建模为各向同性的均匀体，忽略了其皮层与髓部之间的结构差异。虽然这种简化方式在农业材料的相关有限元研究中较为常见（Celik, 2017; Du等人，2019; Li等人，2024），并且在适度冲击条件下能够提供足够的预测精度，但在极端载荷作用下，尤其是处于破裂阶段时，它可能会导致预测精度下降。第四，网格收敛性分析表明，所选的网格尺寸（0.5 mm）已经达到了令人满意的收敛效果（与0.3 mm相比相对误差仅为0.7%），这说明网格尺寸并非导致预测误差的重要因素。最后，本构模型参数是在准静态蠕变条件下校准的，将其外推到动态碰撞载荷条件下会引入额外的不确定性，这一点从模拟结果与实验结果在接触力特性上的差异中可以看出。这些局限性为未来的研究指明了方向。由于观察到了应变率对材料响应的影响，未来有必要开发考虑应力变化率的本构模型，以便更准确地描述材料在准静态和动态载荷作用下的响应特性。引入能够区分皮层与髓部的多层材料模型，也有助于提高在极端载荷条件下的预测精度。此外，本研究建立的粘弹性损伤本构模型还可以通过与离散元方法相结合的形式，通过DEM-FEM耦合仿真技术，实现对莲子在真实采收场景中的损伤预测与分析，比如使用弹簧滚筒分离莲藕根茎的过程。这样的耦合方式还能够为整合上述各种改进措施——包括考虑应力变化率、环境因素以及材料异质性——形成一个完善的仿真工具，用于优化相关机械化加工设备的设计。4 结论本研究通过结合实验与仿真方法，建立了考虑损伤演化的粘弹性本构模型，并相应地开发了有限元仿真方法，以此解决莲子在机械化采收、运输和加工过程中的损伤问题。主要结论如下：(1) 通过在不同载荷水平下的蠕变和恢复试验，发现莲子在受载时的瞬时弹性变形与卸载时的瞬时弹性恢复之间存在显著差异，这表明标准Burgers模型无法准确描述莲子在受载过程中的损伤演化行为。通过分析蠕变参数的变化规律，研究者将损伤因子DE1和Dη1分别引入瞬时弹性模量E1和粘度系数η1中，进而建立了损伤因子随应力和时间演化的数学模型。在额外增加5.54 N载荷（对应应力为0.784 MPa）的验证试验中，所建立的修正版Burgers模型在0.416–1.631 MPa的应力范围内与实验结果吻合良好，瞬时应变、最大应变和塑性应变的相对误差分别为1.16%、1.87%和7.04%。(2) 通过使用五种不同质量的摆体进行冲击损伤试验，研究者采用了指标稳定性评估方法对莲子损伤的内在评价指标进行了筛选。结果表明，冲击能量的相对变化幅度最小，因此确定冲击能量是引发莲子损伤的内在指标，其临界阈值范围为0.0498–0.0564 J。(3) 该修正版Burgers模型被扩展到了三维形式，同时研究人员还用Fortran语言开发了Abaqus VUMAT用户材料子程序，从而实现了对莲子粘弹性-损伤耦合行为的有限元仿真。蠕变模拟结果表明，该本构模型能够准确捕捉莲子的粘弹性变形特性，在不同应力水平下的平均拟合优度为0.86。在碰撞模拟中，研究者基于吸收能量密度（1.386 mJ/mm3）对损伤体积进行了定量预测，不同碰撞强度下的平均相对误差为25.08%，这一结果证明了使用有限元方法对莲子的损伤进行定量分析是可行的。本研究建立的粘弹性损伤本构模型及有限元仿真方法，可为莲子及类似块茎类作物的机械化加工设备的优化设计及损伤控制提供理论依据和技术支持。CRediT作者贡献说明Mingli Geng：数据整理。Guozhong Zhang：写作——审阅与编辑、研究指导、资金获取。Peng Wang：写作——审阅与编辑、方法设计、资金获取、概念构思。Xuejin Nie：结果验证。Xuefeng Tang：初稿撰写、软件应用、实验研究、正式分析。Xuefeng Tang|Peng Wang|Xuejin Nie|Mingli Geng|Guozhong Zhang