<img src="https://www.zjujournals.com/gcsjxb/images/1006-754X/images/logo.png" class="img-responsive">

图1 高速公路清雪铲结构示意图

1—路面；2—铲刃；3—铲体；4—积雪；5—清雪铲方通；6—加强筋；7—铲体焊接板；8—积雪挡板；9—焊接板；10—连接架方通；11—连接杆；12—右连接板；13—左连接板；14—中间连接板；15—三点悬挂装置安装底座。

Fig.1 Structural diagram of highway snow plow blade

根据笔者前期研究所建立的清雪铲除雪过程的理论力学模型和运动模型^[9]，可知：

1）当行进速度为10 m/s时，切削角在[34°, 37°]内增大、行进角在[50°, 53°]内增大，单位时间除雪体积持续提升至峰值0.928 5 m³/s。其中：行进角为铲体长度方向与路面宽度方向的夹角；切削角为铲刃与路面水平面的夹角。

2）当切削角为35.9°、行进角为51.9°时，行进速度在[25, 30] m/s内增大，单位时间除雪体积同步增大至峰值5.74 m³/s。

3）当采用35.9°切削角、51.9°行进角、28.7 m/s行进速度开展除雪作业时，清雪铲的最小抛雪高度可达1.5 m，最小抛雪位移（积雪沿路面宽度方向的移动距离）可达2.5 m。该参数组合为清雪铲的理论最优作业参数。

2 清雪铲除雪联合仿真原理

清雪铲的优化设计与性能验证需应对复杂的积雪条件，但传统物理试验成本高且周期长。为此，本文利用ANSYS/LS-DYNA软件对清雪铲与积雪的相互作用进行精准模拟，为清雪铲的作业参数设计提供科学依据；同时，使用Rocky DEM软件模拟积雪颗粒运动，并分析抛雪高度与抛雪位移。通过与理论分析结果进行对比，来验证联合仿真模型的准确性与可靠性，进而为清雪铲的结构设计和作业参数优化提供依据。

2.1　显式动力学仿真原理

清雪铲在除雪过程中与积雪层发生动态碰撞接触，使得积雪受到持续的冲击与挤压作用并逐步破碎、剥离，形成连续的力学响应及破坏演化过程。清雪铲将积雪剥离路面后，积雪沿着雪铲曲面运动并最终被抛至护栏外，该过程属于非线性结构动力学问题。其中，非线性既包括积雪的非线性，也包括清雪铲与积雪接触的非线性^[10]。

ANSYS/LS-DYNA软件是一款全球领先的显式动力学求解器，其核心优势体现在对高度非线性、瞬态冲击和多物理场耦合问题的高效精准求解能力。因此，本文使用ANSYS/LS-DYNA对清雪铲的除雪作业进行显式动力学仿真，以模拟清雪铲对积雪层的破碎剥离过程^[11]。

在有限单元法中，显式动力学的系统控制方程如下：

M \ddot{u} (t) + C \dot{u} (t) + K u (t) = F_{e x t} (t)

(1)

式中： M 为质量矩阵，用于描述系统的惯性特性； C 为阻尼矩阵，用于描述系统的阻尼效应； K 为刚度矩阵，用于描述系统的弹性特性； u $(t)$ 、 $\dot{u} (t)$ 、 $\ddot{u} (t)$ 分别为系统在时刻t的位移向量、速度向量和加速度向量； F_ext $(t)$ 为外力向量，表示作用于系统的外部载荷。

显式动力学法通过当前时间步的信息来计算下一个时间步的响应，而中心差分法通过离散化牛顿方程来得到显式的更新公式^[11]。设时间步长为 $Δ t$ ，则基于中心差分法的位移与速度的显式更新公式为：

\{\begin{array}{l} u (t + Δ t) = u (t) + Δ t \dot{u} (t) + \frac{1}{2} Δ t^{2} \ddot{u} (t) \\ \dot{u} (t + Δ t) = \dot{u} (t) + Δ t \ddot{u} (t) \end{array}

(2)

联立式(1)与式(2)，可直接求得每个时刻的加速度，具体表达式如下：

\ddot{u} (t) = M^{- 1} (F_{e x t} (t) - C \dot{u} (t) - K u (t))

(3)

显式动力学法的核心优势在于：每个时间步的位移求解仅依赖当前时刻的位移和速度状态，并可直接递推得到各节点的加速度响应，无需多次迭代求解，大幅降低了计算的复杂性与时间成本，计算效率高^[11]。该方法尤其适用于清雪铲与积雪层动态碰撞接触这类瞬态强非线性问题的分析。

2.2　离散单元法仿真原理

在清雪铲的除雪过程中，可将积雪视作大量积雪颗粒黏结而成的离散介质集合体。通过Rocky DEM 软件开展积雪颗粒的离散元仿真，以获取积雪的运动轨迹与抛离特性。

离散单元法将运动中的物体视作离散化的颗粒集合，通过模拟每个颗粒的运动状态及其相互作用来描述整体的宏观行为。该方法的核心在于牛顿第二定律与接触模型的耦合求解^[12]，具体可表示为：

\{\begin{array}{l} m_{i} \frac{d v_{i}}{d t} = \sum_{j = 1}^{N_{c}} F_{i j}^{c} + F_{i}^{g} + F_{i}^{a e r o} \\ I_{i} \frac{d ω_{i}}{d t} = \sum_{j = 1}^{N_{c}} (r_{i j} \times F_{i j}^{t} + M_{i j}^{r}) \end{array}

(4)

式中：m_i 为颗粒i的质量；I_i 为颗粒i的转动惯量； v_i 为颗粒i的平动速度； ω_i 为颗粒i的角速度； $F_{i j}^{c}$ 为颗粒间的接触力，包括法向分量和切向分量； $F_{i}^{g}$ 为颗粒i所受的重力， $F_{i}^{g} = m_{i} g$ ； $F_{i}^{a e r o}$ 为风雪耦合作用下的气动力； r_ij 为颗粒i质心指向接触点的位置矢量； $F_{i j}^{t}$ 为颗粒间的切向接触力； $M_{i j}^{r}$ 为颗粒间的滚动摩擦力矩；N_c为与颗粒i接触的颗粒数量。

Rocky DEM软件采用显式时间积分算法来迭代更新每个颗粒的运动状态：

\{\begin{array}{l} v_{i, t + Δ t} = v_{i, t} + Δ t \frac{F_{i, t}}{m_{i}} \\ x_{i, t + Δ t} = x_{i, t} + Δ t v_{i, t + Δ t} \end{array}

(5)

式中： v_i_,_t 为颗粒i在t时刻的速度； x_i_,_t 为颗粒i在t时刻的位移； F_i_,_t 为t时刻作用于颗粒i的合力（包括积雪颗粒间接触力、积雪颗粒-清雪铲接触力、风雪耦合气动力等），其大小由t时刻的速度与位移决定。

在Rocky DEM软件中建模时，积雪颗粒通常被视为可变形黏弹性球体，其法向接触力的大小可采用Hertz接触模型计算^[12]。在Hertz接触模型中，两接触物体间的法向接触力可表示为：

F_{n}^{e} = \frac{4}{3} E^{*} \sqrt[]{R^{*}} δ_{n}^{3 / 2}

(6)

其中：

E^{*} = {(\frac{1 - ν_{1}^{2}}{E_{1}} + \frac{1 - ν_{2}^{2}}{E_{2}})}^{- 1}

R^{*} = {(\frac{1}{R_{1}} + \frac{1}{R_{2}})}^{- 1}

式中： $F_{n}^{e}$ 为两接触物体间的法向接触力，E^*为等效弹性模量，R^*为等效接触半径，δ_n为法向重叠量（即接触侵入位移），E₁、E₂分别为2个接触物体的弹性模量，R₁、R₂分别为2个接触物体的半径， $ν$ ₁、 $ν$ ₂分别为2个接触物体的泊松比。

在Hertz接触模型中，耗散力包括法向阻尼力和切向阻尼力，其中法向阻尼力的表达式如下：

F_{n}^{d} = 2 \sqrt[]{\frac{5}{6}} β \sqrt[]{k_{n} m^{*}} v_{n}^{r e l}

(7)

其中：

β = \frac{- l n e}{\sqrt[]{l n^{2} e + π^{2}}}

k_{n} = 2 E^{*} \sqrt[]{R^{*} δ_{n}}

式中： $F_{n}^{d}$ 为法向阻尼力，m^∗为两接触物体的等效质量，v_n^rel为两接触物体的法向相对速度大小，β为法向阻尼系数，e为恢复系数，k_n为法向接触刚度。

在此基础上，两接触物体间的切向摩擦力（由法向接触力和法向阻尼力决定）可表示为：

F_{t} = m i n (μ (F_{n}^{e} + F_{n}^{d}), k_{t} δ_{t} + γ_{t} v_{t}^{r e l})

(8)

其中：

k_{t} = 8 G^{*} \sqrt[]{R^{*} δ_{n}}

γ_{t} = 2 \sqrt[]{\frac{5}{6}} β \sqrt[]{k_{t} m^{*}}

式中：F_t为两接触物体间的切向摩擦力，μ为摩擦系数，k_t为切向接触刚度，δ_t为切向相对位移，γ_t为切向阻尼系数，v_t^rel为切向相对速度大小，G^∗为等效剪切模量。

3 清雪铲除雪过程显式动力学仿真与参数优化

利用ANSYS/LS-DYNA软件对清雪铲的除雪作业进行显式动力学仿真分析，具体步骤为：构建并简化除雪三维模型、设置材料属性、划分网格、设置边界条件以及定义远程位移并求解。求解清雪铲与积雪层之间相互作用问题的关键在于准确构建清雪铲与积雪层间的接触模型，并选取合适的接触算法与参数^[13-14]。

3.1　清雪铲除雪仿真模型构建

清雪铲选用强度大、刚度高和耐磨损的结构钢材料，路面采用符合实际工况的沥青混凝土材料，两者的材料参数均满足除雪作业的要求，分别如表1和表2所示。

表1 结构钢材料参数

Table 1 Material parameters of structural steel

参数	数值
密度/(kg/m³)	7 850
杨氏模量/GPa	200
泊松比	0.3
体积模量/GPa	166.67
剪切模量/GPa	76.923

表2 沥青混凝土材料参数

Table 2 Material parameters of asphalt concrete

参数	数值
密度/(kg/m³)	2 400
杨氏模量/GPa	1.5
泊松比	0.35
体积模量/GPa	1.67
剪切模量/GPa	0.556

清雪铲所清除的积雪为新降浮雪，故选择＊MAT_PIECEWISE_LINEAR_PLASTICITY模型，所采用的状态方程为＊EOS_TABULATED_COMPACTION，用于描述积雪的体积压缩特性。积雪模型的材料参数如表3所示。

表3 积雪模型材料参数

Table 3 Material parameters of snow model

参数	数值
密度/(kg/m³)	300
杨氏模量/GPa	10
泊松比	0.3
抗压强度/MPa	0.1
剪切强度/MPa	0.03

清雪铲的结构参数如表4所示，利用SolidWorks软件建立其三维模型，如图2所示。

表4 清雪铲结构参数 (mm)

Table 4 Structural parameters of snow plow blade

参数	数值
铲体长度	3 500
直面板高度	1 500
铲板厚度	100
曲面板大端半径	500
曲面板小端半径	250

图2

图2 清雪铲三维模型

Fig.2 Three-dimensional model of snow plow blade

为便于后续的动力学仿真，对清雪铲三维模型进行适当简化，去除无关结构、倒角和圆角，其简化模型如图3所示。

图3

图3 清雪铲简化模型

Fig.3 Simplified model of snow plow blade

随后，采用四面体及六面体单元对清雪铲除雪模型进行网格划分，如图4所示。网格划分后，除雪模型的节点数量为8 620，单元数量为12 883；平均单元质量为0.765 73，处于0.75~0.81之间，验证了模型单元质量的可靠性。

图4

图4 清雪铲除雪模型网格划分

Fig.4 Grid division for snow removal model of snow plow blade

3.2　清雪铲除雪的动力学仿真模型构建

3.2.1　不同切削角下的动力学仿真模型

在清雪铲的除雪作业中，合理选择切削角不仅可避免因角度过小而导致积雪残留和除雪效率降低，还可避免因角度过大而导致路面和铲板损坏、除雪阻力增加以及除雪成本提高。为进一步确定清雪铲除雪时的最优切削角区间（固定行进角为0°），对切削角的理论最优区间[34°, 37°]进行等分，确定以34°~35°、35°~36°、36°~37°这3个区间作为切削角的3个水平，显式动力学仿真时分别取34.5°、35.5°和36.5°表征3个区间。

为了减少动力学仿真中运动碰撞出现振荡等不稳定现象，以及更准确地模拟动态的除雪过程，将时间步设为4步，分别为40、80、120、160 ms，对应的远程位移为400、800、1 200、1 600 mm；清雪铲与积雪间的接触属性设为侵蚀面-面接触，安全系数为0.9，沙漏系数为0.1，静摩擦系数为0.2，动摩擦系数为0.1。选择拉格朗日算法作为接触求解算法。对清雪铲施加远程位移约束，并将清雪铲和路面设为刚体。不同切削角下清雪铲与路面的约束设置如表5所示。表中：X、Y、Z分别表示沿X、Y、Z轴方向的位移；θ_X 、θ_Y 、θ_Z 分别表示绕X、Y、Z轴的旋转角度。

表5 不同切削角下清雪铲与路面的约束设置

Table 5 Constraint settings of snow plow blade and road surface under different cutting angles

约束	清雪铲	路面
X/mm	0	0
Y/mm	0	0
Z/mm	400、800、1 200、1 600	0
θ_X /(°)	0	0
θ_Y /(°)	0	0
θ_Z /(°)	0	0

根据上述约束条件，利用ANSYS/LS-DYNA软件建立不同切削角下的清雪铲除雪仿真模型，如图5所示。

图5

图5 不同切削角下的清雪铲除雪仿真模型

Fig.5 Snow removal simulation models of snow plow blade under different cutting angles

3.2.2　不同行进角下的动力学仿真模型

根据文献[9]中的清雪铲除雪理论模型，可得最优切削角为35.9°。在此基础上，为进一步确定清雪铲除雪时的最优行进角区间，对行进角的理论最优区间[50°, 53°]进行等分，确定以50°~51°、51°~52°、52°~53°这3个区间作为行进角的3个水平，显式动力学仿真时分别取50.5°、51.5°和52.5°表征3个区间。

同样地，设置4个时间步，分别为40、80、120、160 ms，对应的远程位移为1 000、2 000、3 000、4 000 mm；清雪铲与积雪间的接触属性设为侵蚀面-面接触，安全系数为0.9，沙漏系数为0.1，静摩擦系数为0.2，动摩擦系数为0.1。不同行进角下清雪铲与路面的约束设置如表6所示。

表6 不同行进角下清雪铲与路面的约束设置

Table 6 Constraint settings of snow plow blade and road surface under different travel angles

约束	清雪铲	路面
X/mm	0	0
Y/mm	0	0
Z/mm	1 000、2 000、3 000、4 000	0
θ_X /(°)	0	0
θ_Y /(°)	0	0
θ_Z /(°)	0	0

根据上述约束条件，利用ANSYS/LS-DYNA软件建立切削角为35.9°时不同行进角下的清雪铲除雪仿真模型，如图6所示。

图6

图6 不同行进角下的清雪铲除雪仿真模型（切削角为35.9°）

Fig.6 Snow removal simulation models of snow plow blade under different travel angles (with cutting angle of 35.9°)

3.2.3　不同行进速度下的动力学仿真模型

根据文献[9]中的清雪铲除雪理论模型，可得切削角为35.9°时的最优行进角为51.9°。在此基础上，为进一步确定清雪铲除雪时的最优行进速度区间，对行进速度的理论最优区间[27, 30] m/s进行等分，确定以27~28、28~29、29~30 m/s这3个区间作为行进速度的3个水平，显式动力学仿真时分别取27.5、28.5、29.5 m/s表征3个区间。

同样地，将时间步设为4步，分别为40、80、120、160 ms，不同行进速度对应不同的远程位移；清雪铲与积雪间的接触属性设为侵蚀面-面接触，安全系数为0.9，沙漏系数为0.1，静摩擦系数为0.2，动摩擦系数为0.1。不同行进速度下清雪铲与路面的约束设置如表7所示。

表7 不同行进速度下清雪铲与路面的约束设置

Table 7 Constraint settings of snow plow blade and road surface under different travel speeds

约束	清雪铲			路面
约束	行进速度为27.5 m/s	行进速度为28.5 m/s	行进速度为29.5 m/s	路面
X/mm	0	0	0	0
Y/mm	0	0	0	0
Z/mm	1 100、2 200、3 300、4 400	1 140、2 280、3 420、4 560	1 180、2 360、3 540、4 720	0
θ_X /(°)	0	0	0	0
θ_Y /(°)	0	0	0	0
θ_Z /(°)	0	0	0	0

根据上述约束条件，利用ANSYS/LS-DYNA软件建立切削角为35.9°、行进角为51.9°时不同行进速度下的清雪铲除雪仿真模型，如图7所示（通过调整行进速度即可获得不同条件下的模型）。

图7

图7 不同行进速度下的清雪铲除雪仿真模型（切削角为35.9°，行进角为51.9°）

Fig.7 Snow removal simulation models of snow plow blade under different travel speeds (with cutting angle of 35.9° and travel angle of 51.9°)

3.3　除雪性能仿真结果分析

3.3.1　不同切削角下的除雪性能

等效应力可将复杂的多轴应力等效为单向拉伸应力，本文通过显式动力学仿真获得不同切削角下积雪的剪切应力曲线和等效应力云图，用于表征清雪铲的除雪性能。不同切削角下积雪的剪切应力曲线如图8所示。

图8

图8 不同切削角下积雪的剪切应力曲线

Fig.8 Shear stress curves of snow under different cutting angles

分析图8可得：与切削角为34°~35°时相比，35°~36°和36°~37°切削角下积雪的剪切应力在0—50 ms内快速增大，且均在t=48 ms时达到峰值，该峰值比切削角为34°~35°时的峰值早64 ms出现，说明清雪铲在35°~36°和36°~37°这2个切削角区间下除雪时，积雪表现出更快速的应力积累过程；36°~37°切削角下积雪的剪切应力峰值为46.893 MPa，比35°~36°切削角下的剪切应力峰值45.887 MPa大，说明在切削角为36°~37°的清雪铲除雪过程中，积雪所受的径向压力更大，清雪铲对积雪的破碎效应更显著。35°~36°和36°~37°切削角下积雪的剪切应力在50—160 ms均处于缓慢变化阶段，但35°~36°切削角下的平均剪切应力大于36°~37°切削角下（35.446 MPa> 33.767 MPa），表明切削角为35°~36°的清雪铲对积雪的破碎效果更显著且更稳定。因此，确定35°~36°为清雪铲切削角的最优区间，这样既能保证较高的除雪效率，又能有效降低能耗和减少设备磨损。

图9所示为不同切削角下清雪铲除雪过程中积雪的等效应力云图（限于篇幅，仅展示t=40, 80, 120, 160 ms四个关键时刻的分布云图）。

图9

图9 不同切削角下积雪的等效应力云图

Fig. 9 Equivalent stress cloud maps of snow under different cutting angles

结合图9可得：在35°~36°和36°~37°切削角下，积雪的等效应力在t=48 ms时达到最大，表明此时清雪铲的除雪性能优越，而切削角为34°~35°时积雪的等效应力仍处于快速增大阶段，其峰值出现的时间比其余2个切削角区间晚64 ms；在50—160 ms的稳定除雪阶段，切削角为35°~36°时积雪的平均等效应力比切削角为36°~37°时大，表明切削角为35°~36°时清雪铲对积雪的破碎效果更好。综上，切削角为35°~36°时清雪铲具有较高的除雪效率和较好的除雪稳定性，即35°~36°为切削角的最优区间。

3.3.2　不同行进角下的除雪性能分析

在理论最优切削角35.9°下，通过显式动力学仿真探究行进角对清雪铲除雪性能的影响。不同行进角下积雪的剪切应力曲线如图10所示。

图10

图10 不同行进角下积雪的剪切应力曲线（切削角为35.9°）

Fig.10 Shear stress curves of snow under different travel angles (with cutting angle of 35.9°)

分析图10可得：当行进角分别为50°~51°、51°~52°和52°~53°时，积雪的剪切应力在0—120 ms内总体上呈增大趋势；51°~52°和52°~53°行进角下积雪的剪切应力在t=112 ms时达到峰值，比50°~51°行进角下的峰值早出现32 ms，说明清雪铲在51°~52°和52°~53°这2个行进角区间下除雪时，积雪表现出更快速的应力积累过程；51°~52°行进角下积雪的剪切应力峰值15.935 MPa比52°~53°行进角下的剪切应力峰值14.488 MPa大，这说明在行进角为51°~52°的清雪铲除雪过程中，积雪所受的径向压力更大，清雪铲对积雪的破碎效应更显著。在120—160 ms的除雪阶段，51°~52°和52°~53°行进角下积雪的剪切应力缓慢下降，但行进角为51°~52°时的平均剪切应力比行进角为52°~53°时的大（15.400 MPa>13.351 MPa），表明行进角为51°~52°的清雪铲对积雪的破碎效果更显著且更稳定。因此，确定51°~52°为清雪铲行进角的最优区间。

图11所示为不同行进角下清雪铲除雪过程中积雪的等效应力云图（限于篇幅，仅展示t=40, 80, 120, 160 ms四个关键时刻的分布云图）。

图11

图11 不同行进角下积雪的等效应力云图（切削角为35.9°）

Fig. 11 Equivalent stress cloud maps of snow under different travel angles (with cutting angle of 35.9°)

结合图11可得：在51°~52°和52°~53°行进角下，积雪的等效应力在t=112 ms时达到最大，表明此时清雪铲的除雪性能优越，而行进角为50°~51°时积雪的等效应力仍处于快速增大阶段，其峰值出现的时间比其余2个行进角区间晚32 ms；在120—160 ms的稳定除雪阶段，行进角为51°~52°时积雪的平均等效应力比行进角为52°~53°时大，表明行进角为51°~52°的清雪铲对积雪的破碎效果更好。综上，行进角为51°~52°时清雪铲具有较高的除雪效率和较好的除雪稳定性，即51°~52°为行进角的最优区间。

3.3.3　不同行进速度下的除雪性能分析

在理论最优切削角35.9°、理论最优行进角51.9°的条件下，通过显式动力学仿真探究行进速度对清雪铲除雪性能的影响。不同行进速度下积雪的剪切应力曲线如图12所示。

图12

图12 不同行进速度下积雪的剪切应力曲线（切削角为35.9°，行进角为51.9°）

Fig.12 Shear stress curves of snow under different travel speeds (with cutting angle of 35.9° and travel angle of 51.9°)

分析图12可得：当行进速度分别为27~28、28~29、29~30 m/s时，积雪的剪切应力在0—120 ms内总体上呈增大趋势；28~29 m/s行进速度下积雪的剪切应力在t=112 ms时率先达到峰值，该峰值比27~28 m/s和29~30 m/s行进速度下的峰值分别早出现48 ms和16 ms，说明清雪铲在行进速度为28~29 m/s下除雪时，积雪表现出更快速的应力积累过程；同时，28~29 m/s行进速度下积雪的剪切应力峰值大于其余2个行进速度区间，说明在行进速度为28~29 m/s的清雪铲除雪过程中，积雪所受的径向压力更大，清雪铲对积雪的破碎效应更显著。在120—160 ms的除雪阶段，28~29 m/s行进速度下积雪的剪切应力缓慢减小，其平均剪切应力12.358 MPa大于27~28 m/s行进速度下的11.690 MPa和29~30 m/s行进速度下的11.369 MPa，表明行进速度为28~29 m/s的清雪铲对积雪的破碎效果更显著且更稳定。因此，确定28~29 m/s为清雪铲行进速度的最优区间。

图13所示为不同行进速度下清雪铲除雪过程中积雪的等效应力云图（限于篇幅，仅展示t=40, 80, 120, 160 ms四个关键时刻的分布云图）。

图13

图13 不同行进速度下积雪的等效应力云图（切削角为35.9°，行进角为51.9°）

Fig.13 Equivalent stress cloud maps of snow under different travel speeds (with cutting angle of 35.9° and travel angle of 51.9°)

结合图13可得：当行进速度为28~29 m/s时，积雪的等效应力在t=112 ms时率先达到最大，表明此时清雪铲的除雪性能优越，而行进速度为27~28 m/s时积雪的等效应力仍处于快速增大阶段，其峰值出现的时间比28~29 m/s行进速度区间晚48 ms；在120—160 ms的稳定除雪阶段，行进速度为28~29 m/s时积雪的平均等效应力比29~30 m/s时的大，表明该行进速度区间下清雪铲对积雪的破碎效果更好。综上，行进速度为28~29 m/s的清雪铲具有较高的除雪效率和较好的除雪稳定性，即28~29 m/s为行进速度的最优区间。

4 清雪铲抛雪轨迹离散元仿真分析

利用Rocky DEM软件对清雪铲的除雪过程进行离散元仿真，具体步骤如下：构建简化的除雪三维模型，设置材料属性，设置颗粒与壁面、颗粒间的相互作用参数，约束模型自由度，设置清雪铲的速度与加速度并求解。求解抛雪高度和抛雪位移的关键是构建准确的清雪铲-积雪颗粒接触模型^[15-18]。

4.1　清雪铲除雪的离散元仿真模型构建

基于清雪铲的简化结构，利用SolidWorks软件建立用于离散元仿真的清雪铲除雪模型，如图14所示。

图14

图14 用于离散元仿真的清雪铲除雪模型

Fig.14 Snow removal model of snow plow blade for discrete element simulation

将图14所示的清雪铲除雪模型导入Rocky DEM软件，设动量柔软因子为0.1，Y轴方向的重力加速度为-9.8 m/s²（Y轴方向与重力方向相反）；积雪颗粒的堆积密度为300 kg/m³，杨氏模量为100 MPa，泊松比为0.3。积雪颗粒间的静摩擦系数为0.3，动摩擦系数为0.1，积雪颗粒与壁面之间的静摩擦系数为0.3，动摩擦系数为0.14，恢复系数为0.3；设粒径为0~≤5 mm、5~≤8 mm、8~≤10 mm的积雪颗粒依次占30%、20%、50%，颗粒形状设为球形。积雪颗粒的填充效果如图15所示。其中，积雪填充区域的长、宽、高分别为200、1 200、700 mm，填充质量为500 kg。离散元仿真总时长设为3 s，时间步长为0.01 s，选取t=0.6、1.2、1.8、2.4 s时刻的积雪颗粒分布云图来表征清雪铲除雪性能的离散元仿真效果。选择GPU（graphics processing unit，图形处理单元）为计算模块，以提高仿真效率。

图15

图15 积雪颗粒填充效果

Fig.15 Filling effect of snow particles

在清雪铲除雪的离散元仿真模型中，不同行进速度下清雪铲的约束设置如表8所示。

表8 不同行进速度下清雪铲的约束设置

Table 8 Constraint settings of snow plow blade under different travel speeds

约束	行进速度/(m/s)
约束	28.0	28.5	29.0
X/mm	1 120、2 240、3 360、4 480	1 140、2 280、3 420、4 560	1 160、2 320、3 480、4 640
Y/mm	0	0	0
Z/mm	0	0	0
θ_X /(°)	0	0	0
θ_Y /(°)	0	0	0
θ_Z /(°)	0	0	0

4.2　抛雪轨迹仿真结果分析

由显式动力学仿真结果可知，在理论最优切削角35.9°及理论最优行进角51.9°下，清雪铲的最优行进速度区间为[28, 29] m/s。为节省仿真耗时，仅验证行进速度最优区间对清雪铲抛雪高度和抛雪位移的影响。通过离散元仿真得到该工况下的抛雪高度与抛雪位移（行进速度分别取28、28.5、29 m/s），结果如图16所示。

图16

图16 不同行进速度下的抛雪高度和抛雪位移

Fig.16 Snow throwing height and snow throwing displacement under different travel speeds

分析图16可知，当行进速度为28.5 m/s时，积雪率先到达清雪铲末端，对应时刻为t=0.547 s；当行进速度分别为28、29 m/s时，积雪分别在t=0.587 s和t=0.566 s时到达清雪铲末端。上述对应时刻下积雪的位移均为2.160 m左右。当行进速度为28 m/s时，t=2.0 s和t=2.1 s时刻的抛雪高度分别为2.350 m和1.662 m，抛雪位移分别为5.697 m和5.933 m；当行进速度为28.5 m/s时，t=2.0 s和t=2.1 s时刻的抛雪高度分别为2.658 m和1.983 m，抛雪位移分别为6.234 m和6.506 m；当行进速度为29 m/s时，t=2.0 s和t=2.1 s时刻的抛雪高度分别为2.836 m和2.174 m，抛雪位移分别为5.974 m和6.228 m。目前，我国高速公路标准车道的宽度为3.750 m，除雪要求是将积雪抛离至距离路面5.75 m远、1 m高的护栏外。综上，行进速度为28~29 m/s时清雪铲的抛雪高度与抛雪位移均满足高速公路的除雪要求，验证了28~29 m/s为行进速度最优区间的可靠性。

图17所示为切削角为35.9°、行进角为51.9°时清雪铲除雪过程中的积雪颗粒分布云图。

图17

图17 不同行进速度下的积雪颗粒分布云图（切削角为35.9°，行进角为51.9°）

Fig. 17 Cloud maps of snow particle distribution under different travel speeds (with cutting angle of 35.9° and travel angle of 51.9°)

结合图17可得：在t=0—0.6 s阶段，行进速度为28.5 m/s时积雪颗粒率先到达清雪铲曲面板末端，说明此时清雪铲的抛雪性能最优；在t=0.6—2.4 s的稳定抛雪阶段，行进速度为28.5 m/s时清雪铲均可将积雪抛越至距离路面5.75 m远、1 m高的护栏外；在t=2.4 s时，行进速度为28.5 m/s时清雪铲铲面上的残雪最少，说明该行进速度下清雪铲的清雪效率最高。综上，行进速度为28~29 m/s时清雪铲满足高速公路除雪作业对抛雪高度、抛雪位移的要求，且除雪效率高，验证了该行进速度区间为最优区间的可靠性。

4.3　参数关系拟合与分析

使用MATLAB软件对离散元仿真得到的清雪铲的抛雪高度和抛雪位移数据进行曲线拟合，并与除雪理论模型^[9]的计算结果进行对比，结果如图18所示。

图18

图18 抛雪高度与抛雪位移的拟合曲线

Fig.18 Fitting curves of snow throwing height and snow throwing displacement

分析图18可得：抛雪高度曲线的残差标准差为0.034 m；抛雪高度随着行进速度的增大呈先增大后减小的趋势，当行进速度达到28.5 m/s左右时，抛雪高度达到峰值。由此可知，抛雪高度满足除雪要求的最优行进速度区间为28~29 m/s。抛雪位移曲线的残差标准差为0.174 m，抛雪位移随着行进速度的增大单调递增。由此可知，抛雪位移满足除雪要求的最优行进速度区间为26~30 m/s。综上，可同时使抛雪高度和抛雪位移满足高速公路除雪要求的最优行进速度区间为28~29 m/s，该工况下的最小抛雪高度为1.360 m，最小抛雪位移为11.700 m。

5 基于响应面法的清雪铲作业参数优化与试验研究

5.1　试验设计与响应面模型构建

响应面法^[19]是统计学与应用数学深度融合的试验设计与优化方法，其核心思想是将系统的输出响应拟合为输入变量的函数，并通过响应曲面与等高线图直观表征输出响应与输入变量间的非线性关系。该分析方法可高效解决多因素耦合下的非线性数据处理、工艺参数寻优等问题。

本文以清雪铲的切削角A、行进角B和行进速度C为设计变量，以抛雪高度S₁和抛雪位移S₂为响应指标，开展响应面分析。采用中心复合设计构建17组试验方案，并采用图19所示的除雪试验装置开展相关测试。响应面试验设计方案及结果如表9所示。

图19

图19 清雪铲除雪试验装置

Fig.19 Snow removal test device for snow plow blade

表9 响应面试验设计方案与结果

Table 9 Response surface experimental design schemes and results

序号	切削角A/(°)	行进角B/(°)	行进速度C/(m/s)	抛雪高度S₁/m	抛雪位移S₂/m
1	35	52	29.0	1.415	12.225
2	37	51	28.5	1.415	12.225
3	35	52	28.0	1.415	12.225
4	35	51	28.5	1.415	12.225
5	36	52	28.5	1.450	12.600
6	36	52	28.5	1.450	12.600
7	36	52	28.5	1.450	12.600
8	36	51	28.0	1.450	12.600
9	36	52	28.5	1.450	12.600
10	36	52	28.5	1.450	12.600
11	36	53	28.0	1.450	12.600
12	37	52	29.0	1.415	12.225
13	37	53	28.5	1.415	12.225
14	36	53	29.0	1.450	12.600
15	36	51	29.0	1.415	12.225
16	37	52	28.0	1.380	11.850
17	35	53	28.5	1.380	11.850

采用Design-Expert软件对试验数据进行多元回归分析，构建以A、B、C为变量的二阶多项式模型。二阶多项式模型的表达式如下^[20]：

y = β_{0} + \sum β_{a} x_{a} + \sum β_{a a} x_{a}^{2} + \sum \sum β_{a b} x_{a} x_{b} + ε

(9)

式中：y为响应量，x_a 、x_b 为设计变量，β₀为常数项，β_a 为线性系数，β_aa 为平方项系数，β_ab 为交互项系数，ε为误差项。

通过方差分析对所构建的抛雪高度及抛雪位移的响应面模型进行显著性检验。结果表明：所构建的响应面模型具有较高的拟合精度，对应的P值均小于0.01，说明模型整体显著，可准确反映响应量与影响因素之间的非线性关系。

5.2　因素交互作用与最优作业参数分析

利用Design-Expert软件对除雪试验结果进行处理与拟合^[21]，并通过响应曲面直观分析抛雪高度S₁和抛雪位移S₂与切削角A、行进角B及行进速度C之间的关系，并探究因素两两交互作用的影响，进而确定最优的作业参数。清雪铲的作业参数对抛雪性能的影响如图20所示。

图20

图20 清雪铲作业参数对抛雪性能的影响

Fig.20 Influence of operation parameters on snow throwing performance for snow plow blade

分析图20可得：当行进速度固定为28.5 m/s时，随着行进角从51°增大至53°，切削角增大使得抛雪高度与抛雪位移均发生显著变化，说明行进角与切削角之间存在较强的交互效应。当行进角固定为52°时，随着行进速度从28 m/s增大至29 m/s，切削角增大同样引起抛雪高度与抛雪位移显著波动，说明切削角与行进速度的交互作用也十分显著。当切削角固定为36°时，行进速度与行进角的变化对抛雪高度和抛雪位移的影响均较小，说明行进角与行进速度的交互作用相对不显著。综合上述结果，可确定清雪铲的最优作业参数组合：切削角为36°、行进角为52°、行进速度为28.5 m/s。

6 结束语

本文基于显式动力学与离散单元法，对清雪铲的除雪过程进行了仿真分析；同时，通过响应面法对清雪铲的除雪试验数据进行了拟合，揭示了因素两两交互作用的影响规律，并确定了最优作业参数组合，从而验证了理论分析结果与仿真结果的可靠性，具体结论如下。

1）清雪铲切削角的最优区间为35°~36°。动力学仿真结果显示，切削角为35°~36°时积雪的剪切应力峰值较切削角为36°~37°时提前了64 ms出现，且对应的平均剪切应力明显更高，说明该切削角区间下清雪铲能够实现更高效的积雪切削。

2）清雪铲行进角的最优区间为51°~52°。动力学仿真结果显示，行进角为51°~52°时积雪的剪切应力峰值较行进角为50°~51°时提前了32 ms出现，且对应的平均剪切应力明显更高，表明该行进角区间更有利于积雪脱离铲面并提升清雪铲的抛雪效率。

3）清雪铲行进速度的最优区间为28~29 m/s。动力学仿真结果显示，行进速度为28~29 m/s时积雪的剪切应力峰值较27~28 m/s和29~30 m/s行进速度下分别提前了48 ms和16 ms出现；同时，28~29 m/s行进速度下积雪的平均剪切应力明显更高，表明该行进速度区间更有利于积雪的切削与抛射。

4）利用Rocky DEM软件对切削角为35.9°、行进角为51.9°、行进速度为28~29 m/s工况下的清雪铲除雪过程进行了离散元仿真验证。MATLAB拟合结果显示，在该工况下，清雪铲的最小抛雪高度为1.360 m，最小抛雪位移为11.700 m，可满足将积雪抛至距离路面5.75 m、1 m高护栏外侧的除雪要求。

5）响应面分析结果显示，切削角与行进角、行进速度之间存在显著的交互效应，而行进角与行进速度之间的交互作用相对较弱。为实现最优的除雪性能，通过响应面优化确定了清雪铲除雪作业时的最佳工作参数组合：切削角为36°、行进角为52°、行进速度为28.5 m/s。在该参数条件下，对应的抛雪高度为1.450 m，抛雪位移为12.600 m，满足高速公路除雪作业要求，并验证了理论分析结果与仿真结果的可靠性与工程适用性。

本文链接：https：//www.zjujournals.com/gcsjxb/CN/10.3785/j.issn.1006-754X.2026.05.196

参考文献

原文顺序

文献年度倒序

文中引用次数倒序

被引期刊影响因子

[1]

魏晖昕

晋中高速公路除雪车推雪铲宽度优化

［J］. 公路交通科技， 2022， 39（4）： 166-172.

WEI

H X

Optimizing width of snow shovel for Jinzhong expressway

［J］. Journal of Highway and Transportation Research and Development， 2022， 39（4）： 166-172.

[2]

WAN

X D

， WANG

， HU

， et al.

Research on snow removal method of ultra high voltage DC composite bushing based on hot water jet

［C］//20th International Conference on AC and DC Power Transmission 2024. Shanghai， Jul. 12-15， 2024.

[3]

ZHAO

H R

， GU

Z R

， ZHENG

， et al.

Self-prior guided diffusion model for marine snow removal

［C］//OCEANS 2024-Singapore. Singapore， Apr. 15-18， 2024.

[4]

SHEN

， XU

J G

， SHAO

X T

， et al.

Toward snow removal via the diversity and complexity of snow image

［C］//2022 IEEE International Conference on Image Processing. Bordeaux， Oct. 16-19， 2022.

[5]

ZHANG

Q J

， LI

S L

， WANG

， et al.

Research on cooperative formation operation of unmanned snow removal vehicles on airport runway

［C］//2024 5th International Conference on Artificial Intelligence and Electromechanical Automation. Shenzhen， Jun. 14-16， 2024.

[6]

NGUYEN

P H D

， TRAN

Constraint-based snowplow optimization model for winter maintenance operations

［J］. Transportation Research Part A： Policy and Practice， 2024， 179： 103911.

[7]

HAN

， DONG

Research on a kind of self-propelled sawing equipment for snow and ice removal

［C］//2022 International Conference on Cloud Computing， Big Data and Internet of Things. Wuhan， Oct. 22-23， 2022.

[8]

MOHAMMADI

， MOSCHOPOULOS

Performance investigation of roof-mounted PV systems in snowy conditions considering snow coverage and snow removal process

［C］//2020 IEEE Power & Energy Society General Meeting. Montreal， QC， Aug. 2-6， 2020.

[9]

孔诚琳，刘向东，吴洪山.

基于MATLAB清雪铲高速公路除雪的运动分析

［J/OL］. 机电工程技术， 2025： 1-11（2025-07-15）［2025-09-28］. .

URL [本文引用: 4]

KONG

C L

， LIU

X D

， WU

H S

Kinematic analysis of highway snow removal using snow plow blades based on MATLAB

［J］. Mechanical & Electrical Engineering Technology， 2025： 1-11（2025-07-15）［2025-09-28］. .

URL [本文引用: 4]

[10]

郭扬程，邢蕾，赵明，等.

基于LS-DYNA的弹性结构清雪铲缓冲性能仿真及试验

［J/OL］. 机电工程技术， 2025： 1-8（2025-07-25）［2025-09-28］. .

GUO

Y C

， XING

， ZHAO

， et al.

Simulation and experiment of buffer performance of elastic structure snow removal blade based on LS-DYNA

［J］. Mechanical & Electrical Engineering Technology， 2025： 1-8（2025-07-25）［2025-09-28］. .

[11]

赵明，刘向东，郭扬程，等.

基于ANSYS/LS-DYNA仿真的刷毛除雪性能分析

［J/OL］. 机电工程技术， 2025： 1-12（2025-07-11）［2025-09-28］. .

URL [本文引用: 3]

ZHAO

， LIU

X D

， GUO

Y C

， et al.

Analysis of brush snow removal performance based on ANSYS/LS-DYNA simulation

［J］. Mechanical & Electrical Engineering Technology， 2025： 1-12（2025-07-11）［2025-09-28］. .

URL [本文引用: 3]

[12]

叶涛，杨瑞，蒋钦，等.

基于Rocky DEM的高压辊磨机粉碎效果影响因素分析

［J］. 有色金属（选矿部分）， 2025（3）： 91-98.

[本文引用: 2]

， YANG

， JIANG

， et al.

Analysis of influencing factors of crushing effect of high pressure roller mill based on Rocky DEM

［J］. Nonferrous Metals Mieral Processing Section， 2025（3）： 91-98.

[本文引用: 2]

[13]

何晓宁，李梦竹，常学良，等.

基于LS-DYNA的圆盘式花生切根装置设计与试验

［J］. 农业工程学报， 2025， 41（21）： 32-42.

X N

， LI

M Z

， CHANG

X L

， et al.

Design and experiment of the LS-DYNA-based disk-type peanut root-cutting device

［J］. Transactions of the Chinese Society of Agricultural Engineering， 2025， 41（21）： 32-42.

[14]

颜洪斌，崔峻，袁培银.

基于LS-DYNA的船舶-浮箱碰撞动力特性研究

［J］. 舰船科学技术， 2025， 47（24）： 57-61.

YAN

H B

， CUI

， YUAN

P Y

Research on dynamic characteristics of ship pontoon collision based on LS-DYNA

［J］. Ship Science and Technology， 2025， 47（24）： 57-61.

[15]

李炎，张磊.

基于Rocky DEM的辊式制粉机磨辊结构优化

［J］. 食品与机械， 2024， 40（10）： 93-100.

， ZHANG

Optimization of roll structure in rollermill based on Rocky DEM discrete element simulation

［J］. Food & Machinery， 2024， 40（10）： 93-100.

[16]

金浩，蒋明洋，徐超，等.

基于离散元法的柔性片烟建模及仿真参数标定

［J］. 沈阳农业大学学报， 2023， 54（5）： 563-571.

JIN

， JIANG

M Y

， XU

， et al.

Modeling and simulation parameters calibration of flexible tobacco strips based on discrete element

［J］. Journal of Shenyang Agricultural University， 2023， 54（5）： 563-571.

[17]

曹成茂，丁武阳，刘忠滨，等.

基于Rocky DEM的宁前胡排种器设计与试验

［J］. 农业机械学报， 2023， 54（8）： 53-64.

CAO

C M

， DING

W Y

， LIU

Z B

， et al.

Design and experiment of seed disperser for Ning-guo Radix peucedani based on Rocky DEM

［J］. Transactions of the Chinese Society for Agricultural Machinery， 2023， 54（8）： 53-64.

[18]

秦宽，步坤亭，沈周高，等.

连续复式茶叶理条机优化设计与试验

［J］. 农业机械学报， 2023， 54（3）： 382-393.

QIN

， BU

K T

， SHEN

Z G

， et al.

Optimization design and experiment of continuous compound tea leaf barber

［J］. Transactions of the Chinese Society for Agricultural Machinery， 2023， 54（3）： 382-393.

[19]

MITA

A F

， LABIBA

M J

， MALIHA

S R

， et al.

Optimizing properties of concrete containing stone dust and ceramic tile： a response surface methodology approach

［J］. Cleaner Engineering and Technology， 2025， 29： 101083.

[20]

刘文广，阳浩，孙杰，等.

基于响应面法的高压截止阀开启过程分析与结构优化研究

［J/OL］. 机电工程， 2025： 1-11（2025-09-15）［2025-09-28］. .

LIU

W G

， YANG

， SUN

， et al.

Research on opening process analysis and structural optimization of high-pressure globe valve based on response surface methodology

［J］. Journal of Mechanical & Electrical Engineering， 2025： 1-11（2025-09-15）［2025-09-28］. .

[21]

许丽，汪溪远，崔紫微，等.

非均相混合体系热洗含油污泥及洗脱废液处理研究

［J］. 环境工程技术学报， 2025， 15（6）： 2068-2076.

， WANG

X Y

， CUI

Z W

， et al.

Research on the treatment technology of oil-containing sludge and elution waste liquid from hot washing in a heterogeneous mixed system

［J］. Journal of Environmental Engineering Technology， 2025， 15（6）： 2068-2076.