基于理想弹性模态分析的滚石对树木冲击力计算方法

豆红强; 李心彧; 王浩; 郑彬彬; 郭朝旭

doi:10.12404/j.issn.1671-1815.2505939

科学技术与工程 ›› 2026, Vol. 26 ›› Issue (19) :8311-8319 DOI: 10.12404/j.issn.1671-1815.2505939

建筑科学

基于理想弹性模态分析的滚石对树木冲击力计算方法

豆红强 ¹^,²^,³ ,
李心彧 ¹ ,
王浩 ¹^,²^,³^,^* ,
郑彬彬 ¹ ,
郭朝旭 ²^,³

作者信息 +

Calculation Method for Rockfall Impact Force on Trees Based on Ideal Elastic Modal Analysis

DOU Hong-qiang ¹^,²^,³ ,
LI Xin-yu ¹ ,
WANG Hao ¹^,²^,³^,^* ,
ZHENG Bin-bin ¹ ,
GUO Chao-xu ²^,³

Author information +

文章历史 +

摘要

树木是天然的滚石灾害生态防护屏障,准确计算滚石对树木的冲击力是生态防护林设计的关键。然而,当前有关滚石冲击树木的相关研究仍较为匮乏。为此,基于合理假设,首先构建了滚石冲击树木的简化动力学模型,开展了滚石冲击作用下的树木横向强迫振动模态分析,进一步基于理想弹性模型分析方法联合滚石冲击树木过程中的位移协调方程,推导了滚石对树木峰值冲击力的解析表达式。而后,依托具体的计算工况,将基于理想弹性模态分析方法所得的计算结果与经典Hertz碰撞接触理论、基于ANSYS/LS-DYNA数值仿真方法所得的峰值冲击力进行了对比分析。研究结果表明,所提方法假设基本合理,基于理想弹性模态分析方法所得的滚石对树木冲击力的变化趋势和大小均与基于树木采用弹性本构模型的数值仿真计算结果高度一致,可作为滚石冲击树木峰值冲击力的一种快速估算方法。研究成果可为滚石灾害生态防护工程的设计提供一定的科学依据。

Abstract

Trees are regarded as natural ecological barriers against rockfall hazards. Accurate calculation of the impact force of rockfall on trees is essential for the design of ecological protection forests. However, research on rockfall impact on trees is still insufficient. A simplified dynamic model for rockfall impacting trees was established based on reasonable assumptions. The lateral forced vibration of trees under rockfall impact was analyzed through modal analysis. An analytical expression for the peak impact force was derived by combining the displacement coordination equations during the impact process with the ideal elastic model analysis method. The peak impact forces calculated by the proposed method were compared with those obtained from the classical Hertz contact theory and the ANSYS/LS-DYNA numerical simulation under specific computational conditions. The results show that the assumptions of the proposed method are reasonable. The variation trend and magnitude of the impact force obtained by the ideal elastic modal analysis method are in good agreement with the numerical simulation results which adopt an elastic constitutive model for trees. Thus, the proposed method can serve as a rapid estimation tool for the peak impact force of rockfall on trees. Result provides a scientific basis for the design of ecological protection engineering against rockfall hazards.

Graphical abstract

关键词

滚石灾害 / 树木 / 冲击力 / 模态分析 / 计算方法

Key words

rockfall hazards / trees / impact force / modal analysis / calculation method

引用本文

引用格式 ▾

豆红强, 李心彧, 王浩, 郑彬彬, 郭朝旭. 基于理想弹性模态分析的滚石对树木冲击力计算方法[J]. 科学技术与工程, 2026, 26(19): 8311-8319 DOI:10.12404/j.issn.1671-1815.2505939

DOU Hong-qiang, LI Xin-yu, WANG Hao, ZHENG Bin-bin, GUO Chao-xu. Calculation Method for Rockfall Impact Force on Trees Based on Ideal Elastic Modal Analysis[J]. Science Technology and Engineering, 2026, 26(19): 8311-8319 DOI:10.12404/j.issn.1671-1815.2505939

登录浏览全文

4963

注册一个新账户忘记密码

滚石灾害是一种全球广泛分布的典型山区地质灾害,具有突发性、多发性与高度随机性等特点,严重威胁人民群众的生命财产安全,影响工农业生产与交通运输等关键领域^[1-3]。树木是天然的滚石灾害生态防护屏障。已有现场滚石释放试验表明,树木对滚石灾害具有显著的消能缓冲效应,对比其他典型滚石灾害的防护技术,树木防护能级具有工作范围宽、适应性强的特点,尤其对低冲击能量水平的滚石灾害防护具有良好的补充作用^[4-7]。其中,准确计算滚石对树木的冲击力是滚石生态防护设计的关键。

当前,国内外学者在滚石冲击力与运动参量计算方面的研究已取得了显著进展。在理论研究方面,何思明等^[8]基于Hertz接触理论研究了冲击荷载下土体的初始屈服压力与冲击速度;陈泰江等^[9]基于线性黏弹性接触理论构建了滚石冲击地面的物理模型,其计算得到的最大冲击力与作用时间与全尺试验结果吻合良好。在试验研究方面,Labiouse等^[10]通过滚石冲击实验,监测冲击过程中的加速度变化、冲击持续时间及缓冲层侵彻深度,提出了滚石冲击力的经验计算公式;Zhang等^[11]在不考虑滚石自旋影响的前提下,研究了不同入射角与砂垫层厚度对冲击过程的影响,依据冲量定理建立了法向冲击力及滚石侵彻深度的半经验-半理论公式;Pichler等^[12]模拟滚石冲击砂砾土垫层过程,提出了一种半经验法的冲击力计算模型;Kawahera等^[13]基于自由下落试验数据,提出了滚石冲击力的简化计算方法。在数值模拟方面,陈颖琪等^[14]基于Hertz理论与Thornton弹塑性假设,提出了滚石冲击力修正系数计算方法,解决了传统Hertz理论高估冲击力的问题。刘成清等^[15]基于滚石运动轨迹与接触力学理论,系统分析了坡度、坡面材料对滚石运动的影响,并利用ANSYS/LS-DYNA对被动柔性防护网进行了数值验证。汪峰等^[16]则通过非线性显式动力学方法,系统研究了输电塔在滚石撞击下的动力响应,提出了撞击力峰值拟合公式。但总体而言,现有研究多集中于滚石与地表或工程结构之间的碰撞与冲击力计算,对滚石与树木之间的碰撞及其冲击力计算则涉及较少,准确建立滚石冲击作用下树木的动力学响应模型仍面临较大挑战^[17-18]。

基于此,拟尝试基于合理假设构建滚石-树木的简化冲击动力学模型,开展滚石冲击作用下的树木横向强迫振动模态分析,基于理想弹性模型分析方法联合滚石冲击树木过程中的位移协调方程,推导滚石对树木峰值冲击力的解析表达式,并借助经典Hertz碰撞接触理论和数值仿真方法验证其合理性与准确性,以期为滚石灾害的生态防治和绿色减灾提供理论基础。

1 滚石冲击树木力学模型的基本假设与构建

树木几何形态复杂多样,其结构自上而下又包含冠层、树干、根-土系统等,因此,建立滚石冲击作用下具备完全形态结构的树木动力学模型较为困难。为简化分析,现拟忽略树木的冠层与根-土系统,并将树干视为等截面的完全弹性直杆,其与土体为刚性连接,此时,树干即为等截面的Euler-Bernoulli悬臂梁模型。

与此同时,滚石几何形态各异,但为便于分析,依据最不利原则,将滚石理想化为球形刚体^[19-20]。此外,滚石与树木的碰撞响应会由于碰撞物体和树木结构的力学性质、碰撞物体的质量和速度、碰撞物体冲击角度等不同而产生不同的影响^[5,21]。为便于开展研究且依据最不利原则,假设滚石正面冲击树木。又根据Zhang等^[11]研究,滚石实际旋转动能仅占总能量的约10%,因此可忽略滚石角速度的作用。即假设滚石以正面方向冲击树木,图1所示为简化后的滚石冲击作用下的树木力学模型。树干在承受滚石的冲击荷载过程中,将发生横向强迫振动,为便于后续开展分析,进一步做如下假设。

(1)树干各截面的中心惯性轴在同一平面XOY内。

(2)外荷载(冲击荷载)作用于XOY平面内。

(3)树干在XOY平面做横向振动。

(4)树干的主要变形为弯曲变形。

(5)树干在低频振动时可以忽略剪切变形和截面绕中性轴转动惯量的影响。

2 滚石对树木冲击力的计算方法

2.1 滚石冲击作用下的树木横向强迫振动模态分析

不失一般性,图2(a)为树木简化为等截面悬臂梁模型在任何荷载作用下的示意图,现取距边界x处长度为dx的微段进行受力分析,如图2(b)所示。

取竖向受力平衡,根据达朗贝尔原理,即可得到该悬臂梁在横向弯曲强迫振动下的波动方程为

(1)ρSdx$\frac{{\partial }^{2}\mathit{y}}{\partial {\mathit{t}}^{2}}$+$\left({\mathit{F}}_{\mathrm{s}}+\frac{\partial {\mathit{F}}_{\mathrm{s}}}{\partial \mathit{x}}\mathrm{d}\mathit{x}\right)$-F_s-f(x,t)dx=0

以右截面上任意一点为矩心,根据力矩平衡得

(2)$\left(\mathit{M}+\frac{\partial \mathit{M}}{\partial \mathit{x}}\mathrm{d}\mathit{x}\right)$-M-F_sdx-f(x,t)dx$\frac{\mathrm{d}\mathit{x}}{2}$+ρSdx$\frac{{\partial }^{2}\mathit{y}}{\partial {\mathit{t}}^{2}}\frac{\mathrm{d}\mathit{x}}{2}$=0

将式(2)予以化简,可得到等截面Euler-Bernoulli悬臂梁横向自由振动和强迫振动的动力学方程,即

(3)EI$\frac{{\partial }^{4}\mathit{y}(\mathit{x},\mathit{t})}{\partial {\mathit{x}}^{4}}$+ρS$\frac{{\partial }^{2}\mathit{y}(\mathit{x},\mathit{t})}{\partial {\mathit{t}}^{2}}$=0

(4)EI$\frac{{\partial }^{4}\mathit{y}(\mathit{x},\mathit{t})}{\partial {\mathit{x}}^{4}}$+ρS$\frac{{\partial }^{2}\mathit{y}(\mathit{x},\mathit{t})}{\partial {\mathit{t}}^{2}}$=f(x,t)

式中:EI为弯曲刚度。

结构的振动频率和模态由边界条件决定,悬臂梁在自由振动情况下,将以系统的固有频率进行振动。在此,依据振型叠加法,将悬臂梁在自由振动下的结构反应y(x,t)解耦为与时间相关的函数q(t)和与坐标相关的模态函数ϕ(x)的乘积,此时该悬臂梁的自由振动可表示为

(5)y(x,t)=$\stackrel{\mathit{\infty }}{\sum _{\mathit{n}=1}}$ϕ_n(x)q_n(t)

式(5)中:q_n(t)称为广义坐标下的第n阶振型,在自由振动下,与时间t相关的函数q_n(t)可以表示为q_n(t)=a_nsin(ω_nt+θ_n),a_n为振幅,θ_n为初相位,两者由系统的初始条件确定;ω_n为角频率。

将式(5)代入等截面悬臂梁自由振动动力学方程式(3)中,可得与时间和坐标相关的方程,即

(6)

$\left\{\begin{array}{l}\mathit{q}″\left(\mathit{t}\right)+{\mathit{\omega }}^{2}\mathit{q}\left(\mathit{t}\right)=0\\ \mathit{\varphi }″″\left(\mathit{x}\right)-{\mathit{a}}^{4}\mathit{\varphi }\left(\mathit{x}\right)=0\end{array}\right.$

(7)

$\left\{\begin{array}{l}\mathit{q}\left(\mathit{t}\right)=\mathit{A}\mathrm{c}\mathrm{o}\mathrm{s}\mathit{\omega }\mathit{t}+\mathit{B}\mathrm{s}\mathrm{i}\mathrm{n}\mathit{\omega }\mathit{t}\\ \mathit{\varphi }\left(\mathit{x}\right)={\mathit{A}}_{1}\mathrm{c}\mathrm{o}\mathrm{s}\mathit{a}\mathit{x}+{\mathit{A}}_{2}\mathrm{s}\mathrm{i}\mathrm{n}\mathit{a}\mathit{x}+{\mathit{A}}_{3}\mathrm{c}\mathrm{o}\mathrm{s}\mathrm{h}\mathit{a}\mathit{x}+{\mathit{A}}_{4}\mathrm{s}\mathrm{i}\mathrm{n}\mathrm{h}\mathit{a}\mathit{x}\end{array}\right.$

式中:a 为一个参数;A、B、A_i(i=1~4)为常数,由梁的边界条件所决定;ω为角频率。

由图2可知,等截面悬臂梁的边界条件可表示如下。

①固定端挠度和截面转角为0,即

(8)

$\left\{\begin{array}{l}{\mathit{y}}_{\left|\mathit{x}=0,\mathit{t}\right.}=0\\ \frac{\partial \mathit{y}}{\partial \mathit{x}}\left|{}_{\mathit{x}=0,\mathit{t}}\right.=0\end{array}\right.$

②自由端剪力和弯矩为0,即

(9)

$\left\{\begin{array}{l}\mathrm{M}=\mathrm{E}\mathrm{I}\frac{{\partial }^{2}\mathit{y}}{\partial {\mathit{x}}^{2}}\left|{}_{\mathit{x}=\mathit{l},\mathit{t}}\right.=0\\ {\mathit{F}}_{\mathrm{s}}=\frac{\partial \mathit{M}}{\partial \mathit{x}}=0\end{array}\right.$

经整理可得,等截面自由振动悬臂梁对应的各阶模态函数为

(10)$ \begin{aligned}f_{n}(x)= & C_{n}\left[\cos a_{n} x-\cosh a_{n} x-\frac{\sin a_{n} l-\sinh a_{n} l}{\sin a_{n} l+\sinh a_{n} l} \times\right. \\& \left.\left(\sin a_{n} x-\sinh a_{n} x\right)\right]\end{aligned}$

式(10)中:C_n为归一化系数,可通过振型的正交性条件确定;a_n为与悬臂梁振动阶数n相关的参数,由梁的边界条件和振动特性决定;l为悬臂梁的长度。

进一步的,根据模态函数ϕ_i(x)和ϕ_j(x)的正交性条件

(11)

${\int }_{0}^{\mathit{l}}$

f_n(x)f_m(x)dx=δ_nmρSl

式(11)中:δ_nm为克罗内克δ函数。

将式(4)强迫振动下悬臂梁响应模态予以展开,得

(12)y(x,t)=$\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}$q_n(t)f_n(x)

将模态展开代入强迫振动控制方程,并结合模态函数正交性,化简得到

(13)M_n$\stackrel{··}{{\mathit{q}}_{\mathit{n}}}$+M_n

${\mathit{\omega }}_{\mathit{n}}^{2}$

q_n=${\int }_{0}^{\mathit{l}}$F(x,t)f_n(x)dx

同时,定义广义力为

(14)F_n(t)=${\int }_{0}^{\mathit{l}}$(x,t)f_n(x,t)dx

则模态方程式(13)可表示为

(15)$\stackrel{··}{{\mathit{q}}_{\mathit{n}}}$+${\mathit{\omega }}_{\mathit{n}}^{2}$q_n=$\frac{{\mathit{F}}_{\mathit{n}}\left(\mathit{t}\right)}{{\mathit{M}}_{\mathit{n}}}$

使用杜哈梅积分求解式(15),可得时域函数为

(16)q_n(t)=$\frac{1}{{\mathit{\omega }}_{\mathit{n}}{\mathit{M}}_{\mathit{n}}}{\int }_{0}^{\mathit{t}}$F_n(τ)sin[ω_n(t-τ)]dτ

再代入式(12)得

(17)y(x,t)=$\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}$f_n(x)$\frac{1}{{\mathit{\omega }}_{\mathit{n}}{\mathit{M}}_{\mathit{n}}}$×${\int }_{0}^{\mathit{t}}$F_n(τ)sin[ω_n(t-τ)]dτ

同时,结合模态主质量表达式

(18)M_n=${\int }_{0}^{\mathit{l}}$ρSf${}_{\mathrm{n}}^{2}$(x)dx

整理即可得悬臂梁在强迫振动下的动力响应为

(19)y(x,t)=$\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}\frac{{\mathit{f}}_{\mathit{n}}\left(\mathit{x}\right){\mathit{f}}_{\mathit{n}}\left(\mathit{l}\right)}{{\mathit{\omega }}_{\mathit{n}}^{2}\mathit{\rho }\mathit{S}\mathit{l}}$×${\int }_{0}^{\mathit{t}}$F(τ)sin[ω_n(t-τ)]dτ

式(19)中:f_n(x)为该悬臂梁结构的各阶模态函数;ω_n为悬臂梁第n阶主振动对应的自然角频率。

2.2 建立位移协调方程

就滚石冲击假设为完全弹性的树木结构而言,忽略滚石体的自身振动,那么,滚石撞击位置处滚石与树木之间的侵彻位移可通过变形协调方程表示为

(20)δ(t)=s₁(t)-s₂(t)

式(20)中:s₁(t)和s₂(t)分别为滚石和树木碰撞位置处的位移随时间变化的函数。

以Hertz接触模型^[22]作为滚石与树木碰撞接触模型,则滚石与树木之间的侵彻位移可以表示为

(21)δ(t)=${\left[\frac{\mathit{F}\left(\mathit{t}\right)}{{\mathit{K}}_{\mathrm{h}}}\right]}^{2/3}$

式(21)中:F(t)为碰撞接触力;δ(t)为滚石与树木侵彻位移变化函数;K_h为Hertz碰撞接触刚度,可按下式计算确定。

(22)K_h=$\frac{4}{3}$E^*

$\sqrt{\mathit{R}}$

(23)$\frac{1}{{\mathit{E}}^{\mathrm{*}}}$=$\frac{1-{\mathit{\nu }}_{1}^{2}}{{\mathit{E}}_{1}}$+$\frac{1-{\mathit{\nu }}_{2}^{2}}{{\mathit{E}}_{2}}$

式中:R为滚石的半径;E^*为滚石与树木间的等效弹性模量;E₁和E₂为滚石和树木的弹性模量;ν₁和ν₂为滚石和树木的泊松比。

冲击力F(t)为碰撞过程中的冲击力,m为滚石的质量,那么在时间间隔dτ内滚石的速度减小量为[F(τ)/m]dτ,若碰撞初始时刻为时间t,那么碰撞的位移减小量可以表示为[F(τ)/m]dτ(t-τ),因此,在t时刻,滚石的位移变化量可以表示为

(24)

$\left\{\begin{array}{l}{\mathit{s}}_{1}\left(\mathit{t}\right)={\mathit{v}}_{0}\mathit{t}-\frac{1}{\mathit{m}}{\int }_{0}^{\mathit{t}}\mathit{F}\left(\mathit{\tau }\right)(\mathit{t}-\mathit{\tau })\mathrm{d}\mathit{\tau }\\ {\mathit{s}}_{1}\left(0\right)=0,\mathit{s}\mathit{\text{'}}{}_{1}\left(0\right)={\mathit{v}}_{0}\end{array}\right.$

式(24)中:v₀为滚石碰撞前的冲击速度。

如图3所示为滚石冲击树木下的位移变形示意图。此时结合简化的滚石冲击树木动力学模型,滚石与树木碰撞位置处的位移随时间变化函数可由模态分析计算得出。

联立以上方程,即可获得滚石冲击力变化的积分方程,即

(25)${\left[\frac{\mathit{F}\left(\mathit{t}\right)}{{\mathit{K}}_{\mathrm{H}}}\right]}^{2/3}$=v₀t-$\frac{1}{\mathit{m}}{\int }_{0}^{\mathit{t}}$F(τ)(t-τ)dτ-$\stackrel{\mathit{n}}{\sum _{\mathit{i}=1}}\frac{{\mathit{f}}_{\mathit{n}}\left(\mathit{x}\right){\mathit{f}}_{\mathit{n}}\left(\mathit{l}\right)}{{\mathit{\omega }}_{\mathit{n}}^{2}\mathit{\rho }\mathit{S}\mathit{l}}{\int }_{0}^{\mathit{t}}$F(τ)sin[ω_n(t-τ)]dτ

直接求解包含无穷积分的偏微分方程较为困难,为此可仅针对该简化悬臂梁结构的一阶模态响应进行计算,且由于Hertz接触模型是理想弹性状况下的计算模型,并未考虑滚石冲击过程中的能量损失,因此其滚石冲击力变化时程曲线应体现出对称性,曲线呈脉冲形状。因此,可以假设滚石冲击力变化时程曲线表达式为

(26)F(t)=P_maxsin(ωt)

式(26)中:P_max为滚石冲击力变化时程曲线的峰值冲击力;ω为曲线角频率,假设滚石冲击力在t_max时达到峰值,角速度ω可表示为π/(2t_max)。

而后将理想滚石峰值冲击力变化时程曲线代入位移协调方程中,并取t_max时刻为计算时间节点即可计算滚石碰撞中的冲击力变化。

3 模型验证与讨论

为验证本文所提模型的准确性与合理性,现选取一假想工况开展滚石-树木冲击力的计算对比验证与分析,即分别采用经典Hertz碰撞接触理论、本文所提的理想弹性模态分析方法和基于ANSYS/LS-DYNA数值仿真方法计算不同滚石冲击速度和滚石半径下的滚石冲击树木的峰值冲击力。有必要说明的是,在采用ANSYS/LS-DYNA开展滚石冲击树木的数值仿真时,为丰富并提高计算结果的可比性,针对树木的本构模型分别选取了线弹性模型和ANSYS/LS-DYNA内置的木材模型(横观各向同性材料)。其中,木材本构行为是通过弹性本构关系、失效准则、塑性流动、应变硬化、峰后软化以及应变率效应等进行描述^[23-24],由于涉及计算参数较多且非本文重点研究对象,此处仅给出相应的弹性参数,即木材弹性模量为11 350 MPa,泊松比为0.35,密度为673.1 kg/m³,其他计算参数取值详见文献[18]。

鉴于滚石灾害具有高度随机性,为研究不同滚石冲击速度和滚石半径等工况下的滚石对树木的冲击力,设计如表1所示的碰撞计算方案。此外,假定树木高度为5 m,截面呈等值圆形截面,其直径为0.25 m,滚石与树木发生正面碰撞,碰撞高度为1.5 m。

3.1 对比验证

如图4和图5所示分别为基于经典Hertz接触理论、本文所提的理想弹性模态分析方法和数值仿真计算的不同滚石冲击速度和滚石半径下的滚石冲击树木的峰值冲击力变化规律。由图4和图5可知,前述三种方法所得的滚石对树木的峰值冲击力均随滚石冲击速度和滚石半径的增大而增大,这也符合基本常识。

同时,由图4、图5可知,基于经典Hertz接触理论所得的峰值冲击力最大,且随着滚石速度和半径的增加愈加偏离其他方法的结果。这是由于Hertz接触理论基于弹性接触假设,其计算结果主要取决于滚石速度与接触刚度,且该方法主要适用于冲击体尺寸远小于被冲击体的情况,而此时滚石的直径已与树木的直径基本一致,导致其冲击效应被严重高估,表明经典Hertz接触理论已不适用于滚石与树木的碰撞计算。

值得说明的是,所提的理想弹性模态分析方法所得的峰值冲击力与基于ANSYS/LS-DYNA木材弹性本构模型所得结果基本一致。但是,这两种方法均未考虑树木在滚石冲击过程中的塑性变形和能量吸收效应,因此,其计算得到的峰值冲击力明显高于基于ANSYS/LS-DYNA木材本构模型所得的结果。此外,由于本文所提的理想弹性模态分析方法仅使用第一阶模态进行计算,由式(19)计算所得的树木位移小于基于树木弹性本构模型所得的结果,导致滚石与树木的碰撞接触时间也较短,进而致使由此得到的滚石峰值冲击力相对较大。但总的来说,本文所提的理想弹性模态分析方法可较好地反映滚石冲击树木过程中的峰值冲击力大小,且适用于不同滚石冲击速率和滚石直径。

3.2 讨论分析

3.2.1 基于理想弹性模态分析方法的冲击力时程曲线

如图6和图7所示分别为考虑不同冲击速率和不同滚石半径下基于理想弹性模态分析方法所得的滚石-树木冲击力的时程曲线。可知,由于不考虑碰撞过程中的能量耗散,基于理想弹性模态分析方法所得的滚石冲击力变化时程曲线均为对称曲线,滚石的冲击压缩阶段与弹性回弹阶段互为逆过程,当然这也是式(26)假设的结果。此外,由图亦可知,滚石冲击树木的作用时间并不随滚石冲击速度和滚石直径的变化而发生明显变化,进一步表明滚石冲击树木的动力学响应主要受滚石动能的影响,进而影响滚石峰值冲击力的大小。

3.2.2 基于数值仿真的冲击力时程曲线

由于数值仿真可以给出更为丰富多样的计算结果,该部分拟以半径为0.1 m的球形滚石,冲击速度设定为10 m/s的工况开展数值模拟,其中木材的力学特征分别采用弹性本构模型和木材本构模型予以模拟。以数值仿真计算结果为依托,深入讨论本文所提理想弹性模态分析方法的优劣性。

(1)弹性本构模型。如图8所示为基于ANSYS/LS-DYNA且树木采用弹性本构模型所得的滚石运动参量(冲击速度、位移与冲击力)的变化时程曲线。可知,滚石冲击树木的过程发生于极短时间内,约为1.15 ms(对应于树木-滚石碰撞后脱离,对应冲击力变为0的时刻),这与本文中基于理想弹性模态分析方法所得的冲击时间高度一致。由图8(a)进一步可知,该工况下的滚石冲击力时程曲线呈对称形态,这与本文理想弹性模态分析方法所假设的一致,表明本文假设是合理的。

同时,根据图8(a)和图8(b)可把滚石冲击树木的过程划分为3个阶段:冲击压缩阶段、回弹阶段和稳态阶段。①冲击压缩阶段:滚石与树木初始接触瞬间,滚石位移时程曲线与树木碰撞点位移曲线重合,表明二者相互侵彻。此阶段滚石冲击力迅速攀升至峰值(252 kN),同时其运动速度剧烈衰减,从初始的10 m/s降至约1 m/s,冲击方向上的侵彻位移达到最大值。②回弹阶段:由于树木采用弹性本构模型,滚石随树木发生弹性回弹,其速度恢复至4 m/s。伴随侵彻位移减小,冲击力迅速下降,其时程曲线呈三角函数状,并以峰值力为对称轴呈现对称变化,这也印证该模型下碰撞过程无能量耗散,且滚石冲击压缩阶段与弹性回弹阶段互为逆过程。③稳态阶段:滚石以4 m/s的速度沿反方向运动,碰撞过程结束。综上可知,在此模型中,树木对滚石的拦挡作用主要表现为改变其运动模式,致使滚石速度方向发生反转。由于弹性模型未考虑能量损失,滚石在冲击过程中未发生能量耗散,弹性树木对滚石冲击能量的吸收率约为60%。

(2)木材本构模型。同样的,如图9所示树木采用木材本构模型所得的滚石运动参量(冲击速度、位移与冲击力)的变化时程曲线。可知,当采用木材本构模型时,由于该模型反映了树木在动态加载下的“应变硬化-应变软化”特征,致使其冲击过程中会产生较大的塑性变形,从而起到缓冲作用,具体表现为冲击接触时间增大(约为2.3 ms),滚石冲击过程中伴随有明显的冲击能量损失。具体来讲,在冲击压缩阶段,滚石与树木初始接触瞬间,滚石位移时程曲线与树木碰撞点位移曲线重合,表明二者相互侵彻。此阶段滚石冲击力迅速攀升至峰值(65 kN),同时其运动速度剧烈衰减,从初始的10 m/s降至约2 m/s,表明整个过程中树木对滚石冲击的能量吸收率在80%。在回弹阶段,伴随侵彻位移减小,滚石冲击力迅速下降,其时程曲线呈现非对称形态,且冲击加载阶段的持续时间小于回弹阶段,进一步证实该木材本构模型下的碰撞过程存在显著能量耗散。

对比两种本构模型下的滚石运动参量变化可知,弹性本构模型因忽略碰撞过程中的塑性耗能,导致峰值冲击力(252 kN)被显著高估(木材模型峰值力为65 kN),且木材本构模型下的能量吸收率(80%)高于弹性模型(60%)。可见,从理论上讲,木材本构模型可考虑滚石-树木碰撞中的能量耗散,更能准确反映实际中的滚石冲击行为。然而,由于木材本构模型过于复杂,很难借此推导滚石冲击力的解析解。相反的,虽然基于弹性本构模型所得的计算结果难以如实反映滚石-树木的冲击过程,但其计算结果也相对合理。而本文所提的理想弹性模态分析方法与基于ANSYS/LS-DYNA的弹性本构模型所得结果高度一致,表明本文所提方法亦具有一定的合理性和适用性。

4 工程案例与其应用

某典型滚石灾害隐患点位于福建省泉州市,其斜坡宽200~460 m,坡度介于20°~40°,坡顶底相对高差L为130~150 m,坡面植被覆盖较差,主要为草灌局部稀疏乔木。该坡面分布大量孤、滚石及危岩体,部分滚石无埋深直接散落于坡面上,其余大部分孤石体呈半嵌入坡体状态。

基于现场勘查与设计资料选取直径约为0.3 m的近似球形滚石予以分析,根据斜坡地形特征估算,滚石抵达防护林带时的冲击速度取8 m/s。用于防护的木荷平均直径取0.22 m,平均高度5 m,将其简化为等截面悬臂梁,冲击作用点高度取1.5 m,仅考虑其一阶模态响应,结合位移协调方程式(20)与滚石冲击力变化时程曲线式(26),通过求解积分方程式(25),计算得到在此工况下,滚石对树木的峰值冲击力约为108 kN,由此产生的最大弯矩约为162 kN·m。

为确保坡脚人员和建筑安全,拟采用如图10所示的孤石分解、主/被动防护网、基础嵌补以及截排水等岩土工程措施,并辅以乔木种植形成梯级生态防护屏障进行崩塌滚石灾害防治。为判断滚石冲击力下乔木的防护可行性,初步选取直径为0.22 m的木荷,其树干在冲击点处横截面的抗弯截面模量约为1.05×10^-3m³,抗弯强度一般在90~110 MPa^[25]。取其下限保守估计,则单株木荷所能提供的临界弯矩约为95 kN·m,可见仅采用单株乔木难以抵抗滚石的峰值冲击力。因此,在设计时采取梅花型(株间混交)的方式,种植穴规格统一为0.5 m×0.5 m×0.5 m,株距保持在2.5 m×2.5 m。

5 结论

(1)构建了滚石冲击树木的简化动力学模型,基于理想弹性模态分析法建立了相应的位移协调方程,借此推导了滚石冲击树木峰值冲击力的解析解,并通过对比分析ANSYS/LS-DYNA所得的数值解以及经典Hertz接触理论所得结果,验证了本文方法的合理性与适用性,可作为滚石冲击树木最大冲击力的一种估算方法。

(2)在不同滚石冲击速度和滚石半径等碰撞工况下,基于本文所提的理想弹性模态分析方法和数值仿真计算得到的滚石峰值冲击力变化趋势一致,尤其与树木采用弹性本构模型所得的结果高度一致。而经典Hertz接触理论已不再适用于求解滚石冲击树木过程中的峰值冲击力。

(3)值得说明的是,在构建滚石冲击树木的简化动力学模型时作了诸多假设,如忽略根土系统的作用且将树木简化为等截面的完全弹性悬臂梁等,尤其是本文已论证树木采用木材本构模型更能反映实际情况,这些都有待进一步研究。

参考文献

原文顺序 | 出版日期 | 本文引用

[1]	Ye P, Yu B, Chen W H, et al. Fracture characterization based on data fusion technology and its application in rockfall hazard assessment[J]. Environmental Earth Sciences, 2024, 83(7): 208.

[2]	Graber A, Santi P. Power law models for rockfall frequency-magnitude distributions: review and identification of factors that influence the scaling exponent[J]. Geomorphology, 2022, 418(1): 135-157.

[3]	罗刚, 程谦恭, 沈位刚, 等. 高位高能岩崩研究现状与发展趋势[J]. 地球科学, 2022, 47(3): 913-934.

[4]	Luo Gang, Cheng Qiangong, Shen Weigang, et al. Research status and development trend of the high-altitude extremely-energetic rockfalls[J]. Earth Science, 2022, 47(3): 913-934.

[5]	黄润秋, 刘卫华, 龚满福, 等. 树木对滚石拦挡效应研究[J]. 岩石力学与工程学报, 2010, 29(S1): 2895-2901.

[6]	Huang Runqiu, Liu Weihua, Gong Manfu, et al. Study of trees resistance effect test on rolling rock blocks[J]. Chinese Journal of Rock Mechanics and Engineering, 2010, 29(S1): 2895-2901.

[7]	Caviezel A, Ringenbach A, Demmel S E, et al. The relevance of rock shape over mass implications for rockfall hazard assessments[J]. Nature Communications, 2021, 12(1): 5546.

[8]	Wu Y, Wang Z, Li X, et al. Blocking-baffle mechanism of vegetation to rock avalanche on the forested slope of the east Sichuan-Tibet Mountains[J]. Landslides, 2024: 1-16.

[9]	Dupire S, Bourrier F, Monnet J M, et al. The protective effect of forests against rockfalls across the French Alps: influence of forest diversity[J]. Forest Ecology & Management, 2016, 382: 269-279.

[10]	何思明, 李新坡, 吴永. 滚石冲击荷载作用下土体屈服特性研究[J]. 岩石力学与工程学报, 2008(S1): 2973-2977.

[11]	He Siming, Li Xinpo, Wu Yong. Research on yield property of soil under rock-fall impact[J]. Chinese Journal of Rock Mechanics and Engineering, 2008(S1): 2973-2977.

[12]	陈泰江, 向欣, 章广成. 基于线黏弹性接触的落石冲击地面参量理论研究[J]. 岩土力学. 2022, 43(9): 2410-2420.

[13]	Chen Taijiang, Xiang Xin, Zhang Guangcheng. Characteristic parameters theoretical analysis of rockfall impact on ground based on linear viscoelastic contact theory[J]. Rock and Soil Mechanics, 2022, 43(9): 2410-2420.

[14]	Labiouse V, Heidenreich B. Half-scale experimental study of rockfall impacts on sandy slopes[J]. Natural Hazards and Earth System Sciences, 2009, 9(6): 1981-1993.

[15]	Zhang G C, Tang H M, Xiang X, et al. Theoretical study of rockfall impacts based on logistic curves[J]. International Journal of Rock Mechanics and Mining Sciences, 2015, 78: 133-143.

[16]	Pichler B, Hellmich C, Mang H A. Impact of rocks onto gravel design and evaluation of experiments[J]. International Journal of Impact Engineering, 2005, 31(5): 559-578.

[17]	Kawahara S, Muro T. Effects of dry density and thickness of sandy soil on impact response due to rockfall[J]. Journal of Terramechanics, 2006, 43(3): 329-340.

[18]	陈颖骐, 王全才. 基于Hertz弹性理论和Thornton弹塑性假设的滚石冲击力的修正计算[J]. 科学技术与工程, 2018, 18(13): 37-41.

[19]	Chen Yingqi, Wang Quancai. Study on correction calculation of impact force of rockfall based on Hertz contact theory and Thornton elastoplasticity hypothesis[J]. Science Technology and Engineering, 2018, 18(13): 37-41.

[20]	刘成清, 陈玉满. 基于滚石运动参数的防护网数值模拟[J]. 科学技术与工程, 2020, 20(18): 7422-7428.

[21]	Liu Chengqing, Chen Yuman. Numerical simulation of protective net based on rolling stone motion parameters[J]. Science Technology and Engineering, 2020, 20(18): 7422-7428.

[22]	汪峰, 唐现梓, 黄伟. 山区输电塔滚石撞击响应分析及撞击力计算[J]. 科学技术与工程, 2022, 22(28): 12440-12448.

[23]	Wang Feng, Tang Xianzi, Huang Wei. Response analysis and simplified calculation of impact force of transmission tower in mountain area under rolling stone impact[J]. Science Technology and Engineering, 2022, 22(28): 12440-12448.

[24]	Li L, Lan H. Probabilistic modeling of rockfall trajectories: a review[J]. Bulletin of Engineering Geology and the Environment. 2015, 74(4): 1163-1176.

[25]	豆红强, 郑彬彬, 王浩, 等. 滚石冲击作用下考虑全植株形态的树木拦挡效应与能量吸收特性[J]. 振动与冲击. 2025, 44(15): 259-271.

[26]	Dou Hongqiang, Zheng Binbin, Wang Hao, et al. Tree blocking effect and energy absorption characteristics considering whole plant morphology under rolling stone impact[J]. Journal of Vibration and Shock, 2025, 44(15): 259-271.

[27]	Lu G, Ringenbach A, Caviezel A, et al. Mitigation effects of trees on rockfall hazards: does rock shape matter?[J]. Landslides, 2021, 18(1): 59-77.

[28]	Yan P, Zhang J, Kong X, et al. Numerical simulation of rockfall trajectory with consideration of arbitrary shapes of falling rocks and terrain[J]. Computers and Geotechnics, 2020, 122: 103511.

[29]	Masuya H, Amanuma K, Nishikawa Y, et al. Basic rockfall simulation with consideration of vegetation and application to protection measure[J]. Natural Hazards and Earth System Sciences, 2009, 9(6): 1835-1843.

[30]	Thorton C. Coefficient of restitution for collinear collisions of elastic-perfectly plastic spheres[J]. Transaction of ASME Journal of Applied Mechanics, 1997, 64(2): 383.

[31]	刘建辉, 柏亚双, 王兴宇, 等. 木材非线性本构模型研究进展[J]. 木材科学与技术, 2023, 37(1): 18-24.

[32]	Liu Jianhui, Bai Yashuang, Wang Xingyu, et al. Research progress on nonlinear constitutive model for wood[J]. Chinese Journal of Wood Science and Technology, 2023, 37(1): 18-24.

[33]	Murray Y D, Reid J D, Faller R K, et al. Evaluation of LS-DYNA wood material model 143: FHWA-HRT-04-096[R]. New York: FHWA, 2005.

[34]	周侃侃, 徐漫平, 郭飞燕, 等. 木荷木材物理力学性质及其加工性能研究[J]. 浙江林业科技, 2009(3): 19-23.

[35]	Zhou Kankan, Xu Manping, Guo Feiyan, et al. Study on physical and mechanical properties and machinability of schima superba wood[J]. Journal of Zhejiang Forestry Science and Technology, 2009(3): 19-23.