脆性细长结构受弯曲主导载荷作用发生碎裂的现象屡见不鲜，例如，建筑梁在爆炸和冲击载荷作用下会出现多处弯曲断裂；在微细观尺度上，脆性纤维在横向载荷作用下经常断裂成多段。为了理解这些场景中破坏的产生与发展，对弯曲卸载波的传播过程及弯曲脆断机理的研究必不可少。许多实验现象表明，即使在准静态载荷作用下，脆性细长梁在突然断裂时也可能引发次生断裂，典型的事例最初由Feynman等^[1]在观察意大利面的弯曲断裂现象时发现，由此产生“为什么意大利面条总是断成许多段”的问题，后称作Feynman问题。该问题涉及弯曲卸载波传播和脆断机理两方面。对于弯曲波传播问题，Schindler等^[2]基于Euler-Bernoulli梁理论给出了一个突发断裂导致卸载弯曲应力波的解析解，发现纯弯曲梁的突然断裂将导致毗邻区域的弯矩出现1.43倍初始弯矩的过冲，并通过实验测量到了因弯曲断裂而激发出的卸载弯曲应力波。Audoly等^[3]着重于Feynman问题，基于Euler-Bernoulli梁理论再次推导出初始断裂将导致邻近区域的曲率将超过初始值的结果，并解释了意大利面总是断裂成多段的现象。虽然Feynman问题可以通过Euler-Bernoulli梁理论得到合理的解释，但存在两点不足：(1) Euler-Bernoulli梁理论不能完整说明弯曲波的传播，在该理论框架下，一旦弹性梁突然断裂，所激发的卸载弯曲波将瞬间影响整个梁，基于半无限长梁的自相似解虽然说明了局部曲率高于断裂时刻曲率的过冲现象，但无法进一步研究二次断裂相关特征参数；(2) 采用断裂点弯矩突降，忽略了断裂发生的过程和断裂时间，也无法进一步分析断裂韧性对卸载波和二次断裂的影响。

相较于Euler-Bernoulli简单梁理论，Timoshenko梁理论^[4]进一步考虑了旋转惯性和切应力效应的影响，是一种更符合实际的梁理论。由于Timoshenko梁理论中弯曲波波速为有限值，该理论包含了弯曲波传播的特征时间。龙龙等 ^[5]在Timoshenko梁理论框架下研究了半无限长脆性梁发生突然弯曲断裂和斜坡断裂条件下所产生的卸载弯曲应力波传播问题，其分析结果来自数值反变换技术，计算精度有限。当采用积分变换方法求解，即以卷积公式写出在阶跃弯矩和斜坡弯矩作用下的弯曲波传播问题的解析表达式，并进一步分析半无限长Timoshenko梁的突然卸载问题，发现在主断裂点邻近区域存在一个特定位置，其峰值弯矩幅值的过冲达到1.68。

积分变换方法虽然对一些问题可以给出完整的解析解，但通常只适用于线弹性波，而对广泛存在的非线性波动问题以及波动与断裂相互耦合问题，其应用场景有限^[6]。另一方面，延时卸载的弯矩边界条件虽然为梁的断裂引入了历时过程，但将断裂视为一个已经预定的过程，其延时时间、卸载路径缺乏清晰的物理含义。李凤云^[7]采用基于拉伸-分离规律的内聚力单元模拟了脆性梁在纯弯曲载荷作用下的断裂过程，发现在一个广泛的材料参数和弯曲转动率下，梁的断裂点残余弯矩随断口转角单调变化，这意味着可以采用一个普适的内聚力弯曲断裂模型，唯象地描述脆性梁弯曲断裂过程中残余弯矩与断口开裂角度之间的关系。由于内聚力弯曲断裂规律通常为非线性，对于这一类具有复杂断裂模型的问题，积分变换方法无法使用。

Leonard等^[8]最早将特征线法引入Timoshenko梁的计算；Al-Mousawi等^[9-10]基于Timoshenko梁理论运用特征线法计算了变截面梁上弯曲波的传播，并论证了算法的稳定性和收敛性，以及Timoshenko梁理论在研究弯曲波问题时的准确性。以广义特征理论为基础的特征线计算方法，其优点是物理概念和物理图像十分清晰，可以清楚地揭示波传播的物理过程，差分格式稳定、计算效率高。本文中采用数值方法分析Timoshenko梁中弯曲波的传播过程，以及弯曲波传播与内聚力断裂之间的相互作用。计算方法以基本控制方程的特征线分析为基础。在Timoshenko梁理论框架下，研究脆性梁突然弯曲断裂所产生的卸载弯曲波传播问题，采用特征线方法得到卸载弯曲波传播的响应数值解。利用Timoshenko梁受冲击载荷作用下的解析解，验证特征线法的有限差分数值解的有效性；采用考虑断裂能的内聚力弯曲断裂模型，即断裂点所承受的残余弯矩（内聚力弯矩）是转角的线性递减函数，计算弯曲波传播的数值解，分析不同断裂参数下卸载弯曲波的传播特征，给出脆性梁突然弯曲断裂时二次断裂碎片尺寸的数值预测。

1.   Timoshenko梁理论和特征线法处理

相较于通过Euler-Bernoulli梁对弯曲波问题进行研究，同时考虑了旋转惯性和切应力效应的Timoshenko梁理论更适合描述瞬态响应问题。其中，旋转惯性的引入使Timoshenko梁中卸载弯曲波的传播具有强烈的局部化效应；而切应力效应的引入进一步降低了Timoshenko梁中弯曲扰动的传播，在梁中形成以剪切速度传播的扰动。因此，Timoshenko梁可以提供弯曲断裂的一个特征空间尺度，方便预测二次断裂。

1.1   Timoshenko梁基本方程

Timoshenko梁的理论模型如图1^[4]所示，其中$M$、$Q$分别表示弯矩和剪力，是梁所受的广义载荷；$\psi$、$\gamma$分别表示弯矩、剪力引起的截面转角。此外，为了后续表述方便清晰，规定如下符号：${w}_{M}$、${w}_{Q}$分别表示弯矩、剪力引起的截面横向位移（$y$轴方向），$w$表示总的横向位移；$v$、$\omega$分别表示截面的横向移动速度、转动角速度；$\kappa$表示梁的曲率。由动力学分析可知,以上10个量均为独立变量$x$、$t$的函数，$x$、$t$分别表示梁的轴向坐标、时间。

图  1  Timoshenko梁理论模型

Figure  1.  Model of Timoshenko beam theory

下载: 全尺寸图片 幻灯片

通过动力学、运动学和几何关系分析，Timoshenko梁的控制方程可写成如下仅包含2个未知函数$w$和$\psi$的形式 ^[4]：

以上方程中涉及的材料参数或结构常数有：体积密度$\rho$、弹性模量$E$、剪切模量$G$、剪切修正系数$k$、横截面面积$A$、横截面对中性轴的惯性矩$I$，以及在无量纲化中会使用的回转半径$R=\sqrt{I/A}$；涉及的本构方程包括：$M=EI\kappa$，$Q=GAk\gamma$。

1.2   特征线及其相容关系

采用特征线法^[8]处理Timoshenko梁的基本方程，将其处理为如下矩阵格式：

式中：${{\boldsymbol{W}}}=\left[\begin{array}{c}\begin{array}{c}M\\ \omega \\ Q\end{array}\\ v\end{array}\right]$，${\boldsymbol{B}}=\left[\begin{array}{cccc}0& EI& 0& 0\\ \dfrac{1}{\rho I}& 0& 0& 0\\ 0& 0& 0& -GAk\\ 0& 0& -\dfrac{1}{\rho A}& 0\end{array}\right]$，${\boldsymbol{b}}=\left[\begin{array}{c}\begin{array}{c}0\\ \dfrac{1}{\rho I}Q\\ -GAk\omega \end{array}\\ 0\end{array}\right]$。

在给定梁的参数矩阵$\boldsymbol{B}$情况下，未知函数$\boldsymbol{W}$可通过求解一阶波动方程式(2)获得。由于该方程涉及相互耦合的未知函数，对矩阵$\boldsymbol{B}$进行对角化处理，首先寻找$\boldsymbol{B}$的特征值$\mu$和相应的左特征行向量${\boldsymbol{l}}^{\mathrm{T}}$，使得：

上述存在非零的${\boldsymbol{l}}^{\mathrm{T}}$解的条件是：$\mathrm{det}\left({\boldsymbol{B}}-\mu {\boldsymbol{I}}\right)=0$，因此，得到4个特征值以及相应的4个左特征向量，分别为：

式中：${c}_{0}=\sqrt{\dfrac{E}{\rho }}$, ${c}_{Q}=\sqrt{\dfrac{Gk}{\rho }}$，分别表示纵波波速和等效剪切波速。

用常数的左特征向量${\boldsymbol{l}}_{i}^{\mathrm{T}}(i=1,2,3,4)$左点乘式(2)，注意到${\boldsymbol{l}}_{i}^{\mathrm{T}}{\boldsymbol{B}}={\mu }_{i}{\boldsymbol{l}}_{i}^{\mathrm{T}}$，得到：

这是4个完全解耦的以${\boldsymbol{l}}_{i}^{\mathrm{T}}\cdot {\boldsymbol{W}}$为未知函数的一阶波动方程，每个方程的特征解为：

式中只含有沿方向$\mu$的全导数$\dfrac{\mathrm{d}{\boldsymbol{W}}}{\mathrm{d}t}$。

由此可知，对于Timoshenko梁问题，存在4条特征线，沿着每条特征线，特定基本函数的组合满足相容关系，即：

称传播速度为纵波波速$\pm {c}_{0}$的特征线为主特征线，传播速度为剪切波速$\pm {c}_{Q}$的特征线为次特征线。上述Timoshenko梁方程的特征线及其相容关系的结果与文献[8]中所给一致。对于式(2)的一阶半线性偏微分方程组，其特征线走向与解无关，也与初始条件和边界条件无关，可在求解之前预先确定，并形成数值计算网格。

1.3   无量纲化

选取弹性模量E、密度ρ和截面回转半径R为3个基本特征量，则特征速度为杆波速c₀，特征时间为R/c₀（纵波横向传播回转半径的距离所需时间），特征弯矩为ER³，特征剪切应力为ER²。采用如下参数：

对独立变量和未知函数进行无量纲化，得到简化的特征线及相容关系：

式中：$\lambda$为等效剪切波速${c}_{Q}$与杆中纵波波速${c}_{0}$之比，$\lambda =\dfrac{{c}_{Q}}{{c}_{0}}=\sqrt{\dfrac{Gk}{E}}=\sqrt{\dfrac{k}{2\left(1+\nu \right)}}，\lambda {\text{＜}} 1$；$\nu$为泊松比；$\overline{A}$为梁的横截面形状系数。后文中的计算分析中，将统一使用无量纲数据，定义：$M$为无量纲弯矩，$Q$为无量纲剪力，$v$为无量纲横向移动速度，$\omega$为无量纲转动角速度，$x$为无量纲轴向坐标，$t$为无量纲时间，$T$为无量纲的当前时刻，$A$为无量纲横截面面积。

1.4   内点差分格式

使用沿着特征线上的有限差分格式计算时间空间解析区域的未知函数，差分网格如图2所示，空间步长$\mathrm{\Delta }x$和时间步长$\mathrm{\Delta }t$之间满足主特征线关系：$\mathrm{\Delta }x=\mathrm{\Delta }t$，次特征线$\mathrm{\Delta }x=\lambda \mathrm{\Delta }t$在网格内部。对沿着4条特征线$\dfrac{\mathrm{d}x}{\mathrm{d}t}=\mu \left(\mu =\pm 1, \pm \lambda \right)$的相容关系进行一阶精度差分，给出每个时刻$T+\mathrm{\Delta }t$内部节点位置$P$的4个未知函数与上一时刻T的4个点L、l、r、R的物理量的关系：

图  2  差分网格示意图

Figure  2.  Differential grid

下载: 全尺寸图片 幻灯片

式中：$\mathrm{\Delta }t$表示一个时间步长。令$f$指代未知函数$\mathit{W}={\left[M,\omega ,Q,v\right]}^{\mathrm{T}}$中的任一物理量，f_L、f_O、f_R表示左侧、中心和右侧节点上的量，f_l、f_r表示左侧和右侧的计算辅助点，即次特征线与等时刻线的交点。辅助点的物理量由邻近的2个主内点物理量进行线性插值得到：

式中的线性插值的权重由特征线传播的速度比决定，即$\lambda =\dfrac{{c}_{Q}}{{c}_{0}}$。

由式(10)可知，主特征线$\mu =\pm 1$控制了弯矩$M$和旋转角速度$\omega$之间的线性增量关系，并受积分路径$LP$或$RP$上剪应力${Q}_{L}^{(T)}$或${Q}_{R}^{(T)}$的影响；而次特征线$\mu =\pm \lambda$控制了剪切应力$Q$和平动速度$v$之间的线性增量关系，并受积分路径$lP$或$rP$上旋转角速度${\omega }_{l}^{(T)}$或${\omega }_{r}^{(T)}$的影响。差分方程组(11)的左侧为特征线相容关系式(9)的精确积分，而右侧项为弯矩-剪力耦合项的向前差分格式。该差分格式为1阶精度，联立辅助插值关系式(11)，可以写出以$T$时刻L、O和R点物理量表达的$T+\Delta t$时刻$P$点的未知函数。参考图2(b)，首先定义$T$时刻L和R点物理量的平均和以及平均差：

相应地，根据式(11)，l和r点物理量的平均和以及平均差为：

利用矩阵代数方法直接求解$T+\Delta t$时刻$P$点未知函数的显式表达式，首先将式(10)写作矩阵格式：

式中：

对矩阵$\boldsymbol{A}$求逆得：${\boldsymbol{A}}^{-1}=\left[\begin{array}{cccc}\dfrac{1}{2}& \dfrac{1}{2}& 0& 0\\ \dfrac{1}{2A}& -\dfrac{1}{2A}& 0& 0\\ 0& 0& \dfrac{1}{2}& \dfrac{1}{2}\\ 0& 0& -\dfrac{1}{2A\lambda }& \dfrac{1}{2A\lambda }\end{array}\right]$。

${\boldsymbol{W}}^{(T+\mathrm{\Delta }t)}={\boldsymbol{A}}^{-1}{\boldsymbol{a}}=\left[\begin{array}{c}\dfrac{1}{2}\left({M}_{L}^{(T)}+{M}_{R}^{(T)}\right)-\dfrac{A}{2}\left({\omega }_{R}^{(T)}-{\omega }_{L}^{(T)}\right)-\dfrac{1}{2}\left({Q}_{R}^{(T)}-{Q}_{L}^{(T)}\right)\Delta t\\-\dfrac{1}{2A}\left({M}_{R}^{(T)}-{M}_{L}^{(T)}\right)+\dfrac{1}{2}\left({\omega }_{R}^{(T)}+{\omega }_{L}^{(T)}\right)+\dfrac{1}{2A}\left({Q}_{R}^{(T)}+{Q}_{L}^{(T)}\right)\Delta t\\\dfrac{1}{2}\left({Q}_{r}^{(T)}+{Q}_{l}^{(T)}\right)+\dfrac{1}{2}\lambda A\left({v}_{r}^{(T)}-{v}_{l}^{(T)}\right)- \dfrac{1}{2}{\lambda }^{2}A\left({\omega }_{r}^{(T)}+{\omega }_{l}^{(T)}\right)\Delta t\\\dfrac{1}{2A\lambda }\left({Q}_{r}^{(T)}-{Q}_{l}^{(T)}\right)+\dfrac{1}{2}\left({v}_{r}^{(T)}+{v}_{l}^{(T)}\right)-\dfrac{1}{2}\lambda \left({\omega }_{r}^{(T)}-{\omega }_{l}^{(T)}\right)\Delta t\end{array}\right]$。

1.5   边界点差分格式

数值计算在有限解析区域：$0{\text{≤}} x{\text{≤}} L,t{\text{≥}} 0$内进行。$t=0$时刻各个函数数值由初始条件确定，在任意计算时刻$T+\Delta t$，网格内点根据式(14)直接计算；而边界点的函数值需要利用给定的边界条件、以及边界点物理量与邻接点$T$时刻物理量之间的2个特征线相容关系得出。图3给出左边界点${P}_{L}$和右边界点${P}_{R}$的计算网格。与内部节点类似，边界点的辅助点函数值采用主点插值：

图  3  边界网格示意图

Figure  3.  Boundary differential grid

下载: 全尺寸图片 幻灯片

1.5.1   左边界处理

对于左边界点${P}_{L}$，沿着2条向左传播的特征线$R{P}_{L}$和$r{P}_{L}$（分别对应特征值$\mu =-1,-\lambda$），存在2个相容关系等式，写成矩阵格式：

式中：${\boldsymbol{W}}^{{P}_{L}}={\left[\begin{array}{c}\begin{array}{c}M\\ \omega \\ Q\end{array}\\ v\end{array}\right]}^{{P}_{L}}$为左边界点未知物理量，${\boldsymbol{A}}_{\mathrm{2,4}}=\left[\begin{array}{cccc}1& -A& 0& 0\\ 0& 0& 1& \lambda A\end{array}\right]$为式(15)中$\boldsymbol{A}$矩阵的第2、4行组成的子矩阵，${\boldsymbol{a}}_{\mathrm{2,4}}=\left[\begin{array}{c}{M}_{R}^{(T)}-A{\omega }_{R}^{(T)}-{Q}_{R}^{(T)}\mathrm{\Delta }t\\ {Q}_{r}^{(T)}+\lambda A{v}_{r}^{(T)}-{\lambda }^{2}A{\omega }_{r}^{(T)}\mathrm{\Delta }t\end{array}\right]$为式(15)中$\boldsymbol{a}$向量的第2、4元素组成的子向量。上述矩阵条件不足以求解4个未知物理量，需要增加相应的边界条件完成计算。通常，在给定边界上规定了2个载荷或者位移条件。

考虑梁左端$\left(x=0\right)$弯矩为0，受突加速度${v}_{0}$作用，即${M}^{{P}_{L}}={M}_{{O}_{L}}^{(T)}=0；{v}^{{P}_{L}}={v}_{{O}_{L}}^{(T)}={v}_{0}$，将此条件代入式(17)，可得边界剪切力和角速度计算公式：

对于其他可能的随时间变化的边界条件，只要给出${\boldsymbol{W}}^{{P}_{L}}$中的任意2个未知量或者未知量之间的约束关系，即可算出边界点上的所有未知量。

1.5.2   右边界处理

对于右边界点${P}_{R}$，沿着2条向右传播的特征线$L{P}_{R}$和$r{P}_{R}$（分别对应特征值$\mu =1,\lambda$），存在2个相容关系等式，写成矩阵格式：

式中：${\boldsymbol{W}}^{{P}_{R}}={\left[\begin{array}{c}\begin{array}{c}M\\ \omega \\ Q\end{array}\\ v\end{array}\right]}^{{P}_{R}}$为右边界点未知物理量，${\boldsymbol{A}}_{\mathrm{1,3}}=\left[\begin{array}{cccc}1& A& 0& 0\\ 0& 0& 1& -\lambda A\end{array}\right]$为式(15)中$\boldsymbol{A}$矩阵的第1、3行组成的子矩阵，${\boldsymbol{a}}_{\mathrm{1,3}}=\left[\begin{array}{c}{M}_{L}^{(T)}+A{\omega }_{L}^{(T)}+{Q}_{L}^{(T)}\mathrm{\Delta }t\\ {Q}_{l}^{(T)}-\lambda A{v}_{l}^{(T)}-{\lambda }^{2}A{\omega }_{l}^{(T)}\mathrm{\Delta }t\end{array}\right]$为式(15)中$\boldsymbol{a}$向量的第1、3元素组成的子向量。

考虑右端为固支边界，即${v}^{{P}_{R}}={v}_{{O}_{R}}^{(T)}=0，{\omega }^{{P}_{R}}={\omega }_{{O}_{R}}^{(T)}=0$，代入到式(19)，可得到待求的2个边界载荷的显式表达：

对于其他可能的右边界条件，只要给出${\mathit{W}}^{{P}_{R}}$的4个未知量中任意2个未知量或者未知量之间的约束关系，即可算出边界点上的所有未知量。

2.   数值解及与理论解的对比

为了检验数值差分格式的可靠性，采用数值计算方法分析右端固支的静止不受力的Timoshenko梁，其左边界受突然冲击作用时的响应问题，如图4所示。关注梁的早期响应时，右侧边界尚未对梁的弯曲波造成影响。Boley等 ^[11]采用Laplace积分变换法给出了一系列特定边界条件下的解析解，部分已作为经典解收入专著中。在此选取文献[11]中的问题1：边界突加横向速度、无弯矩（图4(a)）这一条件进行计算。仿照Boley等 ^[11]的分析方法，采用积分变换方法也得到了左边界剪切力为0、突加弯矩（图4(b)）的解析解，在此作为第2个检验算例。

图  4  初始静止、不受力梁左端受突加载荷作用问题

Figure  4.  Left end of unstressed beam subjected to a sudden applied load with initial static state

下载: 全尺寸图片 幻灯片

2.1   静止梁端部受阶跃横向速度作用

选取Boley等^[11]提出的边界突加横向速度、无弯矩的条件进行计算：

式中；$H\left(t\right)$为单位阶跃函数。边界条件示意图如图4(a)所示，计算过程中取波速比$\lambda =\dfrac{{c}_{Q}}{{c}_{0}}=\dfrac{5}{9}$，突加横向速度${v}_{0}=1$，时间步长$\mathrm{\Delta }t=0.001$与空间步长$\mathrm{\Delta }x=\mathrm{\Delta }t=0.001$；计算结果与解析解的对比如图5～7所示。

图  5  本文与文献[11]在$x=0, 5$位置的剪力时程对比图

Figure  5.  Variation with time of the shear force at $x=0, 5$ compared with analytical solution in reference [11]

下载: 全尺寸图片 幻灯片

图  6  本文与文献[11]在$t=5, 10$时刻的剪力波形对比

Figure  6.  Variation with position of the shear force at $t=5, 10$ compared with analytical solution in reference [11]

下载: 全尺寸图片 幻灯片

图  7  本文与文献[11]在$t=5$时刻的速度和弯矩波形对比

Figure  7.  Variation with position of the velocity and the bending moment at $t=5$ compared with analytical solution in reference [11]

下载: 全尺寸图片 幻灯片

在横向速度冲击下，梁的剪切变形占主导地位，分析时更关注剪力峰值和剪切波的传播。图5给出了在左边界x = 0和梁上x = 5位置剪力随时间的变化。比较2个位置的剪力方向，会发现：(1)在边界上为了维持横向速度为常数，边界剪切力逐渐减小；(2)剪切力扰动最早在t = 5时刻到达x = 5位置，这是通过弯矩形成耦合效应，以弹性纵波速度c₀ = 1传播的扰动，且到达时刻剪切力方向是反向的；(3)在x = 5处，可以看到2个清晰的波阵面先后到达此处，首先是以纵波波速c₀为主导的边界扰动，在t = 9时刻突加的边界速度扰动以等效剪切波速c_Q = 5/9传播至此，导致剪力发生突变，这符合剪力和横向速度共同沿次特征线走的计算过程。

图5中所示的虚线是文献[11]中给出的解析解结果。数值计算曲线（实线）的整体走势与解析解（虚线）保持一致。图6给出了在$t=5$时刻和$t=10$时刻剪力的空间分布曲线，其中图6中的虚线是文献[11]中给出的解析解结果。同样可以看出，数值解基本重复了解析解的曲线走向，但在剪力突变的强间断位置，即解析解曲线的陡峭处会出现一定程度的光滑作用，平滑了$t=5,\ 9$时的曲线跳跃。造成这个现象的原因是采用了主特征线网格，对于该网格，CFL (Courant-Friedrichs-Lewy)条件^[12]中的CFL数$\dfrac{{c}_{0}\mathrm{\Delta }t}{\mathrm{\Delta }x}=1$，而剪力传播是沿着由主特征线网格插值而来的次特征线进行的，插值导致CFL数$\dfrac{{c}_{Q}\mathrm{\Delta }t}{\mathrm{\Delta }x} {\text{＜}} 1$。根据CFL条件，CFL数等于1时，数值结果无耗散、无振荡；CFL数小于1时，虽然数值格式的计算稳定性得以保证，但存在由于数值扩散机制而出现的耗散现象。

图7给出了在$t=5$时刻梁上横向速度和弯矩的波形分布，图中速度强间断以恒定的跳跃幅值传播至$x={c}_{Q}t=\dfrac{25}{9}$处，期间速度幅值由于弯曲波传播过程中的色散而产生衰减，对照弯矩的变化，此处的弯矩也存在一个弱间断。

2.2   静止梁端部受阶跃弯矩作用

不同于横向冲击下关注剪力峰值和剪切波的传播，在研究弯曲变形主导下的弯曲断裂时，弯矩峰值和弯曲波的传播更值得关注。考虑图4(b)所示的静止梁端部受突加弯矩作用，左边界条件为：

式中：$H\left(t\right)$为Heaviside函数。计算过程中取波速比$\lambda =\dfrac{{c}_{Q}}{{c}_{0}}=\dfrac{5}{9}$，突加弯矩${m}_{0}=1$，空间步长与时间步长$\mathrm{\Delta }x=\mathrm{\Delta }t=0.001$。在此初边值条件下，运用特征线法求得t=5, 10, 15, 20时刻的弯矩分布波形如图8实线所示，对比给出了解析解结果，如图8虚线所示。可以看出，数值计算结果与解析解完全一致，特征线法准确地捕捉到了强间断波阵面的传播过程，且没有发生数值衰减，这是因为计算过程中弯矩是沿着主特征线传播的，CFL数$\dfrac{{c}_{0}\mathrm{\Delta }t}{\mathrm{\Delta }x}=1$，数值格式无耗散、无振荡。

图  8  不受力梁边界受突加弯矩作用时的弯曲波传播对比

Figure  8.  Bending wave propagations in a Timoshenko beam subjected to a suddenly applied boundary moment compared with analytical solution

下载: 全尺寸图片 幻灯片

在边界突加弯矩作用下，梁的弯曲变形占主导地位，分析时更关注弯矩峰值和弯曲波的传播过程。从图8可以看到，突加的恒定弯矩并没有得到稳定传播，而是出现了极速衰减，甚至出现反向弯矩，由此形成了一个陡峭的三角波阵面，造成这种弯矩快速下降的原因在于剪切力的耦合作用，图9给出了$t=5,10,15,20$时刻的剪切力分布波形。以$t=15$时刻（蓝线）的波形图为例，波速为${c}_{0}=1$的纵波强间断波阵面到达$x=15$时，导致弯矩出现1的峰值，对应的剪切力分布是一个弱间断波阵面，在$x=\dfrac{75}{9}=8.33$位置以横波速度${c}_{Q}=5/9$传播的剪切扰动到达，剪力快速拉升，相应地梁中的弯矩分布也逐渐提高，并在边界点达到1的边界值。由此可见，弯矩和剪力这一对耦合的波在梁的行进中形成了此消彼长的关系。

图  9  静止梁边界受突加弯矩作用时的剪力传播对比

Figure  9.  Shear wave propagations in a Timoshenko beam subjected to a suddenly applied boundary moment compared with analytical solution

下载: 全尺寸图片 幻灯片

上述问题可以转换为意大利面条在准静态纯弯曲作用下折断的问题：考虑一根以恒定速率弯曲（缓慢增加曲率）的脆性梁，在${T}_{0}$时刻发生断裂，此时梁承受的弯矩达到极限值${m}_{0}$，梁一分为二。在断口中心，由于左右对称性，剪切应力始终为0，而弯矩则瞬时下降。考虑$T {\text{＞}} {T}_{0}$时刻断裂激发的卸载弯曲波的传播特征，以$t=T-{T}_{0}$为计算时刻，则梁的初始条件和左端边界条件为：

对此问题进行数值模拟，所得到的不同时刻梁中的弯矩分布图等价于自由静止梁突加弯矩m₀问题的解，即图8(a)所描述的弯矩分布的反相和平移，如图10(a)所示（取${m}_{0}=1.0$）。从图10(a)中结果可以看出，对于Timoshenko梁来说，一旦出现突然断裂，则距离断裂点一段距离后，梁中的弯矩才可能出现过冲$\left({M {\text{＞}} m}_{0}\right)$。图10(b)描绘了一旦发生断裂，卸载弯曲波扫过的不同位置所经历的最大弯矩分布。传统的Euler-Bernoulli梁突然卸载问题的解^[2-3]预测，断口附近任意位置的弯矩过冲均为1.43，而Timoshenko梁的解与此不一样，与断裂点距离不同的位置上产生的过冲峰值不一样，距离越远，过冲值先增后减，特别是在$x=18.01$位置，最大过冲值达到初始值的1.68倍。通过理论分析已经发现上述结论，本文中的数值模拟工作进一步证实此结果。

图  10  初始静止处于纯弯曲状态的梁突然断裂时与卸载波相关的弯矩分布

Figure  10.  Bending moment profiles in a Timoshenko beam suddenly broken under pure bending

下载: 全尺寸图片 幻灯片

3.   Timonshenko梁的内聚力弯曲断裂

从上述计算结果可以看出，考虑切应力效应和旋转惯性的Timoshenko梁理论解决了Euler-Bernoulli简单梁理论中弯曲波速度无限大的问题，揭示了弯曲波以有限的纵向波速${c}_{0}$和等效剪切波速${c}_{Q}$传播的特征。另一方面，材料断裂过程也需要一定的时间，断裂时间长短与材料的断裂韧性以及裂纹与外部载荷相互作用有关。对于在极短时间发生的动态脆断，后者往往以卸载波形式发生。本节把弯曲脆断过程视为一个断口抗弯刚度与转角的内聚力断裂过程，采用特征线分析方法研究卸载弯曲波传播、断裂过程以及断裂时间的相互影响。

3.1   脆性梁的内聚力断裂模型

大部分梁的弯曲断裂是一个横向断口从拉伸面向压缩面传播的过程，直到裂纹面切断整个梁，如图11所示。李凤云^[7]在模拟脆性梁弯曲断裂过程中发现断裂中的梁的残余弯曲强度（即可承受弯矩M）与断口开裂角度$\theta$之间存在一致的函数关系。如果把M视为$\theta$的函数，则该函数通常为递减的曲线关系，反映了断口的内聚力断裂特征。而$M\left(\theta \right)$曲线下方的面积为梁完全破断所消耗的能量。只要给定$M\left(\theta \right)$模型和基本参数，则可以将该条件作为耦合边界条件施加到上述数值计算中，从而计算和分析梁的断裂过程和断裂所激发的卸载弯曲波传播规律。

图  11  耦合裂纹开裂角度与所受弯矩的内聚力断裂模型

Figure  11.  Cohesive fracture model coupling bending moment and cracking opening angle

下载: 全尺寸图片 幻灯片

内聚力弯曲断裂模型参考Kipp等^[13]提出的针对韧性金属的线性内聚力拉伸断裂模型，此模型下，将残余弯矩和开裂角度进行积分，即得断裂过程所消耗的断裂能，具有明确的物理意义。该模型包含裂缝附近的应力释放所导致的弯曲波的传播特性以及裂纹生长过程中的局部能量耗散。

图12给出了3种不同断裂能下的线性弯曲断裂模型，保持断裂弯矩（脆性梁所能承受的最大弯矩）不变，通过控制开裂角度的大小间接控制断裂能，开裂角度越大表示断裂所需的断裂能越大。由于数值计算的需要，计算时，初始时刻的残余弯矩取断裂弯矩（脆性梁所能承受的最大弯矩）的99%，为了保持每种情况下断裂能总体不变，断裂结束时的开裂角度相应增加。参考李凤云^[7]计算的结果，断裂弯矩取${M}_{\mathrm{c}}=0.07$；假设断裂时裂纹对称，则计算时开裂角度取${\theta }_{\mathrm{c}}$的一半。

图  12  不同断裂能下的线性弯曲断裂模型

Figure  12.  Model of linear flexural fracture at different fracture energies

下载: 全尺寸图片 幻灯片

3.2   内聚力断裂过程和卸载弯曲波

计算如图12所示3种不同边界条件下断裂点的残余弯矩$M$随时间$t$的变化，结果如图13所示，每条弯矩时程曲线均经历了初始的缓慢卸载和后期的快速卸载阶段，$M\left(t\right)$曲线类似于抛物线形状；与之线性相关的裂纹开裂角度$\theta \left(t\right)$同样有着初始缓慢增加、后期快速开裂的过程，这与Heisser等^[14]用高速摄影所观察到的意大利面弯曲断裂过程相似，即裂纹在起始后相当长一段时间（数毫秒）保持静止，随后以接近声速的高速传播约10 μs，最终导致断口完全断裂。随着临界开裂角度的减小，断口弯矩完全卸载的时间也相应缩短，也就是说断裂能越小，材料越脆，断裂过程越快。按此趋势，当开裂角度趋近于零时，弯矩的卸载过程等同于突降弯矩的情况。

图  13  断裂点（边界点）的弯矩$M\left(t\right)$时程曲线

Figure  13.  History curves of bending moment $M\left(t\right)$ at the break point (boundary)

下载: 全尺寸图片 幻灯片

图14所示为不同断裂能下特征时间$t=$5, 15, 25, 35时刻的弯矩波形（纵坐标表示归一化的弯矩），显示了内聚力弯曲断裂模型下卸载弯曲波的传播过程。图中波头的瞬时到达位置在数值上与时间一一对应，即如在$t=5$时波头传至$x=5$处。对比同一时刻不同断裂能下的弯矩峰值可知，断裂能越小，弯曲断裂所激发的卸载弯曲波的弯矩波峰值越大，与Euler-Bernoulli 梁自相似解所得到的峰值弯矩1.43相比，Timoshenko梁断裂所激发的峰值弯矩会超过该数值，在$t=35$时刻，开裂角度由小到大对应的归一化弯矩的峰值分别为1.64、1.65、1.63。

图  14  不同断裂能下各个时刻的弯矩波形

Figure  14.  Bending wave propagations at different fracture energies

下载: 全尺寸图片 幻灯片

3.3   二次断裂位置预测

为了预测初始断裂后由于卸载弯曲波导致的二次断裂点位置，需要找到弯矩过冲峰值出现的位置，通过计算足够长的时间，计算得到了4种临界断裂角：${\theta }_{\mathrm{c}}=$0.004，0.008，0.012，0.040下的归一化峰值弯矩${M}_{\mathrm{m}\mathrm{a}\mathrm{x}}/{m}_{0}$的空间分布包络线，如图15所示。4种材料参数所对应的断裂能分别为${M}_{\mathrm{c}}{\theta }_{\mathrm{c}}/2=1.4\times {10}^{-4},2.8\times {10}^{-4}, 4.2\times {10}^{-4}, 1.4\times {10}^{-3}$。 从各点可能经历的峰值弯矩曲线可以看出以下趋势。

图  15  不同断裂能下断裂导致邻近区域弯矩过冲峰值的空间分布包络线

Figure  15.  Envelope of the overshot bending moment in a brittle beam when it breaks with different cohesive fracture energy

下载: 全尺寸图片 幻灯片

(1) 4种断裂能所对应的归一化峰值弯矩首次大于1的位置分别出现在${x}_{0.004}=4.95, {x}_{0.008}= 5.17, {x}_{0.012}=5.47, {x}_{0.04}=5.80$。在这个距离以内，二次断裂不会发生。一般而言，对于一根均匀的脆性杆，一次断裂发生点邻近${x}_{0}=5～6$以内的区域不大可能发生二次断裂。

(2) 超过${x}_{0}$距离，脆性杆将经历超过原始弯矩的弯矩过冲。对于断裂能越小（越脆）的材料，邻近区域所达到的过冲峰值弯矩越大，在某些特殊位置，峰值弯矩达到最大值。表1记录了前3种脆性材料杆的最大弯矩${M}_{\mathrm{p}}$、发生位置${x}_{\mathrm{c}}$以及发生时刻${t}_{\mathrm{c}}$。对于这3种非常脆的材料，Timoshenko梁初次断裂所激发的峰值弯矩（1.64～1.67）超过经典理论预测的1.43，所发生位置位于17.7附近，这些最大峰值弯矩发生点是最容易发生二次断裂的位置。断裂能越小，梁上峰值弯矩极值出现的时间越早，出现位置越靠前，但相近的断裂能下极值出现时间和位置变化不大，极值位置集中出现在$x=18$附近。

表  1  不同断裂能（开裂角度）下峰值弯矩极值的时空坐标

Table  1.  Positions and time of peak bending moment extremes at different fracture energies (cracking angles)

开裂角度 θc	归一化峰值弯矩极值 Mp = Mmax/m0	峰值弯矩 位置xc	峰值弯矩 时间tc
0.004	1.668	17.67	32.95
0.008	1.657	17.68	34.33
0.012	1.639	17.72	35.90

下载: 导出CSV 

| 显示表格

(3) 对于前3种非常脆的材料，材料断裂能对最大峰值弯矩的数值和发生位置影响不大。如果材料的脆性显著降低、即材料断裂能明显加大，将得到经典的基于Euler-Bernoulli梁理论所求解的卸载弯曲波过冲自相似解，即过冲峰值弯矩趋近于1.43。为了更好地观察断裂能大小对卸载弯曲波和二次断裂的影响趋势，在图15中给出了断裂能提高一个量级（开裂角度${\theta }_{\mathrm{c}}=0.040$）后的峰值弯矩${M}_{\mathrm{m}\mathrm{a}\mathrm{x}}/{m}_{0}$的空间分布曲线，此时弯矩过冲的最大值下降到接近1.43，发生位置向更远处移动。

(4) 值得注意的是，在保持临界弯矩不变的条件下，过大临界开裂角度意味着断裂能更大，因此，完全断裂时间延长（如图13所示），从而弱化了Timoshenko梁理论中的转动惯性效应和剪切应力效应，使其行为更接近Euler-Bernoulli梁。从图15结果可以看出，材料韧性增加导致断口断裂时间延长，所产生的弯曲卸载波的过冲位置更远，幅值更接近于Euler-Bernoulli梁1.43倍的预测结果。这个结果也从侧面说明了Timoshenko梁理论更适合描述瞬态响应问题。

综合上述结果，可以预见如果材料越脆，则初次断裂所激发的卸载弯曲波越强烈，Timoshenko梁局部的旋转惯性和切应力效应对所激发弯曲波的影响变大。

对于一根缺陷随机分布的脆性杆而言，首次弯曲断裂出现的地方一定是缺陷最严重的位置，根据上述卸载弯曲波的传播特点，二次断裂出现的位置不可能出现在$x {\text{＜}}5$的区域，因为在此区域弯矩不超过首次断裂的临界弯矩。在$x{\text{＞}} 5$的区域内，首次断裂所激发的卸载弯曲波将导致弯矩过冲，均存在二次断裂的可能。对于非常脆的材料，最大弯矩过冲发生在$x=18$的位置，其数值超过1.63，此位置发生二次断裂的可能性最大。

为直观显示数据的物理意义，取Heisser等^[14]提出的意大利面的材料参数：弹性模量$E=3.8\;\mathrm{G}\mathrm{P}\mathrm{a}$，密度$\rho =1.5\;\mathrm{g}/{\mathrm{c}\mathrm{m}}^{3}$，纵波波速${c}_{0}= 1\;592\;\mathrm{m}/\mathrm{s}$；梁的几何截面取边长为1 mm的正方形，则回转半径$R=\dfrac{1}{2 \sqrt{3}}$ mm。此时二次断裂碎片尺寸的下限数值为1.44 mm $(x {\text{＞}} 5)$，最可能发生二次断裂的位置在5.20 mm$(x=18)$附近。表1中参数有量纲的结果如表2所示。

表  2  意大利面断裂峰值弯矩极值的时空坐标

Table  2.  Positions and time of peak bending moment extremes of spaghetti

开裂角度θc/rad	断裂能Gc/μJ	峰值弯矩位置xc/mm	峰值弯矩时间tc/μs
0.004	12.8	5.10	5.97
0.008	25.6	5.10	6.23
0.012	38.4	5.12	6.51

下载: 导出CSV 

| 显示表格

针对意大利面二次断裂的问题，龙龙^[15]采用高速摄像机记录了2种类型的意大利面的断裂过程，并对断裂碎片进行搜集和统计，碎片尺寸集中在6～13倍的梁截面厚度；基于内聚力单元模型进行的数值模拟结果揭示，二次断裂点介于5～12倍梁截面厚度区域。在本文的模型分析中，如果采用瞬时断裂模型，即完全忽略材料的断裂韧性，则最大峰值弯矩发生在5.2倍梁厚度位置($x=18$)，显著优于传统的Euler-Bernoulli梁自相似结果。如果考虑材料的断裂韧性，则二次断裂点的弯矩峰值位置将进一步远离一次断裂点，如图15所示。当前，材料的弯曲断裂韧性具体数据尚有待测定，这是今后工作的重点。

4.   结　论

本文中基于Timoshenko梁理论，运用特征线法精确计算分析了半无限长脆性梁在准静态弯曲加载至断裂时刻，断裂所激发的卸载弯曲波的传播过程与特点。通过对比突加弯矩和弯矩延时卸载这2种模型下的特征线法数值解与积分变换解析解，初步验证了特征线法在计算Timoshenko梁瞬态断裂问题上的有效性、准确性和便捷性；通过将唯象的耦合残余弯矩和开裂角度的内聚力弯曲断裂模型引入到Timoshenko梁的瞬态波动分析中，使脆性梁的突然断裂有了更加清晰的物理含义，在计算分析此断裂机制下卸载弯曲波的传播过程后，给出了无量纲化的二次断裂碎片尺寸的下限数值为$x {\text{＞}} 5$，最可能发生二次断裂的位置在5.2倍梁厚度位置($x=18$)附近，此处卸载弯曲波所产生的弯矩过冲超过1.67倍的初始弯矩。上述结果显示了Timoshenko梁在卸载弯曲波激发过程中强烈的局部化特征，也为确定二次断裂碎片尺寸的预测提供了依据。