弹性细杆静力学的一般理论由Kirchhoff^[1-2]、Clebsch^[3-4], Love^[5-6]以及Dill^[7]的工作形成。在忽略弹性杆的伸缩和截面的剪切变形的条件下(Kirchhoff弹性杆模型), 用弹性杆中心线的Frenet轴系随弧坐标的运动表达弹性杆的位形, 导出的以弧坐标为自变量, 以曲率和挠率为未知函数的平衡微分方程与刚体动力学的Euler方程, 在数学形式上完全相同, 由此产生“Kirchhoff动力学比拟”思想, 为连续弹性细杆的离散化提供了新的研究思路和方法^[8-15], 开创了用刚体动力学的概念和方法研究弹性细杆静力学的新思路。将Cosserat有向介质理论中的方向子取代Frenet轴系作为截面主轴坐标系的基矢量^[9], 此方向子随弧坐标的运动同样形成弹性杆的位形, 以弯扭度的主轴分量为未知函数列出的平衡微分方程仍具有Euler方程的形式, 在体现“Kirchhoff动力学比拟”思想的同时, 可以在概念和方法上进行更广泛的动力学比拟, 还可以方便地考虑弹性杆存在伸缩和截面剪切变形的一般情况, 称为Cosserat弹性杆精确模型。

Kirchhoff方程要求外力、质量几何以及本构方程沿杆长连续。实际上, 间断或不光滑点的存在是不可避免的, Kirchhoff方程只能分段表达, 导致边界条件增加, 使得计算繁琐。在求解梁的弯曲变形中引入奇异函数, 很好地解决了此问题^[16-17]。

采用间断多项式表述梁的弯矩和挠度方程由Clebsch^[3]提出, Macauley^[18]建议用括号〈〉表示。然而, 他们的建议未引起注意。直到Dirac提出δ函数, Schwartz阐明了δ函数并创立广义函数论后, Pilkey^[19]首次用于求解梁的变形。奇异函数法是将集中力或集中力偶, 以及局部分布力, 用奇异函数统一表达成沿全梁的分布量, 从而避免分段计算的麻烦。王燮山^[17]系统阐述了奇异函数及其在材料力学、高等材料力学、弹性薄板等力学中的应用, 冷坳坳等^[20]用奇异函数研究了船舶推进轴系的变形问题, 吴阿林^[21]用奇异函数给出单跨梁的影响线表达式和连续阶梯梁影响线方程, 徐彬等^[22]讨论了框架结构分析的奇异函数方法。

Kirchhoff弹性杆的特点是受空间力系作用, 位形的空间形态复杂, 因此, 将奇异函数法引入Kirchhoff弹性杆力学很有必要。这样, 可以为涉及的分段或局部定义的几何和物理参数给出一个沿杆长连续的表达式, 给数值计算带来方便。

1 Kirchhoff动力学方程

依据Kirchhoff弹性杆力学理论, 以杆的截面为对象, 建立惯性坐标系$O{\rm{-}}\xi \eta \zeta $以及与截面固结的形心主轴坐标系$P{\rm{-}}xyz$, 沿坐标轴的单位基矢量分别为${\boldsymbol{e}^I}={(\begin{array}{*{20}{c}} {{\boldsymbol{e}_\xi }} & {{\boldsymbol{e}_\eta }} & {{\boldsymbol{e}_\zeta }}\end{array})^{\rm{T}}}$和${\boldsymbol{e}^P}={(\begin{array}{*{20}{c}}{{\boldsymbol{e}_1}} & {{\boldsymbol{e}_2}} & {{\boldsymbol{e}_3}}\end{array})^{\rm{T}}}$, 其中${\boldsymbol{e}_i}(i=1, \; 2, \; 3)$为中心线弧坐标s和时间t的函数, e₃为切向基矢量, 指向弧坐标增加方向。两组基的关系为${\boldsymbol{e}^{\rm{P}}}=\boldsymbol{H}{\boldsymbol{e}^{\rm{I}}}$, 式中H为单位正交阵。$P{\rm{-}}xyz$的位置用截面形心相对惯性系的矢径$\boldsymbol{r} = \overline {OP} $的坐标阵$\boldsymbol{r}={(\begin{array}{*{20}{c}}{\xi (s, t)} & {\eta (s, t)} & {\zeta (s, t)}\end{array})^{\rm{T}}}$和截面姿态的Euler角列阵$\boldsymbol{q}={(\begin{array}{*{20}{c}}\psi & \vartheta & \varphi \end{array})^{\rm{T}}}$描述。截面的运动方程为

$ \boldsymbol{r} = \boldsymbol{r}(s, \;t)\;, \;\;\;\boldsymbol{q} = \boldsymbol{q}(s, \;t)。 $

(1)

设6个位形坐标关于s和t为2阶连续可微。对于除端部外不受约束的自由弹性杆, 6个广义坐标为独立变量, 根据Kirchhoff假定, 其偏导需满足方程^[8]:

$ \boldsymbol{r}' = {\boldsymbol{e}_{\rm{3}}}, $

(2)

式中, 撇号表示对弧坐标s的偏导数。式(2)的投影式为

$ \xi ' = \sin \psi \sin \vartheta, \eta ' =-\cos \psi \sin \vartheta, \zeta ' = \cos \vartheta。 $

(3)

方程(3)是不可积的, 构成内约束而无需约束力。弯扭度ω(s, t)和角速度Ω(s, t)用Euler角表示为

$ \boldsymbol{\omega }= \sum\limits_{i = 1}^3 {{\Xi _i}{{q'}_i}}, \;\;\boldsymbol{\Omega } = \sum\limits_{i = 1}^3 {{\Xi _i}{{\dot q}_i}}, $

(4)

其中, q_i依次为3个Euler角, ${\Xi _i}$为Euler角的矢值函数, 变量顶部的点号表示对t的偏导数。式(4)在截面主轴坐标系中的矩阵式为

$ \boldsymbol{\omega } = \boldsymbol{\Theta q}', \;\;\;\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} = \boldsymbol{\Theta \dot q}, $

(5)

其中, ω=(ω₁ ω₂ ω₃)^T, Ω=(Ω₁ Ω₂ Ω₃)^T, ω_i=ω·e_i, Ω_i=Ω·e_i, 矩阵$\mathit{\boldsymbol{ \boldsymbol{\varTheta} }}$的定义^[8]为

$ \mathit{\boldsymbol{ \boldsymbol{\varTheta} }} = \left( {\begin{array}{*{20}{c}} {\sin \vartheta \sin \phi }&{\cos \phi }&0\\ {\sin \vartheta \cos \phi }&{-\sin \phi }&0\\ {\cos \vartheta }&0&1 \end{array}} \right) $

(6)

假设弹性细杆服从虎克定律, 本构关系表示为

$ {M_i} = {B_i}{\omega _i}\;\quad (i = 1, \;2, \;3), $

(7)

其中, 假定杆的原始形态为直杆, M_i和ω_i为截面内力的主矩和弯扭度的主轴分量; B₁和B₂为关于主轴x和y的抗弯刚度, B₃为关于主轴z的抗扭刚度。弹性细杆动力学方程为

$ {\partial _s}\boldsymbol{F}-\rho A{\partial _t}\boldsymbol{v} + \boldsymbol{f} = \boldsymbol{0}, $

(8a)

$ {\partial _s}\boldsymbol{M} + {\partial _s}\boldsymbol{r} \times \boldsymbol{F}-{\partial _t}(\boldsymbol{J} \cdot \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}) + \boldsymbol{m} = \boldsymbol{0}, $

(8b)

其中, F为截面内力的主矢; ρ为杆的密度; A为截面积; $\boldsymbol{v}=\boldsymbol{\dot r}$; J为截面对质心的惯量并矢, 在主轴坐标系下的坐标阵为$\boldsymbol{J}={\rm{diag (}}\begin{array}{*{20}{c}} {{J_{\rm{1}}}} & {{J_2}} & {{J_3}} \end{array}{\rm{)}}$, 可表示为${J_j}=\rho {I_j}$, I₁和I₂为截面对主轴x和y的惯性矩, I₃为对z轴的极惯性矩, 且有${I_3}={I_1} + {I_2}$。

方程(8)要求质量几何参数ρ, J_j, I_j, A以及外力f, m沿杆长都是连续的。当出现集中质量、集中力和集中力偶以及局部分布载荷等时, 这些参数对弧坐标不连续, 需分段列出动力学方程。求解分段连续微分方程的工作量较大, 用奇异函数表达这些量就可以避免分段列写和求解动力学方程。

2 奇异函数的定义及其微分和积分

奇异函数的定义^[17]为

$ f(x) = {\left\langle {x-{x_i}} \right\rangle ^n}, $

(9a)

$ {\left\langle {x-{x_i}} \right\rangle ^n} = \left\{ {\begin{array}{*{20}{l}} {{{(x-{x_i})}^n}, \;\;\;x \ge {x_i}}\\ {0, \;\;\;\;\;\;\;\;\;\;\;\;\;x < {x_i}} \end{array}} \right.(n \ge 0), $

(9b)

$ {\left\langle {x-{x_i}} \right\rangle ^n} = \left\{ {\begin{array}{*{20}{c}} {\infty, \quad x = {x_i}}\\ {0, \quad x \ne {x_i}} \end{array}} \right.(n < 0), $

(9c)

式中, n=-1时称为d函数、脉冲函数或Dirac函数, n=0时称为单位阶跃函数或Heaviside函数, 奇异函数的微分和积分公式为

$ \frac{{\rm{d}}}{{{\rm{d}}x}}{\left\langle {x-{x_i}} \right\rangle ^n} = \left\{ {\begin{array}{*{20}{l}} {{{\left\langle {x-{x_i}} \right\rangle }^{n-1}}, \;\;\;\;\;\;\;\;\;\;\;\;\;\;n \le 0, }\\ {n{{\left\langle {x - {x_i}} \right\rangle }^{n - 1}}, \;\;\;\;\;\;\;\;\;\;\;\;n > 0, } \end{array}} \right. $

(10)

$ \int_{-\infty }^x {{{\left\langle {x-{x_i}} \right\rangle }^n}} {\rm{d}}x = \left\{ {\begin{array}{*{20}{l}} {{{\left\langle {x-{x_i}} \right\rangle }^{n + 1}}, \;\;\;\;\;\;\;\;\;\;n \le 0, }\\ {\frac{1}{{n + 1}}{{\left\langle {x - {x_i}} \right\rangle }^{n + 1}}, \;\;\;n > 0\;} \end{array}} \right. $

(11)

Kirchhoff弹性杆是以弧坐标为自变量, 因此, 将奇异函数用于表达空间形态的Kirchhoff弹性杆的几何和物理参数时, 自变量x和x_i用弧坐标s和s_i表示。

3 弹性细杆不连续量沿杆长的表达 3.1 弹性细杆分段连续量沿杆长的表达

考虑长为l的弹性细杆, 分为n段, 段长从端面开始依次为l_i, 其中${l_1} +... + {l_n}=l$。弹性细杆分段连续量包括几何量和物理量。如阶梯杆, 其截面尺寸、截面积、截面对主轴的二次矩和转动惯量在每一杆段上为常值; 再如, 分段连续的分布载荷, 等等。分段连续量用Z_i表示, 借助单位阶跃函数(式(11)中n=0), 将定义在弹性杆上的同一性质的分段连续量拓展到全杆, 统一表示为

$ Z = {Z_1}{\left\langle s \right\rangle ^0} + \sum\limits_{i = 2}^n {\left( {{Z_i}-{Z_{i-1}}} \right)} {\left\langle {s-\sum\limits_{j = 1}^{i - 1} {{l_j}} } \right\rangle ^0}, $

(12)

这里, Z_i亦可为矢量, 例如分布力。方向不变的常值矢量要指明是相对惯性参照系$O{\rm{-}}\xi \eta \zeta $还是与截面固结的形心主轴坐标系$P{\rm{-}}xyz$(后者称为随动矢量)。

3.2 弹性细杆的集中量沿杆长的表达

弹性细杆的集中量包括集中质量和集中载荷, 例如, 弹性杆的某一阶梯段长度极短可看做集中质量, 作用于弹性杆侧向的集中力和集中力偶等都是集中载荷。集中量用Y_i表示, 借助脉冲函数(式(9)中n=-1), 将定义在弹性杆上一点的集中量拓展到全杆。集中量沿弹性杆轴线的分布密度为

$ {y_i} = \mathop {\lim }\limits_{s \to {s_i}} \frac{{{Y_i}}}{{s-{s_i}}} = {Y_i}\left\langle {s-{s_i}} \right\rangle {}^{-1}, $

(13)

则弹性杆上n个同一性质的集中量Y_i沿弹性杆轴线的分布密度为

$ y = \sum\limits_{i = 1}^n {{Y_i}\left\langle {s-{s_i}} \right\rangle {}^{-1}}, $

(14)

这里, Y_i可以为矢量。

对于作用在${s_i} \in (0, \; \; l)$处的集中力偶M, 可以考虑在$\left[{{s_i}-\Delta, \; {s_i}} \right]$和$\left[{{s_i}, {s_i} + \Delta } \right]$上作用有反向平行的均布力$M/{\mathit{\Delta }^2}$。容易证明, 反向平行均布力的合成即为集中力偶M。下面, 将此定义在2Δ上的反向平行均布力拓展到全杆上。由式(12)得

$ \begin{array}{l} q(s) = \frac{M}{{{\Delta ^2}}}\left\{ {{{\left\langle {s-({s_i}-\Delta )} \right\rangle }^0}-{{\left\langle {s - {s_i}} \right\rangle }^0}} \right\} - \\ \;\;\;\;\;\;\;\;\frac{M}{{{\Delta ^2}}}\left\{ {{{\left\langle {s - {s_i}} \right\rangle }^0} - {{\left\langle {s - ({s_i} + \Delta )} \right\rangle }^0}} \right\}, \end{array} $

(15)

等号右边的两项可化为

$ \begin{array}{l} \frac{{{{\left\langle {s + \Delta-{s_i}} \right\rangle }^0}-{{\left\langle {s-{s_i}} \right\rangle }^0}}}{\Delta }{\rm{ = }}\frac{{\rm{d}}}{{{\rm{d}}s}}{\left\langle {s - {s_i}} \right\rangle ^0} = {\left\langle {s - {s_i}} \right\rangle ^{ - 1}}, \\ \frac{{{{\left\langle {s - {s_i}} \right\rangle }^0} - {{\left\langle {s - \left( {{s_i} + \Delta } \right)} \right\rangle }^0}}}{\Delta }\\ = \frac{{\rm{d}}}{{{\rm{d}}s}}{\left\langle {s - \Delta - {s_i}} \right\rangle ^0} = {\left\langle {s - \Delta - {s_i}} \right\rangle ^{ - 1}}, \end{array} $

两式相减得

$ \begin{array}{l} \frac{{{{\left\langle {s-{s_i}} \right\rangle }^{-1}}-{{\left\langle {s - \left( {{s_i} - \Delta } \right)} \right\rangle }^{ - 1}}}}{\Delta } = \frac{{\rm{d}}}{{{\rm{d}}s}}{\left\langle {s - \left( {{s_i} - \Delta } \right)} \right\rangle ^{ - 1}}\\ \quad \quad \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = {\left\langle {s - \left( {{s_i} - \Delta } \right)} \right\rangle ^{ - 2}}。 \end{array} $

令Δ→0, 得到集中力偶的线分布集度:

$ q(s) = M\frac{{\rm{d}}}{{{\rm{d}}s}}{\left\langle {s-{s_i}} \right\rangle ^0} = M{\left\langle {s-{s_i}} \right\rangle ^{-2}}。 $

(16)

数值计算时, 集中量的上述表述过于理想化, 并不适合数值计算。集中力是力分布区域极小的抽象, 因此将集中量还原为分布区域很小的分布量, 这样, 可用式(12)来表达在弹性杆弧坐标${s_i} \in (0, \; \; l)$处作用的集中力或集中力偶P_i, 则集度为q_i=P_i/Δ(Δ为小量), 用单位阶跃函数表示为

$ q = {q_i}\left( {{{\left\langle {s-{s_i}} \right\rangle }^0}-{{\left\langle {s-({s_i} + \Delta )} \right\rangle }^0}}。 \right) $

(17)

4 弹性杆平衡位形算例

半径为b的圆截面弹性杆长为l, 扭弯刚度之比为$\lambda = {B_3}{\rm{/}}{B_1}$, 在弹性杆弧坐标为s₁和s₂的截面处分别作用集中力f=fe₁和集中力偶m=m₂e₂, 它们都是随动载荷。计算弹性杆在起始值条件下的平衡位形, 并与无随动载荷情形比较。

令v=0, Ω=0, 式(8)即为Kirchhoff弹性杆的平衡微分方程(仍记为式(8))。下面将变量和参数无量纲化:

$ \begin{array}{l} \boldsymbol{F} = \frac{{{B_1}}}{{{l^2}}}\boldsymbol{\hat F}, \boldsymbol{\omega} = \frac{{\boldsymbol{\hat \omega }}}{l}, \boldsymbol{M} = \frac{{{B_1}}}{l}\boldsymbol{\hat M}, \lambda = \frac{{{B_3}}}{{{B_1}}}, b = \frac{{\hat b}}{l}, \\ s = l\hat s, \xi \;{\rm{ = }}\;l\hat \xi, \;\;\eta \;{\rm{ = }}\;l\hat \eta, \;\;\zeta \;{\rm{ = }}\;l\hat \zeta, \end{array} $

(18)

字母上方“^”表示无量纲量。随动载荷化为分布力并无量纲化为

$ \begin{array}{l} \boldsymbol{f} = \frac{{{B_1}}}{{{l^3}}}{{\hat f}^*}{\boldsymbol{e}_1}, \quad {{\hat f}^*} = \frac{{\hat f}}{{\hat \Delta }}\left( {{{\left\langle {s-{s_1}} \right\rangle }^0}-{{\left\langle {s-\left( {{s_1} + \Delta } \right)} \right\rangle }^0}} \right), \\ \boldsymbol{m} = \frac{{{B_1}}}{{{l^2}}}{{\boldsymbol{\hat m}}^*}{\boldsymbol{e}_2}, \quad {{\hat m}^*} = \frac{{\hat m}}{{\hat \Delta }}\left( {{{\left\langle {s - {s_1}} \right\rangle }^0} - {{\left\langle {s - \left( {{s_1} + \hat \Delta } \right)} \right\rangle }^0}} \right), \end{array} $

(19)

其中$\mathit{\hat \Delta }$为无量纲小量。将式(8a)和(8b)分别向Résal坐标系和形心主轴坐标系$P{\rm{-}}xyz$投影^[23], 得

$ {F_1}^\prime-(\varphi ' + \psi '\cos \theta ){F_2} + {F_3}\psi '\sin \theta + {\hat f^*}\cos \varphi = 0, $

(20a)

$ {F_2}^\prime + (\varphi ' + \psi '\cos \theta ){F_1}-{F_3}\theta ' + {\hat f^*}\sin \varphi = 0, $

(20b)

$ {F_3}^\prime-{F_1}\psi '\sin \theta + {F_2}\theta ' = 0, $

(20c)

$ \begin{array}{l} 2\theta ''-(1-\lambda ){{\psi '}^2}\sin 2\theta + 2\lambda \varphi '\psi '\sin \theta-2{F_2} - {{\hat m}^*}\sin \varphi \\ = 0, \end{array} $

(20d)

$ \begin{array}{l} \psi ''\sin \theta-\lambda \theta '(\varphi ' + \psi '\cos \theta ) + 2\theta '\psi '\cos \theta-{F_1} + {{\hat m}^*}\cos \varphi \\ = 0, \end{array} $

(20e)

$ \varphi ' + \psi '\cos \theta = {\omega _{30}}, $

(20f)

式中省略了所有无量纲记号, F_i为截面主矢在Résal坐标系的投影, ω₃₀为ω₃=ω ·e₃的起始值。

考虑如下参数和起始值, 借助Mathematica软件求数值解。

1) λ=0.5, b=0.02, s₁=0.3, s₂=0.6, s∈(0, 3), Δ=0.01; F₁(0)=0, F₂(0)=0, F₃(0)=2, ψ (0)=0, θ (0)=0.2, φ₃(3)=0, ψ' (0)=2, φ'(0)=0.5。

图 1显示载荷${\hat f^*}=20, \; {\hat m^*}=50$(杆1)以及无此载荷作用(杆2)时的弹性杆的位形, 表明集中力使弹性杆发生了显著的变形。

2) λ=0.5, b=0.01, s₁=0.3, s₂=0.6, s∈(0, p), Δ=0.01; F₁(0)=0, F₂(0)=0, F₃(0)=2, ψ (0)=0, θ(0)=π/2, φ₃(0)=0, θ'(0)=0, φ'(0)=2。

图 2显示:无集中载荷和φ'(0)=0时, 弹性杆为平面圆环(杆3);当φ'(0)=2时, 弹性杆不再是平面圆环(杆4);在集中力${\hat f^*}=20, \; \; {\hat m^*}=50$作用下, 且φ'(0)=2时弹性杆的位形(杆5)。

5 结语

Kirchhoff弹性细杆静力学和动力学方程要求物理和几何参数沿杆全长连续分布, 然而一般情形下, 连续的弹性细杆上常存在不连续的质量几何、作用力以及物理参数。利用奇异函数, 将这些不连续量拓展到全杆定义, 借助数学软件, 可以方便地进行数值模拟, 由此避免了分段计算带来的繁琐和工作量, 提高了计算效率。