海洋平台设计原理大作业文档格式.docx

1、2.2 微幅波理论微幅波理论（ Airy 理论）是应用势函数来研究波浪运动的一种线性波浪理论，是波浪理论中最基本、最重要的内容，也是海洋工程中应用的最为广泛的波浪理论。微幅波理论的波面方程、速度势函数和色散关系式如下：波面方程：速度势：色散关系：式中：d为水深（m）； H为波高（m）； T为波浪周期（s）； k为波数，； 为圆频率，。2.3 Stokes波浪理论为了更准确的描述波浪运动， Stokes 提出了一种有限振幅重力波的高阶理论。他的基本假定是，波浪运动能用小扰动级数表示，并且认为，考虑的量阶越高越接近实际波浪情形。这样，就得到了计入不同量阶的波浪理论，即所谓的二阶、三阶和五阶 Sto

2、kes 波浪理论等。其中二阶 Stokes 波浪理论的波面方程、速度势函数和色散关系式如下：二阶 Stokes 波浪理论波面方程：二阶 Stokes 波浪理论速度势函数：五阶 Stokes 波理论是目前工程计算中应用广泛的波浪，与二阶、三阶Stokes 波浪理论相比，它更能反映波浪的非线性特性。 其波面方程、速度势函数和色散关系式如下：五阶 Stokes 波浪理论波面方程：五阶 Stokes 波浪速度势方程：波高H与波面高度之间符合下列关系：将波面高度代入到上式，得到：色散关系式：其中：，整理得：已知波高H、波周期T、水深d后，由于系数、仅仅是d/L的函数，联立求解非线性方程组即可确定系数和L

3、，然后便可得出Stokes五阶波浪理论中的其他18个系数以及波浪特征参数，由此可以确定该波浪的速度势。2.4 波浪力及波浪力矩在海洋工程实际工程应用中，当物体的尺度与波长相比是微小量的情况下，可忽略物体对波浪运动的影响，这个比值一般定为（其中D是物体的特征长度，如圆柱体则D是直径，L是波长）。的构件，一般称为小尺度构件。对于小尺度构件上的波浪力，通常采用著名的Morison公式计算。自升式平台，无论桩腿是圆柱式还是析架式（可折合成圆柱式计算）都可看作是小尺度构件。Morison方程理论假定，柱体的存在对波浪运动无显著影响，认为波浪对柱体的作用主要是粘滞效应和附加质量效应。取如图所示的坐标系，莫

4、里森公式给出，作用于单个钢桩、高dz上的水平波浪力为：为拖曳力系数；为惯性力系数。莫里森等认为作用于柱体任意高度z处的水平波浪力包括两个分量:一是波浪水质点运动的水平速度，引起的对柱体的作业力一水平拖曳力，另一是水质点运动的水平加速度引起的对柱体的作业力水平惯性力。又认为波浪作用在柱体上的拖曳力的模式与单向定常水流作用在柱体上的拖曳力模式相同，即它与波浪水质点的水平速度的平方和单位柱高垂直于波向的投影面积成正比。不同的是波浪水质点作周期性的往复的振荡运动，水平速度是时正时负，因而对柱体的拖曳力也是时正时负，故在式中，取代替了以保持拖曳力的正负性质。 整个钢桩受到的水平波浪力为：整个钢桩的总水平

5、波力矩（对海底求矩）为：2.5 水流力取水流为剪切流，水流速度沿深度方向的变化分布由挪威船级社（DNV）推荐的公式计算：为速度分布指数，取为1/7。作用于单个钢桩，高dz上的水平波浪力为：3. 牛顿迭代法求解非线性方程组牛顿迭代法（Newtons method）是一种在实数域和复数域上近似求解方程的方法。多数方程不存在求根公式，因此求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。方法使用函数f（x）的泰勒级数的前面几项来寻找方程f（x） = 0的根。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f（x） = 0的单根附近具有平方收敛，而且该法还可以用来求方程的重根、复

6、根，此时线性收敛，但是可通过一些方法变成超线性收敛。另外该方法广泛用于计算机编程中。用牛顿迭代法解非线性方程，是把非线性方程f（x）=0线性化的一种近似方法。把f（x）在点x0的某邻域内展开成泰勒级数:取其线性部分（即泰勒展开的前两项），并令其等于0，即以此作为非线性方程 f（x）=0的近似方程，若 ，则其解为这样，得到牛顿迭代法的一个迭代关系式：4. MATLAB计算实例4.1 程序流程1.在MATLAB中定义函数c（x），s（x），以及A11，B35等，以备后续计算调用；2.定义前面所述的波面高度以及色散关系方程组；3.设置误差范围eps以及最大迭代次数N，利用newton迭代法数值求解2

7、中定义的非线性方程组；4.将3中所解未知数的数值作为自变量带入已定义的A11，B35等函数中计算出所有stokes五阶波中所有参数的数值，进而得到速度势方程；5.通过速度势方程得到速度方程；6.利用morison方程和5中所得速度方程计算波浪力和波浪力矩，用梯形法积分得到整个钢柱的波浪力和力矩；7.计算水流力。4.2 海况等参数计算用平台参数及海况条件等参数设置如下：参数数值水深d/m20波高H/m4.5周期T/s13钢柱直径D/m21.11.8V0/m/s1钢柱1位置X1/m钢柱2位置X2/m10计算时间/s30密度/kg/m310254.3 计算结果1.输入水深、波高、周期之后用newto

8、n迭代法解非线性方程组得到：2.计算30s的波浪力、波浪力矩和水流力，计算之后绘图如下：水平波浪力随时间的变化图：波浪力矩随时间变化图：水流力随时间变化图：5. 总结虽然我们专业的名称是船舶与海洋工程，但是专业必修课全部以船舶为对象的讨论，船舶与海洋工程结构物虽然有很多相似的地方，但也有很多不同之处，所以选修这门“海洋平台设计原理”对于拓展知识面非常有帮助，在这门课程中也了解了很多关于海洋平台的知识，对于未来的学习、工作都很有帮助。这门课程的大作业是自己编程计算海洋平台桩腿的波浪力，在编程的过程中遇到了很多问题，比如如何解非线性方程组、如何用梯形法求积分，在查阅了很多资料之后最终一一解决了这些

9、问题完成了大作业，也体会到了解决问题的快乐。大学四年，个人觉得我们专业在学生的实践能力培养方面做得工作还不够，如果更多的课程能有这样的大作业，而不是科普式的概述作业，对于学生的成长肯定大有裨益。最后非常感谢两位老师和助教的辛勤付出，感谢同学在我完成大作业的过程中给予的帮助。6.附录：MATLAB源代码将所有.m文件复制在matlab工作目录下运行platform.m即可开始计算function s=s（x） s=sinh（2*pi*x）;endfunction c=c（x） c=cosh（2*pi*x）;function A11=A11（x） A11=1/s（x）;function A13=A

10、13（x） A13=-c（x）2*（5*c（x）2+1）/（8*s（x）5）;function A15=A15（x） A15=-（1184*c（x）10-1440*c（x）2-1993*c（x）6+2641*c（x）4-249*c（x）2+18）/（1536*s（x）11）;function A22=A22（x） A22=3/（8*s（x）4）;function A24=A24（x） A24=（192*c（x）8-424*c（x）6-312*c（x）4+480*c（x）2-17）/（768*s（x）10）;function A33=A33（x） A33=（13-4*c（x）2）/（64*s（x

11、）7）;function A35=A35（x） A35=（512*c（x）12+4224*c（x）10-6800*c（x）8-1280*c（x）6+16704*c（x）4-3154*c（x）2+107）/（4096*s（x）13*（6*c（x）2-1）;function A44=A44（x） A44=（80*c（x）6-816*c（x）4+1338*c（x）2-197）/（1536*s（x）10*（6*c（x）2-1）;function A55=A55（x） A55=-（2880*c（x）10-72480*c（x）8+324000*c（x）6-432000*c（x）4+163470*c（x）2

12、-16245）/（61440*s（x）11*（6*c（x）2-1）/（8*c（x）4-11*c（x）2+3）;function B22=B22（x） B22=（2*c（x）2+1）*c（x）/（4*s（x）3）;function B24=B24（x） B24=（272*c（x）8-504*c（x）6-192*c（x）4+322*c（x）2+21）*c（x）/（384*s（x）9）;function B33=B33（x） B33=3*（8*c（x）6+1）/（64*s（x）6）;function B35=B35（x） B35=（88128*c（x）14-208224*c（x）12+70848*c

13、（x）10+54000*c（x）8-218*c（x）6+6264*c（x）4-54*c（x）2-81）/（12288*s（x）12*（6*c（x）2-1）;function B44=B44（x） B44=（768*c（x）10-448*c（x）8-48*c（x）6+48*c（x）4+106*c（x）2-21）*c（x）/（384*s（x）9*（6*c（x）2-1）;function B55=B55（x） B55=（19200*c（x）16-262720*c（x）14+83680*c（x）12+20160*c（x）10-7280*c（x）8+7160*c（x）6-1800*c（x）4-1050*

14、c（x）2+225）/（12288*s（x）10*（6*c（x）2-1）*（8*c（x）4-11*c（x）2+3）;function C1=C1（x） C1=（8*c（x）4-8*c（x）2+9）/（8*s（x）4）;function C2=C2（x） C2=（3840*c（x）12-4096*c（x）10+2592*c（x）8-1008*c（x）6+5944*c（x）4-1830*c（x）2+147）/（512*s（x）10*（6*c（x）2-1）;function f=F（X） syms x y H=4.5; T=13; d=20; g=9.80665; f1=-pi*H/d+（y+y3*

15、B33（x）+y5*（B35（x）+B55（x）; f2=x*tanh（x）*（1+y2*C1（x）+y4*C2（x）-2*pi*d/（g*T2）; f=f1 f2;function df=dF（X） f=F（X）; df=diff（f,x）;diff（f,y）; %df=conj（df%H=4.5;T=13;d=20;g=9.80665;D=2;n=1;C_D=1.1;C_M=1.8;x0=0.17 0.06;eps=0.000000000001;N=100000000000000;density=1025;V0=1;tic;%for i=1:N; f=double（subs（F（x0）,

16、,x0（1） x0（2）; df=double（subs（dF（x0）, Y=x0-f/df; if norm（Y-x0）eps break; end x0=Y;newtond/L（x）（y）%,a11=A11（Y（1）;a13=A13（Y（1）;a15=A15（Y（1）;a22=A22（Y（1）;a24=A24（Y（1）;a33=A33（Y（1）;a35=A35（Y（1）;a44=A44（Y（1）;a55=A55（Y（1）;b22=B22（Y（1）;b24=B24（Y（1）;b33=B33（Y（1）;b35=B35（Y（1）;b44=B44（Y（1）;b55=B55（Y（1）;C1（Y（1）

17、;C2（Y（1）;L=d/Y（1）;lambda=Y（2）;k=2*pi/L;speed=L/T;omega=2*pi/T;syms z t xphi=speed/k*（lambda*a11+lambda3*a13+lambda5*a15）*cosh（k*（z+d）*sin（k*x-2*pi/T*t）+（lambda2*a22+lambda4*a24）*cosh（2*k*（z+d）*sin（2*（k*x-2*pi/T*t）+（lambda3*a33+lambda5*a35）*cosh（3*k*（z+d）*sin（3*（k*x-2*pi/T*t）+lambda4*a44*cosh（4*k*（z+

18、d）*sin（4*（k*x-2*pi/T*t）+lambda5*a55*cosh（5*k*（z+d）*sin（5*（k*x-2*pi/T*t）;ux=diff（phi,ax=diff（ux,teta=1/k*（lambda*cos（k*x-2*pi/T*t）+（lambda2*b22+lambda4*b24）*cos（2*（k*x-2*pi/T*t）+（lambda3*b33+lambda5*b35）*cos（3*（k*x-2*pi/T*t）+lambda4*b44*cos（4*（k*x-2*pi/T*t）+lambda5*b55*cos（5*（k*x-2*pi/T*t）;all=phi ux

19、 ax eta; %time =0:0.1:50;force = zeros（1,length（time）;moment = zeros（1,length（time）;zLength = 100;forceCurrent = zeros（1,length（time）;xValue1=0; %1for i = 1:length（time） tValue = time（i）; eta1=double（subs（eta,x,t,xValue1,tValue）; verticalPosition1 = linspace（-d ,eta1, zLength）; for j=1:zLength ux=do

20、uble（subs（ux,x,z,t,xValue1,verticalPosition1（j）,tValue）; ax=double（subs（ax,x,z,t,xValue1,verticalPosition1（j）,tValue）; z = verticalPosition1（j）; force（i） = force（i）+1/2*C_D*density*D*ux*abs（ux）+C_M*density*pi*D2/4*ax; moment（i） = moment（i） + 1/2*C_D*density*D*ux*abs（ux）*z+C_M*density*pi*D2/4*ax*z; V

21、z = V0*（d+z）/d）*（1/7）; forceCurrent（i） = forceCurrent（i）+1/2*C_D*density*D*Vz*abs（Vz）; uBottom=double（subs（ux,x,z,t,xValue1,-d,tValue）; aBottom=double（subs（ax,x,z,t,xValue1,-d,tValue）; uTop=double（subs（ux,x,z,t,xValue1,eta1,tValue）; aTop=double（subs（ax,x,z,t,xValue1,eta1,tValue）; force（i） = force（i）

22、+1/2*（1/2*C_D*density*D*uBottom*abs（uBottom）+C_M*density*pi*D2/4*aBottom+1/2*C_D*density*D*uTop*abs（uTop）+C_M*density*pi*D2/4*aTop）; moment（i） = moment（i）+1/2*（1/2*C_D*density*D*uBottom*abs（uBottom）*（-d）+C_M*density*pi*D2/4*aBottom*（-d）+1/2*C_D*density*D*uTop*abs（uTop）*eta1+C_M*density*pi*D2/4*aTop*

23、eta1）; VTop = V0*（d+eta1）/d）*（1/7）; forceCurrent（i） = forceCurrent（i）+1/2*1/2*C_D*density*D*VTop*abs（VTop）; % xValue2=10; eta2=double（subs（eta,x,t,xValue2,tValue）; verticalPosition2 = linspace（-d ,eta2, zLength）; ux=double（subs（ux,x,z,t,xValue1,verticalPosition2（j）,tValue）; ax=double（subs（ax,x,z,t,xValue1,verticalPosition2（j）,tValue）; z = verticalPosition2（j）; uTop=double（subs（ux,x,z,t,xValue1,eta2,tValue）; aTop=double（subs（ax,x,z,t,xValue1,eta2,tValue）; moment（i） = moment（i）+1/2*（1/2*C_D*density*D*uBottom*abs（uBottom）*（-d）+C_M*density*pi*D2/4*aBott