基于ABAQUS的小尺度桩柱波浪力计算方法

1、基于ABAQUS的小尺度桩柱波浪力计算方法论文导读:：基于达索公司开发的ABAQUS有限元分析软件，利用ABAQUS提供的用户子程序编制程序，计算分析风电机组在海洋环境中受到的波浪力。在计算中，使用了三维梁单元和壳单元建模，应用线性波模拟海洋环境中的波浪。基于莫里森方程和数值积分的方法，分析桩柱形支撑结构上的波浪载荷，计算支撑结构的位移，应力变化等，得到了较好的结果。论文关键词：ABAQUS，海洋风电机组，莫里森方程，线性波，波浪力1 前言能源是现代社会不可或缺的资源，是科技开展的直接推动力。目前，人类消耗能源的绝大局部是石油，天然气，煤炭等化石能源。然而，科学家早已证明，化石能源是有限的，不

2、是取之不尽，用之不竭的，地球上的煤炭，石油等资源迟早有开发殆尽的一天。相关资料显示，世界上的煤炭还能开采200年，石油和天然气还能开采50年。在这样的时代背景之下，面对着日益严重的资源危机和环境问题的挑战，人类开始寻找在传统能源之外的新型能源可再生而且环保的清洁能源。风能以其结构简单，技术成熟成为众多清洁能源中的佼佼者。目前各个国家，尤其是欧洲国家正在大力研究风能技术，开展风能产业，建设并投入使用了许多陆上和海上风电场。伴随着这股新能源的开发大潮，海上风电场的建设日益成为当前海洋工程领域的热点问题。目前我国的海上风电场建设尚处于起步阶段，国家首个风电场示范工程东海大桥风电场工程正在建设过程中，

3、风电机组的许多关键技术还在进一步探索和研究当中，海洋环境载荷分析便是关键技术之一。在海洋风电场所处的近岸海域，波浪载荷是一类重要的海洋环境载荷。当波浪与结构物发生作用后，在结构物附近，波面将发生复杂的非线性变形,理论分析较为困难，在工程中计算和分析的难度很大。在另一方面线性波，支撑结构是风电机组的根底，支撑结构的质量直接决定了风电机组的平安性，可靠性和稳定性。从投资比例来看，陆上风电机组的支撑结构投资占到总投资额的9%，海上风电机组的支撑结构投资更是占到总投资额的25%，可见其重要程度。综上所述，研究作用于风电机组支撑结构上的波浪力是风电机组开发中的核心和关键技术。因此，正确计算和分析支持结构

4、所承载的波浪载荷不仅对于波浪理论具有重要意义，还具备充分的工程应用价值。2 计算理论和方法工程上对于尺度较小的结构物，目前广泛采用Morison公式来计算结构所承载的波浪载荷。计算的根本过程是：1选择适宜的波浪理论(线性波浪理论或者视需要选择其他高阶波浪理论)；2根据雷诺数Reynolds number以及其他相关参数选择适宜的水动力系数：拖曳力系数和惯性力系数；3应用莫里森公式 Morrisons Equation。2.1 直立桩柱上的波浪力海洋风力发电机组的支撑结构存在多种形式，有单桩式，三角架式，群桩式，重力式，导管架式等多种。归纳起来，单桩式和重力式可以看作是浸入水中的直立桩柱，而三桩

5、式，群桩式可以看作多根不同角度倾斜浸入水中的桩柱的组合。因此对直立和倾斜桩柱的波浪力分析能很好地反映工程实际。图一 小尺度直立桩柱波浪力计算示意图假定有一柱体，直立在水深为d的海底上，波高为H的入射波沿着x正方向传播，圆柱体中轴线与静水面的焦点为坐标轴原点，z轴垂直向上，模型示意图如图一所示。根据挪威船级社DNV标准，位于桩柱任意高度z处的一个单位长度的微段，其收到的水平波浪力按Morison方程计算：(1)Morison公式中的和 随选择的波浪理论不同而异。在工程实际中，要针对柱体所在的海域的水深和设计波的波高，周期等条件综合考虑选用适宜的波浪理论来确定和 ，同时合理选择水动力系数和免费。本

6、文在分析和计算中采用线性波理论。线性波是指波高H与波长L，水深d相比为小量的波，线性波理论也称作Airy理论或正弦波理论。线性波的波面方程为：2)根据线性波浪理论，带入Morison方程，考虑圆柱轴心上x=0，那么单位圆柱上的水平波浪力：(3)对于规那么的圆柱体，可以得到积分式的解析解表达式。对于工程实际中的不规那么结构物，工程师常采用梯形法分段数值积分的方式求得桩柱的波浪力和力矩的数值。计算公式如式(4)所示。(4)2.2 倾斜桩柱上的波浪力图二 倾斜桩腿受力示意图假定有一柱体，倾斜在水深为d的海底上，倾斜角为，波向角为，波高为H的入射波沿着x正方向传播，z轴与静水面的焦点为坐标轴原点，模型

7、示意图如图二所示。就沿x方向传播的二维波浪的情况而论，在柱体上任一点处，水质点的速度可分解为与主轴正交和相切的水质点速度、和加速度、。在一般情况下，对于任意方向的倾斜柱体，与主轴正交的速度和加速度将不在同一空间平面内。因此由拖曳力和惯性力的总力，就必须写为矢量和形式，一般形式下的Morison方程应该写为：(5)将合力投影到三个坐标轴上，那么式矢量表达式(5)可以写为：(6)根据线性波浪理论，可推导出速度和加速度的表达式如式(7) 和式(8)所示：(7a)(7b)(7c)(7d) (8a) (8b)(8c) 为单位向量，所以，在坐标轴上的投影为：(9a)(9b)(9c)将上述表达式带入式(5)

8、中，便可以求的斜桩上某处的波浪力，然后采取分段数值积分的方法求出总的波浪力，具体积分方法与垂直桩柱相同。2.3 ABAQUS UserSubroutine对于本文所述的问题线性波，我们将模型按照梁单元建立整体模型，将波浪载荷以分布载荷的方式加载到单元上。由于在桩柱的高度范围内，波浪力的变化幅值较大，方向亦随时间而改变，普通有限元软件中的分布载荷无法满足计算要求。因此为了解决这个问题，我们引入ABAQUS User Subroutine，即用户子程序来计算波浪载荷。本文中所用到的user subroutine为DLOAD(distributed load)，用于描述分布载荷，还可用来描述海洋工程

9、中的风，浪，流以及浮力等分布载荷。这里特别要说明的是，ABAQUS在计算时，DLOAD子程序只向源程序返回分布载荷的数值大小，而不包括分布载荷的方向。对于DLOAD而言，载荷的方向始终正对着单元的外法线方向，如图三和图四所示。由于梁单元存在无穷多条法线，且倾斜桩柱模型的各个梁单元所受到的载荷方向亦各不相同，难以确定DLOAD的方向，这是计算的难点之一。本文采取的解决方法是，用Matlab编制程序，算出模型各个高度的波浪力的方向余弦，根据得到的方向余弦在INP中给每个梁单元单独规定一个法线方向。图三梁单元的法线图四 壳单元的法线3 计算结果分析支撑结构的有限元分析，就是将实际结构简化离散成为有限

10、元模型，根据实际情况划分单元和节点，运用编写的FORTRAN程序，将波浪载荷等效为分布载荷作用在各个单元上，用有限元软件计算各个节点的位移和应力，最终算出整体结构的应力状态和整体变形，本文分别用垂直和倾斜桩柱模型模拟工程实际中的支撑结构。3.1 水动力系数的选择Morison公式中的两个水动力系数，即拖曳力系数和惯性力系数的取值关系到最后计算结果的正确性。关于这两个系数的取值，国内外的学者对此作了大量的实验研究，在相关的论文和标准中给出了许多建议值供工程实际中参考。通常情况下，如果要得到比拟精确的结果，必须要在实际水域中进行水动力实验，以确定系数取值。表一引自MIT开放课程?DesignPri

11、nciples for Ocean Vehicles?第13章?SeakeepingAnalysis Natural Frequency/Natural Periods?。在本文中，依据?海港水文标准?确定和值。表一 CD和CM的推荐值 Wave Theory CD CM Comments Reference Linear Theory 1.0 0.95 Mean values for ocean wave data on 13-24in cylinders Wiegel, et al(1957) 1.0-1.4 2.0 Recommended design values based on s

12、tatistical analysis of published data Agerschou and Edens(1965) Stokes 3rd order 1.34 1.6 Mean values for oscillatory flow for 2-3in.cylinders Keulegan and Carpenter(1958) Stokes 5rd order 0.8-1.0 2.0 Recommended values based on statistical analysis of published data Agerschou and Edens(1965) 3.2 试验

13、工况本文引用上海交通大学兰雅梅所做的孤立桩柱波浪载荷试验结果，与ABAQUS有限元计算结果进行比照。桩柱模型材料为有机玻璃，测试仪器为四分力天平，模型直径为0.05m，长度为1.1m，试验水池水深为0.6m。中实验测试了不同斜率和不同波向夹角的多种工况，本文仅引述其直立桩柱计算结果。其实验工况如表二所示。表二 单桩规那么波实验工况表 波要素 No H(m) T(s) No H(m) T(s) 1 0.08 1.2 8 0.12 1.7 2 0.08 1.5 9 0.15 1.2 3 0.08 1.7 10 0.15 1.5 4 0.08 1.0 11 0.15 1.7 5 0.12 1.2 1

14、2 0.2 1.7 6 0.12 1.35 13 0.23 1.7 7 0.12 1.5 3.3 垂直桩柱计算结果在ABAQUS中由各节点的坐标生成有限元模型，整个支撑结构由梁单元构建，直立桩柱采用平面梁单元B21，底部节点坐标为( 0.0, -0.6 )，顶部节点坐标为( 0.0, 0.5 )，节点总数为111，单元总数为110。节点(0.0, -0.6 )至节点(0.0, 0.5)之间每0.01米设置一个单元，底部节点边界条件为完全固定。使用ABAQUS/STANDARD隐式求解器求解，步长为0.01s，总长为对应工况的周期。表三给出了如表二所示的13种工况下，ABAQUS计算值，实验值和

15、计算值的比拟。表三 规那么波作用下垂直桩柱上波浪力计算值与实验值比拟 工况 实验值 本文计算值 本文误差 计算值 误差 1 1.615 1.5215 5.79% 1.480 8.36% 2 1.327 1.3637 2.77% 1.329 0.15% 3 1.094 1.26229 15.38% 1.184 8.23% 4 2.447 2.50149 2.23% 2.345 4.17% 5 2.247 2.3388 4.09% 2.252 0.22% 6 2.188 2.2368 1.63% 2.138 2.29% 7 2.027 2.10724 3.96% 2.017 0.49% 8 1.8

16、58 1.9448 4.67% 1.893 0.27% 9 3.083 3.05563 0.89% 2.887 6.36% 10 2.837 2.77336 2.24% 2.634 7.16% 11 2.567 2.61732 1.96% 2.485 3.19% 12 4.66 4.11621 11.67% 4.631 0.62% 13 6.08 5.2961 12.89% 5.985 1.56% 3.4 倾斜桩柱计算结果由于倾斜桩柱情况没有适宜的文献数据可用于比拟，在本文中线性波，仅以假设的参数和工况计算一组数据结构以证明计算方法的可行性。本文参照东海大桥风电场示范工程风电机组支撑结构，设计

17、倾斜桩柱模型的参数如表四所示，计算海况如表五所示。表四 倾斜桩柱模型参数 长度(m) 外径(m) 厚度(m) 倾斜角 波向角 19.8 1.7 0.1 45 30 表五 设计波浪要素 海平面高程(m) 泥面高程(m) 水深(m) 波高(m) 周期(s) 波长(m) Cm Cd 0 -12 12 2 10 99.667 2.0 1.0 倾斜桩柱采用空间梁单元B31，底部节点坐标为( 0.0, 0.0, -12.0)，顶部节点坐标为(12.124, 7.0, 2.0 )。由于每个单元都要独立计算每个时刻的法线方向，工作量较大，故设置单元数量较少，节点总数为14，单元总数为13。同样考虑到速度和加速

18、度的变化梯度较大，将相应区域的单元加密。节点( 0.0, 0.0, -12.0)至节点( 9.526, 5.5, -1.0)之间每2米设置一个单元，节点( 9.526, 5.5, -1.0)至节点( 11.258, 6.5, 1.0)之间每0.4米设置一个单元。为了使DLOAD正确加载，在INP文件中设置了13个单元集(element set)，每个单元集对应一个单元，根据编制的Matlab程序的计算结果，每一步计算都要重新定义法线方向。使用ABAQUS/STANDARD隐式求解器求解，步长为1s，总长为20s。如图五所示，倾斜桩模型各单元的法线方向不在一个平面内，DLOAD在加载时，将为朝向

19、不同的方向，模拟了在二维波中倾斜桩柱所受波浪力的实际情况。图六所示为t=2s时刻倾斜模型的应力云图。图五 t = 2s 时倾斜模型各单元的法线图六 t = 2s 时倾斜模型的应力云图图七和图八为倾斜桩柱模型的总力图和总力矩图，其横坐标为时间t，纵坐标为作用于桩柱上的波浪力F和力矩M。图七 倾斜桩柱模型的总力图图八倾斜桩柱模型的总力矩图图九和图十为倾斜桩柱模型自由端的位移图和转角图，其横坐标为时间t，纵坐标为桩柱自由端的位移U和转角UR。图九 倾斜桩柱模型自由端位移 图十 倾斜桩柱模型自由端转角4 DLOAD应用于壳单元有限元模型采用梁单元来计算桩柱所承受的波浪载荷，可以得到较好的结果，无论是与

20、理论解还是实验值都能较好的吻合。但是梁单元也存在着一些局限性，一方面受到单元类型的限制，应用梁单元建立的模型，必须进行简化，省略实际模型中的很多细节局部。另一方面，对于变截面的桩柱结构，梁单元在模拟上存在一定困难免费。目前工程上所建立的有限元模型，包括船舶和海洋平台等，通常都是使用壳单元和实体单元来模拟真实结构物。壳单元的优点在于划分网格较为方便，可以最大程度的保存实际结构物的细节局部，在处理桩柱截面变化的问题上也简单易行。例如中第四章中建立的某大型风力发电机组的支撑结构，截面为变截面锥筒形，塔架筒壁采用壳单元模拟，塔段间的法兰盘采用实体单元模拟。使用梁单元来分析桩柱，我们可以忽略流场的变化，

21、只考虑单位桩柱上承载的波浪力即可。如果用壳单元来构建模型，必须考虑作用于桩柱外表上波浪载荷的变化，因此在应用壳单元之前，本文首先要结合流体力学的知识对流场进行分析。由流体力学可知，海上结构桩柱直立在水中，忽略流场在z方向的变化，流体均匀地向x正方向流动，我们可将这种情况归化为流体力学中的圆柱绕流问题;。二维圆柱低速定常绕流的流型只与Re数有关。在Re较小时，流场中的惯性力与粘性力相比居次要地位，圆柱上下游的流线前后对称；随着Re的增大，在圆柱反面流线开始别离并出现漩涡；随着Re进一步增大，圆柱后缘上下两侧有涡周期性地轮流脱落，形成规那么排列的涡阵线性波，这种涡阵称为卡门涡街。工程实际中的海洋结

22、构物通常都是出于湍流状态中，由前述可知，在桩柱外表上不存在稳定的压力分布。如果在分析桩柱结构波浪力时考虑湍流情况，那么超出了本文所研究的范围。所以，不失一般地，本文假设在桩柱迎浪的半个圆周上，均匀分布压力P，桩柱在x方向的合力Fx为桩柱所受波浪力，桩柱在y方向Fy合力为0，如图十一所示。下面结合前文所述Morsion公式，推导压力P与合力Fx的关系。图十一 波浪载荷作用于壳单元(10a)(10b)根据式(10)，我们将用于梁单元的Fortran程序进行修改，将原值除以圆柱直径D，那么可将此程序移植到壳单元上，亦满足合力的要求。为了验证该方法的正确性，本文以的模型为例，分别用梁单元和壳单元建立有

23、限元模型，分别在两个模型上运用对应的DLOAD子程序，用ABAQUS进行有限元计算。计算工况为表二所示工况1，即波高为0.08m，周期为1.2s。有限元模型尺寸按中模型尺寸设置。壳单元模型如所示。在本文中以立柱底部固定端的支反力值为例，比拟两种模型的计算结果。图十二 壳单元模型所示是壳单元模型底部36个节点的节点力Fxnode值，本文取前十个分析步的结果为例。将底部所有节点的节点力Fxnode值累加，得到壳单元模型底部的总支反力Fx。从表中可以看出，壳单元和梁单元所得的结果是精确相等，这也验证了采用如本文所假设的载荷分布方式满足目前标准的要求。表六 两种单元类型的计算结果 STEP1 STEP

24、2 STEP3 STEP4 STEP5 STEP6 STEP7 STEP8 STEP9 STEP10 NODE1 3.252E-02 2.330E-02 1.391E-02 4.369E-03 -5.278E-03 -1.500E-02 -2.458E-02 -3.409E-02 -4.356E-02 -5.296E-02 NODE2 3.101E-02 2.219E-02 1.321E-02 4.094E-03 -5.129E-03 -1.442E-02 -2.358E-02 -3.267E-02 -4.173E-02 -5.071E-02 NODE3 2.666E-02 1.900E-02

25、 1.121E-02 3.301E-03 -4.700E-03 -1.276E-02 -2.071E-02 -2.860E-02 -3.645E-02 -4.424E-02 NODE4 1.999E-02 1.412E-02 8.149E-03 2.087E-03 -4.041E-03 -1.022E-02 -1.630E-02 -2.235E-02 -2.836E-02 -3.433E-02 NODE5 1.181E-02 8.133E-03 4.392E-03 5.988E-04 -3.233E-03 -7.089E-03 -1.090E-02 -1.468E-02 -1.844E-02

26、-2.217E-02 NODE6 3.110E-03 1.760E-03 3.944E-04 -9.840E-04 -2.370E-03 -3.760E-03 -5.142E-03 -6.514E-03 -7.876E-03 -9.220E-03 NODE7 -5.068E-03 -4.228E-03 -3.361E-03 -2.470E-03 -1.558E-03 -6.288E-04 2.699E-04 1.161E-03 2.056E-03 2.952E-03 NODE8 -1.174E-02 -9.111E-03 -6.423E-03 -3.681E-03 -8.946E-04 1.9

27、25E-03 4.683E-03 7.420E-03 1.016E-02 1.288E-02 NODE9 -1.609E-02 -1.230E-02 -8.421E-03 -4.472E-03 -4.630E-04 3.590E-03 7.562E-03 1.150E-02 1.544E-02 1.935E-02 NODE10 -1.760E-02 -1.340E-02 -9.119E-03 -4.754E-03 -3.241E-04 4.154E-03 8.543E-03 1.290E-02 1.725E-02 2.157E-02 NODE11 -1.607E-02 -1.230E-02 -

28、8.433E-03 -4.497E-03 -5.004E-04 3.540E-03 7.500E-03 1.143E-02 1.535E-02 1.925E-02 NODE12 -1.171E-02 -9.110E-03 -6.447E-03 -3.730E-03 -9.683E-04 1.827E-03 4.561E-03 7.274E-03 9.985E-03 1.268E-02 NODE13 -5.032E-03 -4.228E-03 -3.397E-03 -2.542E-03 -1.666E-03 -7.725E-04 9.059E-05 9.463E-04 1.807E-03 2.6

29、68E-03 NODE14 3.156E-03 1.761E-03 3.487E-04 -1.076E-03 -2.509E-03 -3.944E-03 -5.372E-03 -6.790E-03 -8.196E-03 -9.584E-03 NODE15 1.187E-02 8.134E-03 4.337E-03 4.890E-04 -3.398E-03 -7.309E-03 -1.117E-02 -1.501E-02 -1.882E-02 -2.260E-02 NODE16 2.005E-02 1.412E-02 8.087E-03 1.963E-03 -4.228E-03 -1.046E-

30、02 -1.661E-02 -2.272E-02 -2.880E-02 -3.482E-02 NODE17 2.672E-02 1.901E-02 1.115E-02 3.166E-03 -4.902E-03 -1.303E-02 -2.104E-02 -2.900E-02 -3.692E-02 -4.478E-02 NODE18 3.108E-02 2.219E-02 1.314E-02 3.953E-03 -5.341E-03 -1.471E-02 -2.394E-02 -3.310E-02 -4.222E-02 -5.127E-02 NODE19 3.259E-02 2.330E-02

31、1.383E-02 4.226E-03 -5.493E-03 -1.529E-02 -2.494E-02 -3.452E-02 -4.406E-02 -5.353E-02 NODE20 3.108E-02 2.219E-02 1.314E-02 3.953E-03 -5.341E-03 -1.471E-02 -2.394E-02 -3.310E-02 -4.222E-02 -5.127E-02 NODE21 2.672E-02 1.901E-02 1.115E-02 3.166E-03 -4.902E-03 -1.303E-02 -2.104E-02 -2.900E-02 -3.692E-02

32、 -4.478E-02 NODE22 2.005E-02 1.412E-02 8.087E-03 1.963E-03 -4.228E-03 -1.046E-02 -1.661E-02 -2.272E-02 -2.880E-02 -3.482E-02 NODE23 1.187E-02 8.134E-03 4.337E-03 4.890E-04 -3.398E-03 -7.309E-03 -1.117E-02 -1.501E-02 -1.882E-02 -2.260E-02 NODE24 3.156E-03 1.761E-03 3.487E-04 -1.076E-03 -2.509E-03 -3.

33、944E-03 -5.372E-03 -6.790E-03 -8.196E-03 -9.584E-03 NODE25 -5.032E-03 -4.228E-03 -3.397E-03 -2.542E-03 -1.666E-03 -7.725E-04 9.059E-05 9.463E-04 1.807E-03 2.668E-03 NODE26 -1.171E-02 -9.110E-03 -6.447E-03 -3.730E-03 -9.683E-04 1.827E-03 4.561E-03 7.274E-03 9.985E-03 1.268E-02 NODE27 -1.607E-02 -1.23

34、0E-02 -8.433E-03 -4.497E-03 -5.004E-04 3.540E-03 7.500E-03 1.143E-02 1.535E-02 1.925E-02 NODE28 -1.760E-02 -1.340E-02 -9.119E-03 -4.754E-03 -3.241E-04 4.154E-03 8.543E-03 1.290E-02 1.725E-02 2.157E-02 NODE29 -1.609E-02 -1.230E-02 -8.421E-03 -4.472E-03 -4.630E-04 3.590E-03 7.562E-03 1.150E-02 1.544E-02 1.935E-02 NODE30 -1.174E-02 -9.111E-03 -6.423E-03 -3.681E-03 -8.946E-04 1.925E-03 4.683E-03 7.420E-03 1.016E-02 1.288E-02 NODE31 -5.068E-03 -4.228E-03 -3.361E-03 -2.470E-03 -1.558E-03 -6.288E-04 2.699E-04 1.161E-03 2.056E-03 2.952E-03 NODE32 3.110