考虑水‑桩‑土相互作用的海上风机单桩基础横向自振频率计算

余云燕，陶婧衍，孔嘉乐，李永鹏; YU Yunyan; TAO Jingyan; KONG Jiale; LI Yongpeng

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

考虑水‑桩‑土相互作用的海上风机单桩基础横向自振频率计算 PDF

- ORCID：
余云燕 ¹
✉
- ORCID：
陶婧衍 ¹
- ORCID：
孔嘉乐 ¹
- ORCID：
李永鹏 ²

1. 兰州交通大学土木工程学院，甘肃兰州 730070； 2. 甘南藏族自治州交通运输局，甘肃甘南 747000

中图分类号： TK83； TU473.1

最近更新：2025-03-06

DOI：10.16385/j.cnki.issn.1004-4523.2025.02.012

摘要

采用Euler‑Beroulli梁理论，提出一种考虑水‑桩‑土相互作用的连续变截面海上风机单桩基础的结构简化模型，运用微分变换法求解简化模型横向振动控制方程。基于塔筒直径、连接段高度与水体附加质量、叶轮‑机舱组件质量和三弹簧刚度对横向自振频率展开研究。研究表明，变截面塔筒底部直径对自振频率的影响大于顶部直径对自振频率的影响；水深较深时，水体附加质量对远海风机结构自振频率的影响不能忽略；风机横向自振频率随叶轮‑机舱组件质量的增大而减小；弹簧刚度对土体模量敏感性大小为：水平弹簧>耦合弹簧>旋转弹簧；自振频率对弹簧刚度敏感性大小为：耦合弹簧>水平弹簧>旋转弹簧；土体模量发生变化时，主要通过水平弹簧和耦合弹簧影响自振频率。

关键词

海上风机; 弹簧地基; 自振频率; 大直径单桩; 微分变换法

海上风电由于其风速高、风向稳定、占地面积小以及噪声低等优点，成为优质的可再生能源。截至目前，中国海上风机累计装机容量约40 GW^［

1］，海上风电容量居全球第一。然而，中国海上风机开发技术总体滞后于欧洲，缺乏适用于中国海洋环境下的海上风机结构设计理论，因此，海上风机结构设计仍面临严峻挑战。

截至2020年，欧洲已建成的海上风机超过80%采用单桩基础，中国已建或在建的海上风电项目，超过70%的风机也采用单桩基础^［

2］。单桩基础因维修便利、施工简单、受力明确等原因，已成为海上风电场建设中应用最广泛的基础形式。

风机在运营过程中会受到涡轮机转动产生的1P频率和叶片扫掠塔筒产生的2P/3P频率影响。一旦发生频率偏移，系统自振频率接近1P或2P/3P频率，会导致结构动态响应放大，产生过大位移，将对结构稳定性产生不利影响。为了避免系统共振，通常采用“软‑刚”模式^［

3］使自振频率介于1P和2P/3P频率带之间。因此，设计过程中需要对大直径单桩基础风机结构的横向自振频率进行精确计算。

在风机前期设计中，考虑桩‑土相互作用可以准确求解结构横向自振频率。目前，单桩式风机结构的桩‑土相互作用可以采用表观固定法、耦合弹簧法和分布弹簧法模拟。其中，席义博等^［

4］通过有限元数值模拟，采用表观固定法建立考虑基础刚度弱化的等截面风机模型，并进行不同基础刚度条件下风机的动力响应分析。ZAAIJER^{［参考文献 5

百度学术}5］推导了底端刚性固定模型、非耦合弹簧模型和耦合弹簧模型的结构动力特性，分析表明把桩‑土相互作用等效为泥线处的耦合弹簧模型结果更精确。BHATTACHARYA等^{［参考文献 6

百度学术}6］、ADHIKARI等^{［参考文献 7

百度学术}7］和SHADLOU等^{［参考文献 8

百度学术}8］将上部结构考虑为集中质量和等截面杆系统，桩‑土相互作用简化为泥面处施加的水平和旋转弹簧，建立柔性地基条件下计算风机结构自振频率的解析公式，并与模型试验和有限元计算结果进行对比。ARANY等^{［参考文献 9‑11}9‑11］在BHATTACHARYA的理论基础上将地基作用简化为水平、旋转和耦合弹簧模型，推导了海上风机自振频率的闭合解。杨春宝等^{［参考文献 12

百度学术}12］采用多段均匀梁组合方式将风机塔筒等效为变截面梁，并将桩‑土相互作用简化为线性分布弹簧，基于传递矩阵和有限元模态分析方法分别计算有阻尼系统的风机自振频率。

研究表明水对结构动力响应有一定影响^［

13］，部分学者对水‑结构相互作用开展了一系列分析。其中，赵密等^{［参考文献 14‑15}14‑15］发现在海上桩基设计中考虑水‑桩‑土相互作用是非常必要的，并基于杨春宝的理论，建立附加水质量的解析计算方法。许成顺等^{［参考文献 16

百度学术}16］采用OpenSees有限元平台，探讨水体附加质量对钢管桩自振频率的影响，认为水体附加质量对钢管桩自振频率影响显著。任文渊^{［参考文献 17

百度学术}17］采用ADINA进行数值模拟，建立风机、基础、水、海床的三维计算模型，分析了风电基础结构的静力和考虑流固耦合的动力响应，结果表明水层显著降低结构的高阶频率。田树刚等^{［参考文献 18

百度学术}18］采用有限元方法，考虑流固耦合效应，对近海风机支撑体系的自振频率和振型进行计算分析。

综上所述，以往研究通常将风机塔桶截面简化为等截面梁，忽略风机塔筒变截面特性，致使频率计算结果不精确；或将风机塔桶结构划分为多段均匀梁，此方法虽然能计算风机各阶自振频率，但求解过程中系数冗杂。针对上述问题，本文将风机塔筒简化为连续变截面梁，引入微分变换法（differential transform method，DTM）将非线性控制方程变换成代数方程进行求解^［

19］，解决了四阶变系数控制方程难以求解的计算难点。通过DTM建立质量块‑连续变截面梁‑等截面梁组合模型，在考虑水‑桩‑土相互作用的基础上，提出一种大直径单桩海上风机系统横向自振频率的求解方法，并进一步研究塔筒直径、连接段高度、叶轮‑机舱组件质量和三弹簧刚度等参数对风机结构自振频率的影响规律。

1 方程的建立与求解

1.1 地基模型简化

由桩‑土相互作用提供给桩顶的刚度简化为水平 $K_{L}$ ，耦合 $K_{L R}$ 和旋转 $K_{R}$ 三个独立弹簧，可用如下刚度矩阵表示：

(\begin{array}{l} F_{x} \\ M \end{array}) = (\begin{matrix} K_{L} & K_{L R} \\ K_{L R} & K_{R} \end{matrix}) (\begin{array}{l} u \\ θ \end{array})

（1）

式中， $K_{L}$ 、 $K_{L R}$ 和 $K_{R}$ 分别为柱顶的水平、耦合和旋转弹簧刚度； $F_{x}$ 和 $M$ 分别为柱顶的水平力和弯矩； $u$ 和 $θ$ 分别为柱顶的位移和转角。对于均质土体、线性土体以及非线性土体的三弹簧刚度计算式均不相同，并且弹簧刚度大小与桩径、桩长、桩的模量和土体泊松比、土体模量等参数有关。文献［

8］分别给出了柔性桩和刚性桩的弹簧刚度公式，结果如表1和2所示。

表1 柔性桩的弹簧刚度公式^{［参考文献 8

百度学术}8］

Tab.1 Spring stiffness formula for flexible pile^{［参考文献 8

百度学术}8］

土体类型		$K_{L}$	$K_{L R}$	$K_{R}$
均质		$2.9 {(\frac{E_{p}}{E_{s o}})}^{0.186} f_{(ν s)} E_{s o} D_{p}$	$- 1.2 {(\frac{E_{p}}{E_{s o}})}^{0.5} f_{(ν s)} E_{s o} D_{p}^{2}$	$1.5 {(\frac{E_{p}}{E_{s o}})}^{0.73} f_{(ν s)} E_{s o} D_{p}^{3}$
非均质	非线性	$2.03 {(\frac{E_{p}}{E_{s o}})}^{0.27} f_{(ν s)} E_{s o} D_{p}$	$- 1.17 {(\frac{E_{p}}{E_{s o}})}^{0.52} f_{(ν s)} E_{s o} D_{p}^{2}$	$1.42 {(\frac{E_{p}}{E_{s o}})}^{0.76} f_{(ν s)} E_{s o} D_{p}^{3}$
非均质	线性	$1.58 {(\frac{E_{p}}{E_{s o}})}^{0.34} f_{(ν s)} E_{s o} D_{p}$	$- 1.07 {(\frac{E_{p}}{E_{s o}})}^{0.567} f_{(ν s)} E_{s o} D_{p}^{2}$	$1.38 {(\frac{E_{p}}{E_{s o}})}^{0.78} f_{(ν s)} E_{s o} D_{p}^{3}$

注： $E_{p}$ 为桩弹性模量； $E_{s o}$ 为土体模量； $f_{(ν s)} = 1 + 0.6 |ν_{s} - 0.25|$ ， $ν_{s}$ 为土体泊松比； $D_{p}$ 为桩径。

表2 刚性桩的弹簧刚度公式^{［参考文献 8

百度学术}8］

Tab.2 Spring stiffness formula for rigid pile^{［参考文献 8

百度学术}8］

土体类型		$K_{L}$	$K_{L R}$	$K_{R}$
均质		$6.4 {(\frac{L_{p}}{D_{p}})}^{0.162} f_{(ν s)} E_{s o} D_{p}$	$- 7.1 {(\frac{L_{p}}{D_{p}})}^{1.56} f_{(ν s)} E_{s o} D_{p}^{2}$	$13.2 {(\frac{L_{p}}{D_{p}})}^{2.5} f_{(ν s)} E_{s o} D_{p}^{3}$
非均质	非线性	$5.33 {(\frac{L_{p}}{D_{p}})}^{1.07} f_{(ν s)} E_{s o} D_{p}$	$- 7.2 {(\frac{L_{p}}{D_{p}})}^{2} f_{(ν s)} E_{s o} D_{p}^{2}$	$13 {(\frac{L_{p}}{D_{p}})}^{3} f_{(ν s)} E_{s o} D_{p}^{3}$
非均质	线性	$4.7 {(\frac{L_{p}}{D_{p}})}^{1.53} f_{(ν s)} E_{s o} D_{p}$	$- 7.1 {(\frac{L_{p}}{D_{p}})}^{2.5} f_{(ν s)} E_{s o} D_{p}^{2}$	$12.7 {(\frac{L_{p}}{D_{p}})}^{3.45} f_{(ν s)} E_{s o} D_{p}^{3}$

注： $L_{p}$ 为桩长。

1.2 塔筒模型简化与控制方程建立

风机单桩基础结构简化模型如图1所示，考虑叶轮‑机舱组件质量产生的轴向力作用，塔筒简化为连续变截面梁，连接段考虑为等截面梁，水体附加质量添加至连接段上，基础简化为三弹簧模型。基于Bernoulli‑Euler梁理论建立塔筒横向振动方程^［

20］：

\begin{array}{l} \frac{\partial^{2}}{\partial x^{2}} [E I (x) \frac{\partial^{2} y_{t} (x, t)}{\partial x^{2}}] + P \frac{\partial^{2} y_{t} (x, t)}{\partial x^{2}} + \\ ρ A (x) \frac{\partial^{2} y_{t} (x, t)}{\partial t^{2}} = 0, L_{c} < x < L_{c} + L_{t} \end{array}

（2）

式中， $y_{t} (x, t)$ 为梁（即塔筒）横向位移； $x$ 为沿梁轴线坐标； $t$ 为时间； $I (x)$ 和 $A (x)$ 分别为 $x$ 截面处塔筒的惯性矩和面积； $L_{t}$ 为塔筒高度； $L_{c}$ 为连接段高度； $P$ 为修正后的风机轴向力， $P = (M_{R} + C_{m} m_{t}) g$ ^［

9］，其中

M_{R}

为风机转子和轮毂集中质量，

C_{m}

为质量矫正因子，取值为0.24，

m_{t}

为塔筒质量；

E

和

ρ

分别为风机结构的弹性模量和材料密度。

图1 风机结构模型

Fig.1 Wind turbine structure model

建立连接段横向振动方程：

\begin{array}{l} \frac{\partial^{2}}{\partial x^{2}} [E I_{c} \frac{\partial^{2} y_{c} (x, t)}{\partial x^{2}}] + P \frac{\partial^{2} y_{c} (x, t)}{\partial x^{2}} + \\ ρ A_{c} \frac{\partial^{2} y_{c} (x, t)}{\partial t^{2}} + m_{a} \frac{\partial^{2} y_{c} (x, t)}{\partial t^{2}} = 0,0 < x < L_{c} \end{array}

（3）

式中， $y_{c} (x, t)$ 为梁（即连接段）横向位移； $I_{c}$ 和 $A_{c}$ 分别为连接段的截面惯性矩和截面面积； $m_{a}$ 为水引起的附加质量。

图1中 $t_{t}$ 为塔筒壁厚； $D_{t}$ 为塔筒顶部截面直径； $D_{b}$ 为塔筒底部截面直径，引入截面变化系数 $γ_{D} = \frac{D_{t}}{D_{b}} - 1$ ，则塔筒直径、截面惯性矩和截面面积可分别表示为：

\{\begin{array}{l} D (x) = D_{b} (1 + γ_{D} \frac{x}{L_{t}}) \\ \begin{array}{l} I (x) = \frac{1}{8} π D^{3} (x) t_{t} - \frac{3}{8} π D^{3} (x) t_{t}^{2} + \\ \frac{1}{2} π D (x) t_{t}^{3} - \frac{1}{4} π t_{t}^{4} \end{array} \\ A (x) = π [D (x) t_{t} - t_{t}^{2}] \end{array}

（4）

本文将结构‑水相互作用产生的动水压力视为部分水体质量与结构加速度的乘积，部分水体质量称附加质量 $m_{a}$ ，表达式^［

14］如下：

\begin{array}{l} m_{a} = m_{0} (0.4327 e^{- 5.844 \frac{D_{c}}{2 L_{c}}} + 0.5369 e^{- 0.0781 \frac{D_{c}}{2 L_{c}}}), \\ m_{0} = ρ_{w} π D_{c}^{2} (0.6 e^{- 0.93 \frac{D_{c}}{L_{c}}} + 0.403 e^{- 0.156 \frac{D_{c}}{L_{c}}}) \end{array}

（5）

式中， $m_{0}$ 为结构刚性运动引起的附加质量； $D_{c}$ 为连接段直径； $ρ_{w}$ 为水体密度。

根据力平衡和位移协调，连接段底部边界条件表示为：

\{\begin{array}{l} E I_{c} \frac{\partial^{2} y_{c} (0, t)}{\partial x^{2}} - K_{R} \frac{\partial y_{c} (0, t)}{\partial x} + K_{L R} y_{c} (0, t) = 0 \\ E I_{c} \frac{\partial^{3} y_{c} (0, t)}{\partial x^{3}} + K_{L R} \frac{\partial y_{c} (0, t)}{\partial x} - K_{L} y_{c} (0, t) = 0 \end{array}

（6）

塔筒顶部边界条件表示为：

\{\begin{array}{l} E I_{t} \frac{\partial^{2} y_{t} (L_{t}, t)}{\partial x^{2}} = 0 \\ E I_{t} \frac{\partial^{3} y_{t} (L_{t}, t)}{\partial x^{3}} + P \frac{\partial y_{t} (L_{t}, t)}{\partial x} - M_{R} \frac{\partial^{2} y_{t} (L_{t}, t)}{\partial t^{2}} = 0 \end{array}

（7）

式中， $I_{t}$ 为塔筒顶部截面惯性矩。

1.3 微分变换法

引入DTM将式（2）和（3）这类变系数微分方程转换为代数递推方程，将式（6）和（7）等边界条件转化为代数方程组，通过求解代数方程的解来获得微分方程的解，解决了数学求解问题。微分变换及微分逆变换分别表示为^［

21］：

\bar{F} (k) = \frac{1}{k!} {(\frac{d^{k} f (r)}{d r^{k}})}_{r = r_{0}}

（8）

f (r) = \sum_{k = 0}^{\infty} \bar{F} (k) (r - r_{0})^{k}

（9）

式中， $k$ 表示函数 $f (r)$ 对自变量 $r$ 的求导阶数； $r_{0}$ 取为任何正整数。

表3列出了本文计算中用到的DTM变换定理。

表3 DTM基本变换定理

Tab.3 Basic transformation theorem of DTM

原函数	变换函数
$f (r) = y (r) \pm z (r)$	$F (k) = Y (k) \pm Z (k)$
$f (r) = c y (r)$	$F (k) = c Y (k)$
$f (r) = \frac{d^{n} y (r)}{d x^{n}}$	$F (k) = (k + 1) (k + 2) \dots Y (k + n)$
$f (r) = y (r) z (r)$	$F (k) = \sum_{l = 0}^{k} G (l) H (k - l)$
$f (r) = r^{n}$	$F (k) = δ (k - n) = \{\begin{matrix} 1, k = n \\ 0, k \neq n \end{matrix}$

注：原函数用小写字母表示，变换后的函数用对应的大写字母表示； $c$ 为非零常数； $δ$ 为狄克拉函数。

1.4 频率方程建立

当梁作自由振动时，其解有如下形式：

\{\begin{array}{l} y_{t} (x, t) = Y_{t} (x) e^{i ω t} \\ y_{c} (x, t) = Y_{c} (x) e^{i ω t} \end{array}

（10）

式中， $Y (x)$ 为梁的横向振动函数； $ω$ 为梁横向自由振动圆频率； $i = \sqrt[]{- 1}$ 。

引入无量纲系数 $ξ = \frac{x}{L_{t}}$ 和 $ε = \frac{x}{L_{c}}$ ，并代入式（2）、（6）和（7）中，无量纲处理后有：

\begin{array}{l} I (ξ) {Y_{t}^{″}}^{″} (ξ) + 2 I^{'} (ξ) Y_{t}^{‴} (ξ) + I^{″} (ξ) Y_{t}^{″} (ξ) + \\ P_{t} Y_{t}^{″} (ξ) - Ω^{2} A (ξ) Y_{t} (ξ) = 0 \end{array}

（11）

{Y_{c}^{″}}^{″} (ε) + P_{c} Y_{c}^{″} (ε) - Φ^{2} Y_{c} (ε) - M^{2} Y_{c} (ε) = 0

（12）

\{\begin{array}{l} Y_{c}^{″} (0) - η_{R} Y_{c}^{'} (0) + η_{L R} Y_{c} (0) = 0 \\ Y_{c}^{‴} (0) + η_{L R} Y_{c}^{'} (0) - η_{L} Y_{c} (0) = 0 \end{array}

（13）

\{\begin{array}{l} Y_{t}^{″} (1) = 0 \\ Y_{t}^{‴} (1) + \frac{P_{t}}{I_{t}} Y_{t}^{'} (1) + α Ω^{2} Y_{t} (1) = 0 \end{array}

（14）

连接段顶部和塔筒底部边界为：

\{\begin{array}{l} Y_{t} (0) = A_{1} Y_{c} (1) \\ Y_{t}^{'} (0) = A_{2} Y_{c}^{'} (1) \\ Y_{t}^{″} (0) = A_{3} Y_{c}^{″} (1) \\ Y_{t}^{‴} (0) = A_{4} Y_{c}^{‴} (1) \end{array}

（15）

式（11）~（15）中的无量纲参数为：

$P_{t} = \frac{P L_{t}^{2}}{E}$ ； $Ω = ω \sqrt[]{\frac{m_{t} L_{t}^{4}}{E I_{t}}}$ ； $P_{c} = \frac{P L_{c}^{2}}{E I}$ ； $Φ = ω \sqrt[]{\frac{ρ A_{c} L_{c}^{4}}{E I_{c}}}$ ； $M = ω \sqrt[]{\frac{m_{a} L_{c}^{4}}{E I_{c}}}$ ； $η_{R} = \frac{K_{R} L_{c}}{E I_{c}}$ ； $η_{L R} = \frac{K_{L R} L_{c}^{2}}{E I_{c}}$ ； $η_{L} = \frac{K_{L} L_{c}^{3}}{E I_{c}}$ ； $α = \frac{M_{R}}{m_{t} L_{t}}$ ； $A_{1} = \frac{L_{c}}{L_{t}}$ ； $A_{2} = \frac{I_{c}}{I_{t}}$ ； $A_{3} = \frac{I_{c} L_{t}}{I_{t} L_{c}}$ ； $A_{4} = \frac{I_{c} L_{t}^{2}}{I_{t} L_{c}^{2}}$ ； $A (ξ) = π [D (ξ) t_{t} - t_{t}^{2}]$ ； $I (ξ) = \frac{1}{8} π D^{3} (ξ) t_{t} - \frac{3}{8} π \cdot D^{2} (ξ) t_{t}^{2} + \frac{1}{2} π D (ξ) t_{t}^{3} - \frac{1}{4} π t_{t}^{4} 。$ 根据DTM变换法则，控制方程式（11）和（12）可转换为：

\begin{array}{l} \sum_{l = 0}^{k} {Q (l) (k - l + 4) (k - l + 3) (k - l + 2) (k - l + 1) W (k - l + 4)} + 2 {(l + 1) Q (l + 1) \cdot \\ (k - l + 3) (k - l + 2) (k - l + 1) W (k - l + 3)} + \sum_{l = 0}^{k} {(l + 2) (l + 1) Q (l + 2) \cdot \\ (k - l + 2) (k - l + 1) W (k - l + 2)} + (k + 2) (k + 1) W (k + 2) + Ω^{2} \sum_{l = 0}^{k} {P (l) W (k - l)} = 0 \end{array}

（16）

\begin{array}{l} (k + 4) (k + 3) (k + 2) (k + 1) D (k + 4) + \\ (k + 2) (k + 1) D (k + 2) - \\ Φ^{2} D (k) - M^{2} D (k) = 0 \end{array}

（17）

连接段底部边界条件表示为：

\{\begin{array}{l} D (2) = \frac{1}{2} [η_{R} D (1) - η_{L R} D (0)] \\ D (3) = \frac{1}{6} [- η_{L R} D (1) + η_{L} D (0)] \end{array}

（18）

塔筒底部与连接段顶部边界表示为：

\{\begin{array}{l} W (0) = A_{1} \sum_{k = 0}^{N} D (k) \\ W (1) = A_{2} \sum_{k = 0}^{N} (k + 1) D (k + 1) \\ W (2) = \frac{1}{2} A_{3} \sum_{k = 0}^{N} (k + 2) (k + 1) D (k + 2) \\ W (3) = \frac{1}{6} A_{4} \sum_{k = 0}^{N} (k + 3) (k + 2) (k + 1) D (k + 3) \end{array}

（19）

式中， $N$ 为泰勒级数的展开项数。

塔筒顶部边界条件表示为：

\{\begin{array}{l} \sum_{k = 0}^{N} (k + 2) (k + 1) W (k + 2) = 0 \\ \begin{array}{l} \sum_{k = 0}^{N} (k + 3) (k + 2) (k + 1) W (k + 3) = \\ - v (k + 1) W (k + 1) - α Ω^{2} W (k) \end{array} \end{array}

（20）

式中， $Q (k)$ 、 $W (k)$ 、 $P (k)$ 和 $D (k)$ 分别为 $I (ξ)$ 、 $Y_{t} (ξ)$ 、 $A (ξ)$ 和 $Y_{c} (ε)$ 的微分变换式。

根据边界条件式（18），令：

D (0) = a

，

D (1) = b

（21）

将式（18）和（21）代入连接段位移函数递推式（17）中，所得结果代入连接段与塔筒的边界条件式（19）中，得到塔筒位移函数初值并代入塔筒位移函数递推式（16），最后将结果代入式（20）中，得到频率特征方程为：

(\begin{matrix} T_{11}^{N} (ω) & T_{12}^{N} (ω) \\ T_{21}^{N} (ω) & T_{22}^{N} (ω) \end{matrix}) (\begin{matrix} a \\ b \end{matrix}) = (\begin{matrix} 0 \\ 0 \end{matrix})

（22）

式中， $T_{11}^{N}$ 、 $T_{12}^{N}$ 、 $T_{21}^{N}$ 和 $T_{22}^{N}$ 为迭代 $N$ 次求出的含有自振频率 $ω$ 的多项式，要使上式有非零解，则系数矩阵行列式为零，即

|\begin{matrix} T_{11}^{N} (ω) & T_{12}^{N} (ω) \\ T_{21}^{N} (ω) & T_{22}^{N} (ω) \end{matrix}| = 0

（23）

利用MATLAB求解式（23），得到梁自振频率。

2 算例分析

2.1 有效性验证

以文献［

11］和［22］中的风机规格为例，具体参数如表4所示，将参数代入上述模型中，求解风机系统一阶横向自振频率。在DTM中，展开项数N的取值决定自振频率的收敛性和精度，N越大计算精度越高，但计算耗时越长。为了在计算效率和精度之间取得平衡，绘制不同展开项数N下的频率变化曲线，结果如图2所示。由图2可知，随着展开项数增加，自振频率在初始阶段突然下降，随后缓慢上升并逐渐趋于稳定；当

N = 14

时，Lely A2风机和Irene Vorrink风机的计算误差

|ω^{N + 1} - ω^{N}|

均为0.5‰，小于容许误差，即计算精度满足要求；当

N \geq 14

时，频率变化率基本为0，即自振频率收敛。因此，本文取展开项数14次进行频率计算。

表4 风机结构输入参数^{［参考文献 11

百度学术}11，22］

Tab.4 Wind turbine structure input parameters^{［参考文献 11

百度学术}11，22］

规格参数/单位	Lely A2	Irene Vorrink
$M_{R}$ /t	32	35.7
$m_{t}$ /t	31.44	37
$L_{t}$ /m	41.5	51
$D_{t}$ /m	1.9	1.7
$D_{b}$ /m	3.2	3.5
$L_{c}$ /m	4.6	3.8
$D_{c}$ /m	3.7	3.5
$t_{t}$ /m	12	14
$E$ /GPa	210	210
$ρ$ /（kg·m^-3）	7850	7850
$K_{L}$ /（GN·m^-1）	0.52	0.58
$K_{L R}$ /GN	-2.74	-3.25
$K_{R}$ /（GN·m·rad^-1）	23.63	29.67

图2 频率收敛性曲线

Fig.2 Frequency convergence curve

为验证算法正确性和地基模型简化方法正确性，将本文计算频率和ARANY计算结果^［

11］分别与实测频率对比并计算相对误差，结果见表5。由表5可知，本文三弹簧DTM算法自振频率最大误差为0.79%，与实测频率基本吻合，而文献［11］采用的传递矩阵法计算结果误差较大，说明在求解风机横向自振频率问题上，采用DTM方法有效且精度高。对比表5中三弹簧模型频率与固接模型频率可知，三弹簧模型误差小且精度满足要求，故将地基简化为三弹簧模型（柔性地基）的方式能有效模拟海上单桩风机系统的桩‑土相互作用。

表5 一阶横向自振频率计算结果与误差

Tab.5 Calculation results and errors of the first-order lateral natural frequency

风机

实测

频率/Hz

本文三弹簧结果

ARANY三弹簧结果^{[参考文献 11

百度学术}11]

本文固接结果

ARANY固接结果^{[参考文献 11

百度学术}11]

频率/Hz

误差/%

频率/Hz

误差/%

频率/Hz

误差/%

频率/Hz

误差/%

Lely A2

0.634

0.639

0.79

0.643

1.41

0.708

11.67

0.713

12.46

Irene Vorrink

0.546

0.549

0.55

0.552

1.10

0.584

6.96

0.583

6.78

2.2 参数分析

在海上风机结构的设计和研究过程中，进行自振频率敏感性分析有助于指导设计优化、评估系统稳定性，这对于确保风机安全运行至关重要。本文基于Lely A2风机尺寸，定量分析了塔筒直径、连接段高度、叶轮‑机舱组件质量和三弹簧刚度等因素变化对风机系统一阶横向自振频率的影响规律。

2.2.1 塔筒直径变化的影响

为定量分析塔筒直径变化对自振频率影响，定义塔筒顶部直径与底部直径比值为锥率 $β_{D}$ ，绘制塔筒直径变化对自振频率的影响曲线，如图3所示。其中，0.59为Lely A2风机锥率， $β_{D}$ 变化产生三种工况：工况一，塔筒顶部直径 $D_{t}$ 保持1.9 m不变，塔筒底部直径 $D_{b}$ 从3.2 m增加至4.4 m，即 $β_{D} < 0.59$ ；工况二，塔筒顶部直径 $D_{t}$ 和底部直径 $D_{b}$ 分别以1.9 m和3.2 m为基准同时变化，即 $β_{D} = 0.59$ ；工况三，塔筒顶部直径 $D_{t}$ 从1.9 m增加至3.1 m，塔筒底部直径 $D_{b}$ 保持3.2 m不变，即 $β_{D} > 0.59$ 。

图3 塔筒直径变化对风机自振频率的影响曲线

Fig.3 Impact curve of tower diameter variation on wind turbine natural frequency

由图3可知，工况一中风机自振频率从0.639 Hz增加到0.955 Hz；工况二中风机自振频率从0.639 Hz增加到0.862 Hz；工况三中风机自振频率从0.639 Hz增加到0.655 Hz。对比三种工况可以看出，当其余参数保持不变，随着底部直径 $D_{b}$ 增大，自振频率显著增加；随着顶部直径 $D_{t}$ 增大，自振频率变化幅度很小。在塔筒直径两个控制截面中，自振频率对底部直径十分敏感，而顶部直径变化对自振频率影响很小。因此，在风机设计时，只要确保塔筒顶部直径满足一定要求即可，关键在于选取塔筒底部直径。

2.2.2 连接段高度的影响

随着能源紧缺日益严峻，海上风机必将向深海和大型化方向发展，要求桩基础和上部结构能在更深的水体环境中稳定运行^［

23］，这直接影响到连接段高度和水体附加质量的大小。当考虑水体附加质量时，随连接段高度增加，水体附加质量增大，而连接段高度增加和水体附加质量增大均导致自振频率降低，二者之间存在叠加效应。为厘清连接段高度和有无水体附加质量对自振频率的影响，现绘制未考虑水体附加质量和考虑水体附加质量两种情况下连接段高度变化对自振频率的影响曲线，如图4所示。

图4 连接段高度与水体附加质量对风机自振频率的影响曲线

Fig.4 Impact curve of connection section height and water‑added mass on wind turbine natural frequency

由图4可知，以连接段高4.6 m为基准，逐渐增加至44.6 m。未考虑水体附加质量时，自振频率减小了42.2%；考虑水体附加质量的情况下，自振频率减小了49.5%。当连接段高度较小时，水体附加质量对自振频率影响不显著；而连接段高度超过20 m后，两条曲线呈现分离趋势，水体附加质量对系统自振频率的影响逐渐显现。因此，在风机设计阶段，近海风机在水深较浅时可以忽略水体附加质量影响，但远海风机需考虑水体附加质量对自振频率偏移的影响。

2.2.3 叶轮‑机舱组件质量的影响

考虑不同功率涡轮机质量，即叶轮‑机舱组件质量大小对自振频率的影响，结果如图5所示。由图5可知，随着叶轮‑机舱组件质量增大，自振频率减小。当叶轮‑机舱组件质量增大至50 t后，曲线斜率也逐渐减小，说明叶轮‑机舱组件质量对风机自振频率的影响程度逐渐降低。因此，在设计过程中为保证风机安全性，通过合理增重可以降低叶轮‑机舱组件质量对自振频率的影响。

图5 叶轮‑机舱组件质量对自振频率的影响曲线

Fig.5 Impact curve of impeller-cabin component mass on natural frequency

2.2.4 三弹簧刚度的影响

本文采用三弹簧地基模型，桩‑土相互作用通过改变三种弹簧刚度大小影响风机自振频率，为厘清三弹簧刚度变化对自振频率的影响机理，进行如下理论分析。现已知Lely A2风机场地最上层为软黏土，其下为致密砂层，由文献［

24］可知土体模量在5~50 MPa之间。根据表1公式，计算出不同土体模量下三种弹簧刚度大小，由此绘制土体模量变化对弹簧刚度的影响曲线，如图6所示。

图6 土体模量对弹簧刚度的影响曲线

Fig.6 Impact curve of soil modulus on spring stiffness

由图6可知，随土体模量从5 MPa增加至50 MPa，水平弹簧刚度 $K_{L}$ 由0.27 GN/m增加到1.67 GN/m；耦合弹簧刚度 $K_{L R}$ 由-1.17 GN减小到-3.71 GN；旋转弹簧刚度 $K_{R}$ 从21.24 GN·m/rad增加到39.56 GN·m/rad。分析表明，土体模量对水平弹簧刚度影响最大，耦合弹簧刚度次之，旋转弹簧刚度影响最小。

图7描绘了风机横向自振频率随三种弹簧刚度变化规律。当其他参数保持不变，随着水平弹簧刚度 $K_{L}$ 增加，自振频率迅速增加，当 $K_{L}$ 增加至0.95 GN/m后，自振频率增长幅度基本保持不变。说明 $K_{L}$ 越小对自振频率影响越大，而 $K_{L}$ 增大至一定程度后，对自振频率影响不大。当其他参数保持不变，减小耦合弹簧刚度 $K_{L R}$ 时，自振频率先缓慢减小，当 $K_{L R}$ 达到-2.63 GN后，自振频率急剧减小，共减小了45.5%。说明 $K_{L R}$ 较大时对自振频率影响较小，随 $K_{L R}$ 减小频率变化幅度急剧增大，耦合弹簧刚度 $K_{L R}$ 对自振频率有衰减作用。随旋转弹簧刚度 $K_{R}$ 增加，频率从0.62 Hz增加至0.68 Hz，增大了9.7%，且 $K_{R}$ 增大至32 GN·m/rad后，自振频率增大趋势变缓慢；说明旋转弹簧刚度 $K_{R}$ 较小时对自振频率影响较大；相反， $K_{R}$ 较大时对自振频率影响较小。上述分析表明：耦合弹簧刚度对自振频率偏移影响最显著，其次是水平弹簧刚度，最后为旋转弹簧刚度。

图7 弹簧刚度对自振频率的影响曲线

Fig.7 Impact curve of spring stiffness on natural frequency

综上，三种弹簧刚度对土体模量敏感性顺序为：水平弹簧>耦合弹簧>旋转弹簧；风机系统横向自振频率对三种弹簧刚度敏感性顺序为：耦合弹簧>水平弹簧>旋转弹簧。综合分析可知，采用三弹簧地基模拟桩‑土相互作用时，土体模量发生变化，主要通过三种弹簧中的水平弹簧和耦合弹簧影响风机系统横向自振频率。

3 结论

本文基于Euler‑Bernoulli梁振动理论，采用微分变化法，充分考虑塔筒变截面特性和水‑桩‑土间的相互作用，提出单桩式海上风机的横向自振频率求解计算方法，进一步分析了塔筒直径、连接段高度和水体附加质量、叶轮‑机舱组件质量对单桩式海上风机横向自振频率的影响，并厘清地基模量通过三弹簧刚度对自振频率的影响机理，得出以下结论：

（1）本文算法频率计算结果与实测频率吻合良好，说明采用DTM求解四阶变系数微分方程原理简单，结果精度高，适用性强，为解决此类问题提供了一种有效方法。

（2）底部直径对横向自振频率影响程度更高；风机横向自振频率对连接段高度更敏感，当连接段高度超过20 m后，水体附加质量对风机自振频率影响不能忽略；随叶轮‑机舱组件质量增大，风机自振频率逐渐减小，当叶轮‑机舱组件质量增加至50 t后，其敏感性逐渐下降。

（3）随土体模量增加，水平和旋转弹簧刚度增加，耦合弹簧刚度减小，其中土体模量对水平弹簧刚度影响程度最大，耦合弹簧刚度次之，旋转弹簧刚度最小；风机自振频率随水平弹簧刚度和旋转弹簧刚度增加而增大，随耦合弹簧刚度减小而减小，其中耦合弹簧刚度对自振频率偏移影响最显著，水平弹簧刚度次之，旋转弹簧刚度较小。综上，土体模量变化主要通过水平弹簧和耦合弹簧影响自振频率。

参考文献

王立忠，洪义，高洋洋，等. 近海风电结构台风环境动力灾变与控制［J］. 力学学报， 2023， 55（3）： 567-587. [百度学术]

WANG Lizhong， HONG Yi， GAO Yangyang， et al. Dynamic catastrophe and control of offshore wind power structures in typhoon environment［J］. Chinese Journal of Theoretical and Applied Mechanics， 2023， 55（3）： 567-587. [百度学术]

赖踊卿. 软黏土地基海上风机大直径单桩水平受荷特性与分析模型［D］. 杭州：浙江大学， 2021. [百度学术]

LAI Yongqing. Modelling of lateral behaviour of large-diameter monopiles supporting offshore wind turbines in soft clay［D］. Hangzhou： Zhejiang University， 2021. [百度学术]

DET NORSKE VERITAS. Design of offshore wind turbine structures： DNV-OS-J101［S］. Høvik： DNV， 2005. [百度学术]

席义博，龚优华，潘嘉宁，等. 单桩基础弱化对海上风机动力响应的影响［J］. 中南大学学报（自然科学版）， 2023， 54（1）： 229-237. [百度学术]

XI Yibo， GONG Youhua， PAN Jianing， et al. Influence of monopile foundation weakening on dynamic response of offshore wind turbine［J］. Journal of Central South University（Science and Technology）， 2023， 54（1）： 229-237. [百度学术]

ZAAIJER M B. Foundation modelling to assess dynamic behaviour of offshore wind turbines［J］. Applied Ocean Research， 2006， 28（1）： 45-57. [百度学术]

BHATTACHARYA S， ADHIKARI S. Experimental validation of soil-structure interaction of offshore wind turbines［J］. Soil Dynamics and Earthquake Engineering， 2011， 31（5-6）： 805-816. [百度学术]

ADHIKARI S， BHATTACHARYA S. Dynamic analysis of wind turbine towers on flexible foundations［J］. Shock and Vibration， 2012， 19（1）： 37-56. [百度学术]

SHADLOU M， BHATTACHARYA S. Dynamic stiffness of monopiles supporting offshore wind turbine generators［J］. Soil Dynamics and Earthquake Engineering， 2016， 88： 15-32. [百度学术]

ARANY L， BHATTACHARYA S， ADHIKARI S， et al. An analytical model to predict the natural frequency of offshore wind turbines on three-spring flexible foundations using two different beam models［J］. Soil Dynamics and Earthquake Engineering， 2015， 74： 40-45. [百度学术]

ARANY L， BHATTACHARYA S， HOGAN S J， et al. Dynamic soil-structure interaction issues of offshore wind turbines［C］//Proceedings of the 9th International Conference on Structural Dynamics. Porto， Portugal： 2014： 3611-3617. [百度学术]

ARANY L， BHATTACHARYA S， MACDONALD J H G， et al. Closed form solution of eigen frequency of monopile supported offshore wind turbines in deeper waters incorporating stiffness of substructure and SSI［J］. Soil Dynamics and Earthquake Engineering， 2016， 83： 18-32. [百度学术]

杨春宝，王睿，张建民. 单桩基础型近海风机系统自振频率实用计算方法［J］. 工程力学， 2018， 35（4）： 219-225. [百度学术]

YANG Chunbao， WANG Rui， ZHANG Jianmin. Numerical method for calculating system fundamental frequencies of offshore wind turbines with monopole foundations［J］. Engineering Mechanics， 2018， 35（4）： 219-225. [百度学术]

SPYRAKOS C C， XU C. Soil-structure-water interaction of intake-outlet towers allowed to uplift［J］. Soil Dynamics and Earthquake Engineering， 1997， 16（2）： 151-159. [百度学术]

赵密，常逸夫，王丕光，等. 考虑水-桩-土相互作用的单桩式海上风机结构系统自振频率解析解［J］. 震灾防御技术， 2020， 15（4）： 659-669. [百度学术]

ZHAO Mi， CHANG Yifu， WANG Piguang， et al. An analytical natural frequency solution of monopile offshore wind turbine considering pile-soil interaction［J］. Technology for Earthquake Disaster Prevention， 2020， 15（4）： 659-669. [百度学术]

赵密，苏成坤，王丕光，等. 地震作用下水-结构-土动力相互作用分析［J］. 工程力学， 2022， 39（3）： 51-63. [百度学术]

ZHAO Mi， SU Chengkun， WANG Piguang， et al. Analysis of water-structure-soil dynamic interaction under earthquakes［J］. Engineering Mechanics， 2022， 39（3）： 51-63. [百度学术]

许成顺，孙毅龙，翟恩地，等. 海上风电单桩基础自振频率及参数影响分析［J］. 太阳能学报， 2020， 41（12）： 297-304. [百度学术]

XU Chengshun， SUN Yilong， ZHAI Endi， et al. Offshore turbine monopile foundation natural frequency and parameter impact analysis［J］. Acta Energiae Solaris Sinica， 2020， 41（12）： 297-304. [百度学术]

任文渊. 近海风电基础流固耦合三维动力响应分析［D］. 杨凌：西北农林科技大学， 2011. [百度学术]

REN Wenyuan. Analysis on 3D FSI dynamic response of offshore wind turbine foundation［D］. Yangling： Northwest A&F University， 2011. [百度学术]

田树刚，陈清军. 考虑流固耦合的近海风电支撑体系自振特性分析［J］. 力学季刊， 2014， 35（3）： 473-481. [百度学术]

TIAN Shugang， CHEN Qingjun. Analysis of dynamic characteristics of offshore wind power foundation with consideration of fluid-structure interaction［J］. Chinese Quarterly of Mechanics， 2014， 35（3）： 473-481. [百度学术]

赵家奎.微分变换及其在电路中的应用［M］. 武汉：华中理工大学出版社， 1988. [百度学术]

崔灿，蒋晗，李映辉. 变截面梁横向振动特性半解析法［J］. 振动与冲击， 2012， 31（14）： 85-88. [百度学术]

CUI Can， JIANG Han， LI Yinghui. A semi-analytical method for vibration characteristic of variable cross-section beam［J］. Journal of Vibration and Shock， 2012， 31（14）： 85-88. [百度学术]

ARIKOGLU A， OZKOL I. Solution of boundary value problems for integro-differential equations by using differential transform method［J］. Applied Mathematics and Computation， 2005， 168（2）： 1145-1158. [百度学术]

DAMGAARD M， IBSEN L B， ANDERSEN L V， et al. Cross-wind modal properties of offshore wind turbines identified by full scale testing［J］. Journal of Wind Engineering and Industrial Aerodynamics， 2013， 116： 94-108. [百度学术]

黄震. 海上风机超大直径单桩基础的桩土相互作用研究［D］. 大连：大连理工大学， 2022. [百度学术]

HUANG Zhen. An investigation of pile-soil interaction of super-large diameter monopile foundation for offshore wind turbine［D］. Dalian： Dalian University of Technology， 2022. [百度学术]

常士骠，张苏民. 工程地质手册［M］. 北京：中国建筑工业出版社， 2007. [百度学术]

CHANG Shipiao， ZHANG Sumin. Engineering Geological Handbook［M］. Beijing： China Architecture & Building Press， 2007. [百度学术]

考虑水‑桩‑土相互作用的海上风机单桩基础横向自振频率计算 PDF

摘要

关键词

1 方程的建立与求解

1.1 地基模型简化

1.2 塔筒模型简化与控制方程建立

1.3 微分变换法

1.4 频率方程建立

2 算例分析

2.1 有效性验证

2.2 参数分析

3 结 论

参考文献

3 结论