北京大学学报自然科学版   2016, Vol. 52 Issue(4): 687-691

文章信息

薛纭, 翁德玮
XUE Yun, WENG Dewei
Kirchhoff弹性杆不连续量的奇异函数表达
Discontinuous Quantities of Kirchhoff Elastic Rod Expressed by Singularity Function
北京大学学报(自然科学版), 2016, 52(4): 687-691
Acta Scientiarum Naturalium Universitatis Pekinensis, 2016, 52(4): 687-691

文章历史

收稿日期: 2015-10-10
修回日期: 2016-03-16
网络出版日期: 2016-07-12
Kirchhoff弹性杆不连续量的奇异函数表达
薛纭1, 翁德玮2     
1. 上海应用技术大学机械工程学院, 上海 201418;
2. 上海大学, 上海市应用数学和力学研究所, 上海 200072
摘要: 弹性细杆静力学和动力学的Kirchhoff方程要求在外力、质量几何以及本构方程的间断或不光滑点处分段表达, 这不利于数值计算。根据计算梁弯曲变形的奇异函数法, 将奇异函数引入Kirchhoff方程, 将弹性杆分段定义的量拓展为沿全杆的连续函数。借助Mathematica软件, 对存在侧向集中载荷的弹性杆进行数值模拟, 结果表明, 引入奇异函数可以避免分段导致的繁琐计算, 提高计算效率。
关键词: Kirchhoff弹性杆     不连续量     奇异函数     局部载荷     平衡位形    
Discontinuous Quantities of Kirchhoff Elastic Rod Expressed by Singularity Function
XUE Yun1, WENG Dewei2     
1. School of Mechanical Engineering, Shanghai Institute of Technology, Shanghai 201418;
2. Shanghai Institute of Applied Mathematics and Mechanics, Shanghai University, Shanghai 200072
Corresponding author: XUE Yun, E-mail: xy@sit.edu.cn
Abstract: Kirchhoff equation of thin elastic rod statics and dynamics must be written in piecewise at section with discontinuous or nonsmooth quantities such as external forces, mass geometry and physical parameters, which leads to inconvenience to the numerical calculation. According to singular function method in calculation of beam bending deformation, these discontinuous or nonsmooth quantities of the rod are expressed by singular function and become continues quantities alone centerline of the rod. Numerical simulations of equilibrium configuration of the rod acted by lateral concentrated load are made by means of Mathematics software. Results explain that introducing singular function to express discontinuous or nonsmooth quantities can avoid complicated calculation and improve the computational efficiency.
Key words: Kirchhoff elastic rod     discontinuous quantities     singularity function     local loads     equilibrium configuration    

弹性细杆静力学的一般理论由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的函数, e3为切向基矢量, 指向弧坐标增加方向。两组基的关系为${\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个位形坐标关于st为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)

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

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

其中, ω=(ω1 ω2 ω3)T, Ω=(Ω1 Ω2 Ω3)T, ωi=ω·ei, Ωi=Ω·ei, 矩阵$\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)

其中, 假定杆的原始形态为直杆, Miωi为截面内力的主矩和弯扭度的主轴分量; B1B2为关于主轴xy的抗弯刚度, B3为关于主轴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}$, I1I2为截面对主轴xy的惯性矩, I3为对z轴的极惯性矩, 且有${I_3}={I_1} + {I_2}$

方程(8)要求质量几何参数ρ, Jj, Ij, 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弹性杆的几何和物理参数时, 自变量xxi用弧坐标ssi表示。

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

考虑长为l的弹性细杆, 分为n段, 段长从端面开始依次为li, 其中${l_1} +... + {l_n}=l$。弹性细杆分段连续量包括几何量和物理量。如阶梯杆, 其截面尺寸、截面积、截面对主轴的二次矩和转动惯量在每一杆段上为常值; 再如, 分段连续的分布载荷, 等等。分段连续量用Zi表示, 借助单位阶跃函数(式(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)

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

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

弹性细杆的集中量包括集中质量和集中载荷, 例如, 弹性杆的某一阶梯段长度极短可看做集中质量, 作用于弹性杆侧向的集中力和集中力偶等都是集中载荷。集中量用Yi表示, 借助脉冲函数(式(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个同一性质的集中量Yi沿弹性杆轴线的分布密度为

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

这里, Yi可以为矢量。

对于作用在${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)$处作用的集中力或集中力偶Pi, 则集度为qi=Pi/Δ(Δ为小量), 用单位阶跃函数表示为

$ 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}$, 在弹性杆弧坐标为s1s2的截面处分别作用集中力f=fe1和集中力偶m=m2e2, 它们都是随动载荷。计算弹性杆在起始值条件下的平衡位形, 并与无随动载荷情形比较。

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)

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

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

1) λ=0.5, b=0.02, s1=0.3, s2=0.6, s∈(0, 3), Δ=0.01; F1(0)=0, F2(0)=0, F3(0)=2, ψ (0)=0, θ (0)=0.2, φ3(3)=0, ψ' (0)=2, φ'(0)=0.5。

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

图 1. 集中载荷对弹性杆位形的影响 Figure 1. Effect of concentrated load on the equilibrium configuration of the elastic rod

2) λ=0.5, b=0.01, s1=0.3, s2=0.6, s∈(0, p), Δ=0.01; F1(0)=0, F2(0)=0, F3(0)=2, ψ (0)=0, θ(0)=π/2, φ3(0)=0, θ'(0)=0, φ'(0)=2。

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

图 2. 集中载荷以及起始扭率对弹性杆位形的影响 Figure 2. Effect of concentrated load and starting twist on the equilibrium configuration of the elastic rod

5 结语

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

参考文献
[1] Kirchhoff G. Ueber das Gleichgewicht und die Bewegung eines unendlich duennen elastischen Stabes. J Rein Angew Math , 1859, 56 : 285–313 .
[2] Kirchhoff G. Vorlesungen üeber Mathematische Physic: Mechanik. Leipzig: B G Teubner, 1876 .
[3] Clebsch A. Théorie der elasticit t fester k rper. Leipzig: B G Teubner, 1682 .
[4] Clebsch A. Théorie de l'Elasticité des corps solides. New York: Johnson Reprint Corp, 1966 .
[5] Love A E H. A treatise on mathematical theory of elasticity. Cambridge: Cambridge Univ Press, 1892 .
[6] Love A E H. A treatise on mathematical theory of elasticity. 4th ed. New York: Dover, 1927 .
[7] Dill E H. Kirchhoff's theory of rods. Archive for History of Exact Science , 1992, 44 (1) : 1–23 DOI:10.1007/BF00379680 .
[8] 刘延柱. 弹性细杆的非线性力学: DNA力学模型的理论基础. 北京: 清华大学出版社, 2006 .
[9] Svetlitsky V A. Dynamics of rods. Berlin: Springer, 2005 .
[10] Xue Yun, Liu Yanzhu, Chen Liqun. The Schr dinger equation for a Kirchhoff elastic rod with noncircular cross section. Chinese Physics , 2004, 13 (6) : 794–797 DOI:10.1088/1009-1963/13/6/002 .
[11] Xue Yun, Shang Huilin. Jourdain principle of a super-thin elastic rod dynamics. Chinese Physics Letters , 2009, 26 (7) : 074501–3 DOI:10.1088/0256-307X/26/7/074501 .
[12] 薛纭, 翁德玮, 陈立群. 精确Cosserat弹性杆动力学的分析力学方法. 物理学报 , 2013, 62 (4) : 044601.
[13] 刘延柱, 薛纭. 受圆柱面约束螺旋杆伸展为直杆的动力学分析. 力学学报 , 2011, 43 (6) : 1151–1156.
[14] 薛纭, 曲佳乐, 陈立群. Cosserat生长弹性杆动力学的Gauss最小拘束原理. 应用数学和力学 , 2015, 36 (7) : 700–709.
[15] Moulton D E, Lessinnes T, Goriely A. Morphoelastic rods. Part Ⅰ: a single growing elastic rod. Journal of the Mechanics and Physics of Solids , 2013, 61 : 398–427 DOI:10.1016/j.jmps.2012.09.017 .
[16] 刘鸿文. 材料力学Ⅱ. 4版. 北京: 高等教育出版社, 2004 : 12 .
[17] 王燮山. 奇异函数及其在力学中的应用. 北京: 科学出版社, 1993 .
[18] Macauley W H. Note on deflection of beams. Mess Math , 1919, 48 : 129–130 .
[19] Pilkey W D. Clebsch's method for beam deflections. J Eng Edu , 1964, 54 (5) : 170–173 .
[20] 冷坳坳, 杨庆超, 楼京俊, 等. 基于奇异函数的船舶推进轴系校中计算研究. 武汉理工大学学报:交通科学与工程版 , 2015, 39 (2) : 401–404.
[21] 吴阿林. 连续阶梯梁影响线方程的奇异函数法研究. 力学与实践 , 2004, 26 (2) : 70–72.
[22] 徐彬, 梁启智. 框架结构分析的奇异函数方法. 力学与实践 , 2001, 23 (4) : 47–49.
[23] 刘延柱. 高等动力学. 北京: 高等教育出版社, 2001 .