引言

近年来，随着能源危机和环境问题的凸显，海上新能源的开发和利用倍受关注。风能作为清洁、可再生新能源，对改善能源结构和环境具有重要意义，已成为许多国家可持续发展战略的重要组成部分。我国在《能源发展“十三五”规划》中提出大力发展海上风电，引领能源发展进入创新驱动新阶段，建设清洁低碳、安全高效的现代能源体系。我国东南沿海海域辽阔，开发海上风电场可有效改善东部沿海地区电力供应问题，具有广阔发展前景。风机基础作为风电机组支撑体系，将上部荷载传递给地基，其作用尤为重要。单桩基础因具有结构简单、占用面积小、承载力高、沉降量小且均匀等特点，在国内外已建工程中占比近80%（Wang等，2018d）。

目前，国内外学者主要采用数值法和解析法对单桩式海上风机结构系统自振频率开展研究，在数值模拟方面，Zuo等（2018）和Wang等（2018a）采用非线性Winkler地基模型，考虑桩-土相互作用，通过数值模拟方法研究了单桩式海上风机结构系统自振频率及其在波浪荷载、风荷载作用下的动力响应，结果表明，桩-土相互作用对风机结构系统自振频率和动力响应的影响不可忽略；杨春宝等（2018）基于有阻尼系统运动方程，考虑桩-土相互作用和风机塔筒变截面特性，采用有限元法提出了单桩式海上风机结构系统自振频率求解方法，并对比实际工程验证该方法精度；Andersen等（2012）提出了求解黏土中单桩式风机系统一阶自振频率简化方法，采用非线性p-y曲线模型，考虑桩-土相互作用，基于有限差分法求解桩帽等效水平和旋转弹簧刚度，并进一步求解单桩式风机系统一阶自振频率。在解析理论方面，主要基于欧拉-伯努利梁理论，将风机塔筒和桩简化为连续等截面梁，将上部结构简化为集中质量，如Van der Tempel等（2002）给出了刚性地基条件下计算风机结构自振频率解析公式；Adhikari等（2011，2012）和Bhattacharya等（2011）将桩-土相互作用简化为泥面处施加水平弹簧和旋转弹簧，给出了柔性地基条件下计算风机结构自振频率的解析公式，并将理论计算值与模型试验结果进行对比，结果表明理论计算值较实测值偏大。在进一步的研究中，将风机塔筒和桩简化为分段等截面梁，并设耦合弹簧模拟水平弹簧和转动弹簧间的耦合效应，如Arany等（2015）推导出求解风机结构自振频率的解析方法；余璐庆等（2018）基于Arany等（2015）的研究，系统分析了基础刚度和轴向荷载等对结构前四阶频率的影响。

需指出的是，为方便计算，以上理论方法均将塔筒等效为等截面梁，未考虑塔筒变截面特性。为此，本文考虑水-桩-土相互作用、风机塔筒变截面特性及轴力和转动惯量等因素的影响，求解单桩式海上风机结构系统自振频率，旨在为单桩式海上风机结构设计和运营评估提供简便、快速、有效的系统频率计算方法。

1.   控制方程

将风机结构简化为图 1所示的悬臂梁，其中桩为等截面梁，塔筒为变截面梁，将塔顶机舱、轮毂和叶片等效为与悬臂梁刚性连接的集中质量，其转动惯量为J_r，将底部桩-土相互作用简化为线性弹簧，将水简化为附加质量。

图 1  风机结构简化模型

Figure 1.  Simplified turbine model

下载: 全尺寸图片 幻灯片

欧拉-伯努利梁横向自由振动运动方程如下：

式中，$\omega $为频率；z为结构沿竖向高度的空间坐标；$\phi (z)$为梁横向位移；${H_{\rm{s}}} = {h_{{\rm{f}}}} + {h_{\rm{p}}} + {h_{\rm{t}}}$，为结构总高度，其中h_f、h_p、h_t分别为桩在土中高度、桩在水中高度、塔筒高度；ρ_f、ρ_t分别为桩和塔筒密度；A_f、A_t分别为桩和塔筒截面面积；I_f、I_t分别为桩和塔筒截面惯性矩；E为结构弹性模量；$P$为海上风机轴向荷载；${m_{\rm{a}}}$为动水压力引起的附加质量；${k_{i, {{\rm{f}}}}}$为第i层土弹簧刚度。

考虑轴向荷载时，塔筒质量大于机舱、轮毂等上部质量（Bush等，2009），轴向荷载取值直接影响结构自振频率和动力响应。Jydsk（2002）给出了塔筒在高度z处轴向荷载表达式：

式中，${M_{\rm{r}}}$为上部结构集中质量；g为重力加速度。

桩-水相互作用产生的动水压力可视为部分水体质量与结构加速度的乘积，部分水体质量称为附加质量，桩-水相互作用对水中结构的地震响应具有较大影响，Wang等（2018b，2018c）基于刚性柱法解析解给出了水体附加质量简化公式，王丕光等（2019）针对柔性圆柱体在水中的振动问题，提出了代替水-结构动力相互作用的桩-水附加质量模型，该附加质量表达式为：

式中，${m_0}$为结构刚性运动引起的附加质量；${\rho _{\rm{w}}}$为水体密度；r为圆柱半径。

将底部桩-土相互作用简化为线性弹簧，Chiou等（2020）给出了不同土体参数下侧向弹簧刚度表达式：

式中，D为桩直径；G_{i, s}和${\nu _{i, {\rm{s}}}}$分别为第i层土剪切模量和泊松比；$\eta $为系数，取值为0.4—1.2。

桩底和塔顶边界条件为：

2.   频率方程

崔灿等（2012）将变截面梁等效为多段均匀梁，将塔筒单元离散为相互连接的若干段。根据土层分层特性将土层中的桩分为N₁段，水中的桩作为1段，将变截面塔筒分为N₂段，每段均视为等截面梁，即将整个结构划分为N（N=N₁+N₂+1）段，根据式（1），第i段梁单元运动方程为：

式中，$ E{I}_{i}=\left\{ \begin{array}{l}E{I}_{{\rm{f}}}{\rm{, }}{\rm{\hspace{0.17em}}}i\le {N}_{1}+1\\ \frac{1}{{l}_{i}}{\displaystyle {\int }_{{\rm{ }}{z}_{i+1}}^{{\rm{ }}{z}_{i}}E{I}_{\rm{t}}(z){\rm{d}}z}{\rm{, }}{\rm{\hspace{0.17em}}}i{\rm{ }}＞{N}_{1}+1\end{array} \right.$；$ {\rho }_{i}{I}_{i}=\left\{ \begin{array}{l}{\rho }_{{\rm{f}}}{I}_{{\rm{f}}}{\rm{, }}{\rm{\hspace{0.17em}}}i\le {N}_{1}+1\\ \frac{1}{{l}_{i}}{\displaystyle {\int }_{{z}_{i+1}}^{{z}_{i}}{\rho }_{\rm{t}}{I}_{\rm{t}}(z){\rm{d}}z}{\rm{, }}{\rm{\hspace{0.17em}}}i{\rm{ }}＞{N}_{1}+1\end{array} \right.$；$ {k}_{i}=\left\{ \begin{array}{l}{k}_{i, {\rm{f}}}{\rm{, }}{\rm{\hspace{0.17em}}}i\le {N}_{1}\\ 0{\rm{, }}{\rm{\hspace{0.17em}}}i{\rm{ }}＞{N}_{1}\end{array} \right.$；${m_i} = $ $ \left\{ \begin{array}{l}{\rho }_{{\rm{f}}}{A}_{{\rm{f}}}{\rm{, }}{\rm{\hspace{0.17em}}}i{\rm{ }}＜{N}_{1}\rm{ }\\ {\rho }_{{\rm{f}}}{A}_{{\rm{f}}}{\rm{+}}{m}_{\rm{a}}{\rm{, }}{\rm{\hspace{0.17em}}}{N}_{1}\le i\le {N}_{1}+1\\ \frac{1}{{l}_{i}}{\displaystyle {\int }_{{\rm{ }}{z}_{i+1}}^{{\rm{ }}{z}_{i}}{\rho }_{\rm{t}}{A}_{\rm{t}}(z)}{\rm{ d}}z{\rm{, }}{\rm{\hspace{0.17em}}}i{\rm{ }}＞{N}_{1}+1\end{array} \right.$

若${k_i} - {m_i}{\omega ^2} \leqslant 0$，求解特征方程得到的振型函数为：

式中，${\delta _i} = \sqrt {{{(0.25\nu _i^2 - \beta _i^{})}^{0.5}} + 0.5\nu _i^{}} $；${\varepsilon _i} = \sqrt {{{(0.25\nu _i^2 - \beta _i^{})}^{0.5}} - 0.5\nu _i^{}} $；$\nu _i^{} = (P + {\rho _i}{I_i}{\omega ^2})/E{I_i}$；$\beta _i^{} = ({k_i} - {m_i}{\omega ^2})/E{I_i}$；${l_i} = {z_i} - {z_{i - 1}}{{\rm{ }}}$。

若$ 0{\rm{ }}＜{\rm{ }}{k}_{i}-{m}_{i}{\omega }^{2}\le {(P+{\rho }_{i}{I}_{i}{\omega }^{2})}^{2}/4E{I}_{i}$，求解特征方程得到的振型函数为：

式中，${\delta _i} = \sqrt {{{(0.25\nu _i^2 - \beta _i^{})}^{0.5}} + 0.5\nu _i^{}} $；${\varepsilon _i} = \sqrt {0.5\nu _i^{} - {{(0.25\nu _i^2 - \beta _i^{})}^{0.5}}} $。

若$ {k}_{i}-{m}_{i}{\omega }^{2}{\rm{ }}＞{\rm{ }}{(P+{\rho }_{i}{I}_{i}{\omega }^{2})}^{2}/4E{I}_{i}$，求解特征方程得到的振型函数为：

式中，${\gamma _i} = (\sqrt {1 - {c_i}}){d_i}$；${\theta _i} = (\sqrt {1 + {c_i}}){d_i}$；${c_i} = (P + {\rho _i}{I_i}{\omega ^2})/\sqrt {4E{I_i}({k_i} - {m_i}{\omega ^2})} $；${d_i} = (({k_i} - {m_i}{\omega ^2})/$$4E{I_i}{)^{0.25}}$。

第i段梁和（i+1）段梁在连接点${z_i}$处的位移、转角、弯矩、剪力连续，可得以下关系：

将式（8）、式（9）、式（10）分别代入式（11），可得以下关系：

第i段梁和第（i+1）段梁待定系数分别为：

式中，系数${t_{ij}}$、${s_{ij}}$见附录。

由式（11）得${{\boldsymbol{A}}_{({\rm{N}})}}$和${{\boldsymbol{A}}_{(1)}}$关系为：

式中，${\boldsymbol{Z}} = {{\boldsymbol{Z}}_{({\rm{N}} - 1)}}{{\boldsymbol{Z}}_{({\rm{N}} - 2)}} \cdots {{\boldsymbol{Z}}_{(2)}}{{\boldsymbol{Z}}_{(1)}}$，${{\boldsymbol{Z}}_{(i)}} = {{\boldsymbol{T}}_{(i + 1)}}^{ - 1}{{\boldsymbol{T}}_{(i)}}$。

将式（8）、式（9）、式（10）代入式（6）中前两式底部边界条件，将式（8）代入式（12）再代入式（6）中后两式顶部边界条件，可得以下方程：

式中，$ {\boldsymbol{\varGamma }} = \left[ {\begin{array}{*{20}{c}} {{{\boldsymbol{\varGamma }} _1}} \\ {{{\boldsymbol{\varGamma }} _2}} \end{array}} \right]$。

当$ {k_i} - {m_i}{\omega ^2} \leqslant {(P + {\rho _i}{I_i}{\omega ^2})^2}/4E{I_i}$时，则有：

当$ {k}_{i}-{m}_{i}{\omega }^{2}{\rm{ }}＞{\rm{ }}{(P+{\rho }_{i}{I}_{i}{\omega }^{2})}^{2}/4E{I}_{i}$时，则有：

式中，${\chi _1}{\rm{ = }}\frac{{{\omega ^2}{J_{\rm{r}}}}}{{E{I_{\rm{N}}}}}$；${\chi _2}{\rm{ = }}\frac{{{M_{\rm{r}}}{\omega ^2}}}{{E{I_{\rm{N}}}}}$；${\chi _3}{\rm{ = }}\frac{{P + {\rho _{\rm{N}}}{I_{\rm{N}}}{\omega ^2}}}{{E{I_{\rm{N}}}}}$。

对于式（18），向量${{\boldsymbol{A}}_{(1)}} = {\{ {A_1}{{\rm{ }}}{B_1}{{\rm{ }}}{C_1}{{\rm{ }}}{D_1}\} ^{\rm{T}}}$的系数不能同时为零，而有非零解的条件是系数矩阵行列式为零，即：

根据式（23），利用MATLAB中fsolve函数可求解梁自振频率，代入式（8）、式（9）、式（10）可求解出相应的振型函数。

3.   自振频率分析

以NREL-5MW基准风机为研究对象，进行结构自振频率分析，参照Jonkman等（2009）给出的基准风机工程参数，其中风机塔筒高87.6m，桩长36m，水深30m，底端直径6m，壁厚0.027m，顶端直径3.87m，壁厚0.019m，上部质量总和约3.5×10⁵kg。基础为大直径单桩，截面属性与塔筒底部一致。钢材弹性模量为210GPa，密度为7850kg/m³，考虑油漆、焊接、法兰等的质量，将塔筒密度提高为8.5×10³kg/m³。

3.1   方法验证

将风机结构参数代入本文方法求解结构自振频率，同时在ABAQUS软件中建立梁单元模型进行方法验证。当分段数足够多时，结构自振频率应趋于稳定和精确，如图 2所示。当塔筒分段数增至5段后，一阶频率变化基本趋于稳定，增至10段后分段数对结构自振频率的影响较小。因此，考虑计算精度和效率，本文取N₂=10。

图 2  频率随分段数变化曲线

Figure 2.  The curve of frequency varing with the number of segments

下载: 全尺寸图片 幻灯片

将刚性地基条件下的自振频率解析解与数值模拟结果进行对比验证，结果见表 1。由表 1可知，解析解与数值模拟结果吻合良好，误差＜1%，前三阶频率均具有较好的精度。为方便参数定量分析，本文假定土体为均匀土层。弹簧刚度取100MPa时，考虑桩-土相互作用得到的解析解与数值模拟结果的误差见表 2。由表 2可知，本文方法求解结果与数值模拟结果基本一致，一阶频率精度较好，误差为1%左右，三阶频率误差＜3%。将本文解析解与Arany等（2015）推导的Welney 1S 3.6MW风机结构自振频率进行对比，结果见表 3。由表 3可知，一阶自振频率和实测频率误差＜0.5%，达到了较好的精度，而Arany等（2015）计算结果误差较大，可能因为未考虑塔筒变截面特性。

表 1  NREL-5MW基准风机刚性地基模型验证

Table 1.  Verification of NREL-5MW rigid foundation model

频率验证	频率/rad·s-1
	一阶	二阶	三阶
解析解	5.09	23.93	59.94
数值解	5.12	24.14	60.25
误差/%	0.54	0.85	0.52

下载: 导出CSV 

| 显示表格

表 2  NREL-5MW基准风机桩-土相互作用模型验证

Table 2.  Verification of NREL-5MW pile-soil interaction model

频率验证	频率/rad·s-1
	一阶	二阶	三阶
解析解	1.456	10.963	29.210
数值解	1.443	10.725	28.493
误差/%	0.902	2.213	2.515

下载: 导出CSV 

| 显示表格

表 3  Welney 1S 3.6MW风机桩-土相互作用模型一阶自振频率对比

Table 3.  Verification of Welney 1S 3.6MW pile-soil interaction model

频率验证	实测值	Arany等（2015）	解析解
频率/rad·s-1	0.350	0.331	0.3484
误差/%	—	5.429	0.457

下载: 导出CSV 

| 显示表格

综上所述，本文解析计算方法具有较高精度，计算公式简洁，便于实际应用，可有效替代有限元数值模拟复杂的建模与参数分析过程。

3.2   影响因素分析

为方便研究桩基础埋深、上部质量、转动惯量、桩-水相互作用对结构一阶自振频率的影响，本文定义无量纲量${l_0}$、${\sigma _0}$、${\xi _0}$：

式中，$a$为桩径，取值为3m；M_r、J_r按实际工程取为3.5×10⁵kg和3.6×10⁶kg·m²；h、M、J分别为基础埋深、上部质量和转动惯量的变量。

3.2.1   桩基础埋深的影响

不同弹簧刚度下结构一阶自振频率随${l_0}$的变化曲线如图 3所示，由图 3可知，随着桩基础埋深的增大，自振频率逐渐减小，达到临界埋深后曲线趋于平缓；随着弹簧刚度的增大，梁刚度增大，结构一阶自振频率相应增大。

图 3  桩基础埋深对结构一阶自振频率的影响

Figure 3.  Influence of pile length on natural vibration frequency

下载: 全尺寸图片 幻灯片

3.2.2   上部质量的影响

不同弹簧刚度下结构一阶自振频率随${\sigma _0}$的变化曲线如图 4所示，由图 4可知，随着上部质量的增加，结构一阶自振频率显著减小，且上部质量对地基刚度较大的结构影响更显著。

图 4  上部质量对结构一阶自振频率的影响

Figure 4.  Influence of top mass on natural vibration frequency

下载: 全尺寸图片 幻灯片

3.2.3   转动惯量的影响

结合Jonkman等（2009）给出的海上风机结构机舱质量和轮毂转动惯量参数，定义考虑上部结构转动惯量与未考虑上部结构转动惯量的一阶自振频率比值为R₁，不同弹簧刚度下上部结构转动惯量对结构一阶自振频率的影响如图 5所示。由图 5可知，上部结构转动惯量对结构一阶自振频率的影响较小，可忽略不计。

图 5  上部结构转动惯量对结构一阶自振频率的影响

Figure 5.  Influence of rotary inertia on natural vibration frequency

下载: 全尺寸图片 幻灯片

梁在振动过程中产生横向弯曲变形，使梁发生横向平移，同时发生转动，定义考虑截面转动惯量与未考虑截面转动惯量的一阶自振频率比值为R₂，不同弹簧刚度下截面转动惯量对结构一阶自振频率的影响如图 6所示，由图 6可知，截面转动惯量对结构一阶自振频率的影响很小，可忽略不计。

图 6  截面转动惯量对结构一阶自振频率的影响

Figure 6.  Influence of section inertia on natural vibration frequency

下载: 全尺寸图片 幻灯片

3.2.4   桩-水相互作用的影响

桩-水相互作用是影响海上风机结构设计的重要因素，定义考虑桩-水相互作用与未考虑桩-水相互作用的一阶自振频率比值为R₃，结构一阶自振频率差异如图 7所示，由图 7可知，考虑桩-水相互作用时结构一阶自振频率较小，桩-水相互作用的影响对柔性地基结构更明显。

图 7  桩-水相互作用对结构一阶自振频率的影响

Figure 7.  Influence of hydraulic added on natural vibration frequency

下载: 全尺寸图片 幻灯片

4.   结论

本文基于欧拉-伯努利梁理论和传递矩阵方法，考虑水-桩-土相互作用和塔筒变截面梁特性，建立单桩式海上风机结构系统自振频率解析计算方法，系统的研究了桩基础埋深、上部质量、转动惯量和桩-水相互作用对结构一阶自振频率的影响，主要结论如下：

（1）随着桩基础埋深的增大，风机结构一阶自振频率逐渐减小，到达临界埋深后减小趋势趋于平缓。

（2）随着上部质量的增大，风机结构一阶自振频率显著减小，一方面是因为上部质量的增大增加了风机结构整体质量；另一方面是因为上部质量的增大增加了轴向荷载，降低了梁刚度。

（3）上部结构转动惯量和截面转动惯量对风机结构一阶自振频率的影响较小，在实际工程中可忽略不计。

（4）水-结构相互作用会降低结构一阶自振频率，这种影响在土体刚度较小时更明显，因此在海上风机结构设计中有必要考虑水-结构相互作用的影响。

附录

当$ {k_i} - {m_i}{\omega ^2} \leqslant 0$时：

${t_{11}} = \sin {\delta _i}{l_i}$，${t_{12}} = \cos {\delta _i}{l_i}$，${t_{13}} = \sin {\varepsilon _i}{l_i}$，${t_{14}} = \cos {\varepsilon _i}{l_i}$，${t_{21}} = {\delta _i}\cos {\delta _i}{l_i}$，${t_{22}} = - {\delta _i}\sin{\delta _i}{l_i}$，${t_{23}} = $ ${\varepsilon _i}\cos {\varepsilon _i}{l_i}$，${t_{24}} = - {\varepsilon _i}\sin {\varepsilon _i}{l_i}$，${t_{31}} = - E{I_i}{\delta _i}^2\sin {\delta _i}{l_i}$，${t_{32}} = - E{I_i}{\delta _i}^2\cos {\delta _i}{l_i}$，${t_{33}} = - E{I_i}{\varepsilon _i}^2\sin {\varepsilon _i}{l_i}$，${t_{34}} = $ $ - E{I_i}{\varepsilon _i}^2\cos {\varepsilon _i}{l_i}$，${t_{41}} = - E{I_i}{\delta _i}^3\cos {\delta _i}{l_i}$，${t_{42}} = E{I_i}{\delta _i}^3\sin {\delta _i}{l_i}$，${t_{43}} = - E{I_i}{\varepsilon _i}^3\cos {\varepsilon _i}{l_i}$，${t_{44}} = E{I_i}{\varepsilon _i}^3\sin {\varepsilon _i}{l_i}$，${s_{11}} = {s_{13}} = {s_{22}} = {s_{24}} = {s_{31}} = {s_{33}} = {s_{42}} = {s_{44}} = 0$，${s_{12}} = {s_{14}} = 1$，${s_{21}} = {\delta _{i + 1}}$，${s_{23}} = {\varepsilon _{i + 1}}$，${s_{32}} = - E{I_{i + 1}}\delta _{i + 1}^2$，${s_{34}} = - E{I_{i + 1}}\varepsilon _{i + 1}^2$，${s_{41}} = - E{I_{i + 1}}\delta _{i + 1}^3$，${s_{43}} = - E{I_{i + 1}}\varepsilon _{i + 1}^3$。

当$ 0＜{k}_{i}-{m}_{i}{\omega }^{2}\le {(P+{\rho }_{i}{I}_{i}{\omega }^{2})}^{2}/4E{I}_{i}$时：

${t_{11}} = \sin {\delta _i}{l_i}$，${t_{12}} = \cos {\delta _i}{l_i}$，${t_{13}} = \sinh {\varepsilon _i}{l_i}$，${t_{14}} = \cosh {\varepsilon _i}{l_i}$，${t_{21}} = {\delta _i}\cos {\delta _i}{l_i}$，${t_{22}} = - {\delta _i}\sin{\delta _i}{l_i}$，${t_{23}} = {\varepsilon _i}\cosh {\varepsilon _i}{l_i}$，${t_{24}} = {\varepsilon _i}\sinh {\varepsilon _i}{l_i}$，${t_{31}} = - E{I_i}{\delta _i}^2\sin {\delta _i}{l_i}$，${t_{32}} = - E{I_i}{\delta _i}^2\cos {\delta _i}{l_i}$，${t_{33}} = E{I_i}{\varepsilon _i}^2\sinh {\varepsilon _i}{l_i}$，${t_{34}} = E{I_i}{\varepsilon _i}^2\cosh {\varepsilon _i}{l_i}$，${t_{41}} = - E{I_i}{\delta _i}^3\cos {\delta _i}{l_i}$，${t_{42}} = E{I_i}{\delta _i}^3\sin {\delta _i}{l_i}$，${t_{43}} = E{I_i}{\varepsilon _i}^3\cosh {\varepsilon _i}{l_i}$，${t_{44}} = E{I_i}{\varepsilon _i}^3$ $\sinh {\varepsilon _i}{l_i}$，${s_{11}} = {s_{13}} = {s_{22}} = {s_{24}} = {s_{31}} = {s_{33}} = {s_{42}} = {s_{44}} = 0$，${s_{12}} = {s_{14}} = 1$，${s_{21}} = {\delta _{i + 1}}$，${s_{23}} = {\varepsilon _{i + 1}}$，${s_{32}} = $ $ - E{I_{i + 1}}\delta _{i + 1}^2$，${s_{34}} = E{I_{i + 1}}\varepsilon _{i + 1}^2$，${s_{41}} = - E{I_{i + 1}}\delta _{i + 1}^3$，${s_{43}} = E{I_{i + 1}}\varepsilon _{i + 1}^3$。

当$ {k}_{i}-{m}_{i}{\omega }^{2}{\rm{ }}＞{\rm{ }}{(P+{\rho }_{i}{I}_{i}{\omega }^{2})}^{2}/4E{I}_{i} $时：

${t_{11}} = \cosh ({\gamma _i}{l_i})\cos ({\theta _i}{l_i})$，${t_{21}} = {\gamma _i}{t_{13}} - {\theta _i}{t_{12}}$，${t_{12}} = \cosh ({\gamma _i}{l_i})\sin ({\theta _i}{l_i})$，${t_{22}} = {\gamma _i}{t_{14}} + {\theta _i}{t_{11}}$，${t_{13}} = $ $\sinh ({\gamma _i}{l_i})\cos ({\theta _i}{l_i}){{\rm{ }}}$，${t_{23}} = {\gamma _i}{t_{11}} - {\theta _i}{t_{14}}$，${t_{14}} = \sinh ({\gamma _i}{l_i})\sin ({\theta _i}{l_i})$，${t_{24}} = {\gamma _i}{t_{12}} + {\theta _i}{t_{13}}$，${t_{31}} = {r_i}{t_{11}} - 2{\gamma _i}{\theta _i}E{I_{{\rm{f}}}}{t_{14}}$，${t_{32}} = {r_i}{t_{12}} + 2{\gamma _i}{\theta _i}E{I_{{\rm{f}}}}{t_{13}}$，${t_{33}} = {r_i}{t_{13}} - 2{\gamma _i}{\theta _i}E{I_{{\rm{f}}}}{t_{12}}{{\rm{ }}}$，${t_{34}} = {r_i}{t_{14}} + 2{\gamma _i}{\theta _i}E{I_{{\rm{f}}}}{t_{11}}$，${t_{41}} = ({r_i}{\gamma _i} - 2{\gamma _i}{\theta _i}^2E{I_{{\rm{f}}}}){t_{13}} - $ $({r_i}{\theta _i} + 2{\gamma _i}^2{\theta _i}E{I_{{\rm{f}}}}){t_{12}}$，${t_{42}} = ({r_i}{\gamma _i} - 2{\gamma _i}{\theta _i}^2E{I_{{\rm{f}}}}){t_{14}} + ({r_i}{\theta _i} + 2{\gamma _i}^2{\theta _i}E{I_{{\rm{f}}}}){t_{11}}$，${t_{43}} = ({r_i}{\gamma _i} - 2{\gamma _i}{\theta _i}^2E{I_{{\rm{f}}}}){t_{11}} - ({r_i}{\theta _i} + $ $2{\gamma _i}^2{\theta _i}E{I_{{\rm{f}}}}){t_{14}}$，${t_{44}} = ({r_i}{\gamma _i} - 2{\gamma _i}{\theta _i}^2E{I_{{\rm{f}}}}){t_{12}} + ({r_i}{\theta _i} + 2{\gamma _i}^2{\theta _i}E{I_{{\rm{f}}}}){t_{13}}$，${s_{12}} = {s_{13}} = {s_{14}} = {s_{21}} = {s_{24}} = {s_{32}} = {s_{33}} = {s_{41}}$ $ = {s_{44}} = 0$，${s_{11}} = 1$，${s_{22}} = {\theta _{i + 1}}$，${s_{23}} = {\gamma _{i + 1}}$，${s_{31}} = {r_{i + 1}}$，${s_{34}} = 2{\gamma _{i + 1}}{\theta _{i + 1}}E{I_{{\rm{f}}}}$，${s_{42}} = {r_{i + 1}}{\theta _{i + 1}} + 2\gamma _{i + 1}^2{\theta _{i + 1}}E{I_{{\rm{f}}}}$，${s_{43}} = {r_{i + 1}}{\gamma _{i + 1}} - 2{\gamma _{i + 1}}\theta _{i + 1}^2E{I_{{\rm{f}}}}$，${r_i} = P + {\rho _i}{I_i}{\omega ^2} + ({\gamma _i}^2 - {\theta _i}^2)E{I_{{\rm{f}}}}$。