海洋平台设计原理大作业

1、SHANGHAI JIAO TONG UNIVERSITY海洋平台设计原理课程大作业姓名：王志强学号：5130109174专业：船舶与海洋工程目录1.引言22. 波浪理论22.1 波浪理论概述22.2 微幅波理论22.3 Stokes波浪理论32.4 波浪力及波浪力矩42.5 水流力53. 牛顿迭代法求解非线性方程组54. MATLAB计算实例64.1 程序流程64.2 海况等参数64.3 计算结果65. 总结86.附录：MATLAB源代码81.引言海洋覆盖着地球 3/4 的面积，海底蕴藏着丰富的油气资源，海洋已成为21 世纪人类最重要的能源基础之一。 1947 年在美国建成了世界上第一座钢结

2、构平台， 50 多年以来海底油气的开发和利用越来越受到各国重视。而开发海底油气资源，首先必须设计海洋结构物。波浪荷载是海洋结构物的主要控制荷载之一，要设计安全可靠的海洋结构物，就必须考虑波浪作用的影响。目前人们对波浪与海洋结构物相互作用的研究主要通过三种手段进行：其一是通过现场观测研究；其二是用流体力学或数理统计或能量平衡方法，在某种假设基础上，把自然界的波浪归结为某一模式， 用数学分析的方法进行研究；其三是模拟实验的方法。随着电子计算机的发展和普及，波浪的数值模拟得到了迅速的发展，它弥补了实验室模拟的不足， 而且易于实现、成本低廉， 同时也弥补了纯数学演算的抽象和失真。以数值模拟的波浪数据作

3、为输入可计算海上和海岸建筑物或船体等的响应，又由于数值模拟的可控性更强，可通过输入得到海上和海岸建筑物等长期（甚至数百年）响应的某些重要特征，如最大响应和某些临界值等。 20 世纪 80 年代以来，波浪的数值模拟与物理模拟相结合，即计算机控制下的物理模拟，已成为波浪研究的更有力的手段。随着社会经济的增长，人类对海洋的认识不断提高，利用海洋资源的能力不断增强，对海洋空间的探索也不断扩大。越来越多的领域需要对波浪进行模拟。特别在海洋工程领域，波浪的模拟已成为研究波浪特性、波浪作用的一个重要手段，因此在“海洋平台设计原理”这门课程中我也尝试采用莫里森公式计算多桩腿的波浪和水流作用力力矩。2. 波浪理

4、论2.1 波浪理论概述在海洋工程设计中，常采用的波浪理论有如下三种：(1) 微幅波理论；(2) Stokes波浪理论（二阶近似、 三阶近似、五阶近似）；(3) 流函数理论。对于微幅波理论和 Stokes 波浪理论，要计算水质点的速度和加速度，须首先知道波长，而波长需通过求解波长方程获得。微幅波理论和 Stokes二阶波浪理论波长方程相同，均是一元非线性方程。三阶和五阶 Stokes 波浪理论的波长方程均由色散关系式和附加方程组成，它们都是二元非线性方程组，由于该二元方程组表达式过于复杂， 需要进行数值分析求解。流函数理论直接假设了波面方程和水质点速度的形式， 其波剖面参数和速度参数需通过优化方

5、法获得。在这里我们主要介绍微幅波理论和Stokes波浪理论（五阶）。2.2 微幅波理论微幅波理论（ Airy 理论）是应用势函数来研究波浪运动的一种线性波浪理论，是波浪理论中最基本、最重要的内容，也是海洋工程中应用的最为广泛的波浪理论。微幅波理论的波面方程、速度势函数和色散关系式如下：波面方程：=H2coskx-t速度势： =gH2chk(z+d)chkdsin(kx-t)色散关系：2=gkthkd式中：d为水深（m）； H为波高（m）； T为波浪周期（s）； k为波数，k=2/L； 为圆频率，=2/T。2.3 Stokes波浪理论为了更准确的描述波浪运动， Stokes 提出了一种有限振幅重

6、力波的高阶理论。他的基本假定是，波浪运动能用小扰动级数表示，并且认为，考虑的量阶越高越接近实际波浪情形。这样，就得到了计入不同量阶的波浪理论，即所谓的二阶、三阶和五阶 Stokes 波浪理论等。其中二阶 Stokes 波浪理论的波面方程、速度势函数和色散关系式如下：二阶 Stokes 波浪理论波面方程： 二阶 Stokes 波浪理论速度势函数：五阶 Stokes 波理论是目前工程计算中应用广泛的波浪，与二阶、三阶Stokes 波浪理论相比，它更能反映波浪的非线性特性。 其波面方程、速度势函数和色散关系式如下：五阶 Stokes 波浪理论波面方程：k=coskx-t+2B22+4B24cos2k

7、x-t+3B33+5B35cos3kx-t+4B44cos4kx-t+5B55cos5kx-t五阶 Stokes 波浪速度势方程：kc=A11+3A13+5A15chkz+dsinkx-t+2A22+4A24ch2kz+dsin2kx-t+3A33+5A35ch3kz+dsin3kx-t+4A44ch4kz+dsin4kx-t+5A55ch5kz+dsin5kx-t波高H与波面高度之间符合下列关系：H=|=0-|=将波面高度代入到上式，得到：Hd=1dL+3B33+5B35+B55色散关系式：kc2=C02(1+2C1+4C2)其中：C02=gthkd，c=k=L/T，整理得：dL0=dLth

8、2dL1+2C1+4C2其中：L0=gT22。已知波高H、波周期T、水深d后，由于系数B33、B35、B55、C1、C2仅仅是d/L的函数，联立求解非线性方程组即可确定系数和L，然后便可得出Stokes五阶波浪理论中的其他18个系数以及波浪特征参数，由此可以确定该波浪的速度势。2.4 波浪力及波浪力矩在海洋工程实际工程应用中，当物体的尺度与波长相比是微小量的情况下，可忽略物体对波浪运动的影响，这个比值一般定为D/L0.2(其中D是物体的特征长度，如圆柱体则D是直径，L是波长)。D/L0.2的构件，一般称为小尺度构件。对于小尺度构件上的波浪力，通常采用著名的Morison公式计算。自升式平台，无

9、论桩腿是圆柱式还是析架式(可折合成圆柱式计算)都可看作是小尺度构件。Morison方程理论假定，柱体的存在对波浪运动无显著影响，认为波浪对柱体的作用主要是粘滞效应和附加质量效应。取如图所示的坐标系，莫里森公式给出，作用于单个钢桩、高dz上的水平波浪力为：dFH=fHdz12=CDDuxuxdz+CMD24uxtdz式中：CD为拖曳力系数；CM为惯性力系数。莫里森等认为作用于柱体任意高度z处的水平波浪力包括两个分量:一是波浪水质点运动的水平速度ux，引起的对柱体的作业力一水平拖曳力，另一是水质点运动的水平加速度ax引起的对柱体的作业力水平惯性力。又认为波浪作用在柱体上的拖曳力的模式与单向定常水流

10、作用在柱体上的拖曳力模式相同，即它与波浪水质点的水平速度的平方和单位柱高垂直于波向的投影面积成正比。不同的是波浪水质点作周期性的往复的振荡运动，水平速度是时正时负，因而对柱体的拖曳力也是时正时负，故在式中，取uxux代替ux了以保持拖曳力的正负性质。 整个钢桩受到的水平波浪力为：FH=-ddFH=-dfHdz=-d12CDDuxuxdz+-dCMD24uxtdz整个钢桩的总水平波力矩（对海底求矩）为：MH=-dzfHdz=-dzfD+fIdz=-d12CDDuxuxzdz+-dCMD24uxtzdz2.5 水流力取水流为剪切流，水流速度沿深度方向的变化分布由挪威船级社(DNV)推荐的公式计算：

11、Vz=V0(d+zd)式中：为速度分布指数，取为1/7。作用于单个钢桩，高dz上的水平波浪力为：dFC=fCdz=12CDDVzVzdzFC=-ddFC=-dfCdz=-d12CDDVzVzdz3. 牛顿迭代法求解非线性方程组牛顿迭代法（Newton's method）是一种在实数域和复数域上近似求解方程的方法。多数方程不存在求根公式，因此求精确根非常困难，甚至不可能，从而寻找方程的近似根就显得特别重要。方法使用函数f(x)的泰勒级数的前面几项来寻找方程f(x) = 0的根。牛顿迭代法是求方程根的重要方法之一，其最大优点是在方程f(x) = 0的单根附近具有平方收敛，而且该法还可以用来

12、求方程的重根、复根，此时线性收敛，但是可通过一些方法变成超线性收敛。另外该方法广泛用于计算机编程中。用牛顿迭代法解非线性方程，是把非线性方程f(x)=0 线性化的一种近似方法。把f(x)在点x0的某邻域内展开成泰勒级数: 取其线性部分（即泰勒展开的前两项），并令其等于0，即  以此作为非线性方程  f(x)=0 的近似方程，若   ，则其解为这样，得到牛顿迭代法的一个迭代关系式：4. MATLAB计算实例4.1 程序流程1.在MATLAB中定义函数c（x），s（x），以及A11，B35等，以备后续计算调用；2.定义

13、前面所述的波面高度以及色散关系方程组；3.设置误差范围eps以及最大迭代次数N，利用newton迭代法数值求解2中定义的非线性方程组；4.将3中所解未知数的数值作为自变量带入已定义的A11，B35等函数中计算出所有stokes五阶波中所有参数的数值，进而得到速度势方程；5.通过速度势方程得到速度方程；6.利用morison方程和5中所得速度方程计算波浪力和波浪力矩，用梯形法积分得到整个钢柱的波浪力和力矩；7.计算水流力。4.2 海况等参数计算用平台参数及海况条件等参数设置如下：参数数值水深d/m20波高H/m4.5周期T/s13钢柱直径D/m2CD1.1CM1.8V0/m/s1钢柱1位置X1/

14、m0钢柱2位置X2/m10计算时间/s30密度/kg/m310254.3 计算结果1.输入水深、波高、周期之后用newton迭代法解非线性方程组得到：=0.07930800925395002.计算30s的波浪力、波浪力矩和水流力，计算之后绘图如下：水平波浪力随时间的变化图：波浪力矩随时间变化图：水流力随时间变化图：5. 总结虽然我们专业的名称是船舶与海洋工程，但是专业必修课全部以船舶为对象的讨论，船舶与海洋工程结构物虽然有很多相似的地方，但也有很多不同之处，所以选修这门“海洋平台设计原理”对于拓展知识面非常有帮助，在这门课程中也了解了很多关于海洋平台的知识，对于未来的学习、工作都很有帮助。这门

15、课程的大作业是自己编程计算海洋平台桩腿的波浪力，在编程的过程中遇到了很多问题，比如如何解非线性方程组、如何用梯形法求积分，在查阅了很多资料之后最终一一解决了这些问题完成了大作业，也体会到了解决问题的快乐。大学四年，个人觉得我们专业在学生的实践能力培养方面做得工作还不够，如果更多的课程能有这样的大作业，而不是科普式的概述作业，对于学生的成长肯定大有裨益。最后非常感谢两位老师和助教的辛勤付出，感谢同学在我完成大作业的过程中给予的帮助。6.附录：MATLAB源代码将所有.m文件复制在matlab工作目录下运行platform.m即可开始计算function s=s(x) s=sinh(2*pi*x)

16、;endfunction c=c(x) c=cosh(2*pi*x);endfunction A11=A11(x) A11=1/s(x);endfunction A13=A13(x) A13=-c(x)2*(5*c(x)2+1)/(8*s(x)5);endfunction 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);endfunction A22=A22(x) A22=3/(8*s(x)4);endfunction A24=A24(x) A24=(192*c(

17、x)8-424*c(x)6-312*c(x)4+480*c(x)2-17)/(768*s(x)10);endfunction A33=A33(x) A33=(13-4*c(x)2)/(64*s(x)7);endfunction 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);endfunction A44=A44(x) A44=(80*c(x)6-816*c(x)4+1338*c(x)2-197)/(1536*s(x

18、)10*(6*c(x)2-1);endfunction A55=A55(x) A55=-(2880*c(x)10-72480*c(x)8+324000*c(x)6-432000*c(x)4+163470*c(x)2-16245)/(61440*s(x)11*(6*c(x)2-1)/(8*c(x)4-11*c(x)2+3);endfunction B22=B22(x) B22=(2*c(x)2+1)*c(x)/(4*s(x)3);endfunction 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(

19、x)9);endfunction B33=B33(x) B33=3*(8*c(x)6+1)/(64*s(x)6);endfunction B35=B35(x) B35=(88128*c(x)14-208224*c(x)12+70848*c(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);endfunction 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

20、)2-1);endfunction 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*c(x)2+225)/(12288*s(x)10*(6*c(x)2-1)*(8*c(x)4-11*c(x)2+3);endfunction C1=C1(x) C1=(8*c(x)4-8*c(x)2+9)/(8*s(x)4);endfunction C2=C2(x) C2=(3840*c(x)12-4096*c(x)10+2592*c(x)8-100

21、8*c(x)6+5944*c(x)4-1830*c(x)2+147)/(512*s(x)10*(6*c(x)2-1);endfunction f=F(X) syms x y H=4.5; T=13; d=20; g=9.80665; f1=-pi*H/d+(y+y3*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;endfunction df=dF(X) f=F(X); df=diff(f,'x');diff(f,'y'); %df=conj(d

22、f');end%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;density=1025;V0=1;tic;%¸Ã²¿·Ö²ÎÊýÐèÒªÌáÇ°¸ø¶¨%for i=1:N; f=double(subs(F(x0),'x' 'y&

23、#39;,x0(1) x0(2); df=double(subs(dF(x0),'x' 'y',x0(1) x0(2); Y=x0-f/df; if norm(Y-x0)<eps break; end x0=Y;end%¸Ã²¿·ÖΪnewtonµü´ú·¨¼ÆËãd/L(x)£¬ºÍ¦Ë(y)%È·

24、;¶¨¸÷¸ö²ÎÊý,Ϊ·½±ãºóÐøµ÷ÓÃ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

25、);b35=B35(Y(1);b44=B44(Y(1);b55=B55(Y(1);C1(Y(1);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

26、+d)*sin(3*(k*x-2*pi/T*t)+lambda4*a44*cosh(4*k*(z+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,'x');ax=diff(ux,'t');eta=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)+la

27、mbda4*b44*cos(4*(k*x-2*pi/T*t)+lambda5*b55*cos(5*(k*x-2*pi/T*t);all=phi ux ax eta; %½«º¯Êý»ã×Ü£¬·½±ãµ÷ÓÃ%Á¦µÄ¼ÆËãtime =0:0.1:50;force = zeros(1,length(time);mom

28、ent = zeros(1,length(time);zLength = 100;forceCurrent = zeros(1,length(time);xValue1=0; %λÖÃ1µÄ×ø±êfor 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

29、=double(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

30、; Vz = V0*(d+z)/d)*(1/7); forceCurrent(i) = forceCurrent(i)+1/2*C_D*density*D*Vz*abs(Vz); end 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) = f

31、orce(i)+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*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 = linsp