海洋平臺(tái)設(shè)計(jì)原理大作業(yè)

1、SHANGHAI JIAO TONG UNIVERSITY海洋平臺(tái)設(shè)計(jì)原理課程大作業(yè)姓名：王志強(qiáng)學(xué)號(hào)：5130109174專業(yè)：船舶與海洋工程目錄1.引言22. 波浪理論22.1 波浪理論概述22.2 微幅波理論22.3 Stokes波浪理論32.4 波浪力及波浪力矩42.5 水流力53. 牛頓迭代法求解非線性方程組54. MATLAB計(jì)算實(shí)例64.1 程序流程64.2 海況等參數(shù)64.3 計(jì)算結(jié)果65. 總結(jié)86.附錄：MATLAB源代碼81.引言海洋覆蓋著地球 3/4 的面積，海底蘊(yùn)藏著豐富的油氣資源，海洋已成為21 世紀(jì)人類最重要的能源基礎(chǔ)之一。 1947 年在美國(guó)建成了世界上第一座鋼結(jié)

2、構(gòu)平臺(tái)， 50 多年以來海底油氣的開發(fā)和利用越來越受到各國(guó)重視。而開發(fā)海底油氣資源，首先必須設(shè)計(jì)海洋結(jié)構(gòu)物。波浪荷載是海洋結(jié)構(gòu)物的主要控制荷載之一，要設(shè)計(jì)安全可靠的海洋結(jié)構(gòu)物，就必須考慮波浪作用的影響。目前人們對(duì)波浪與海洋結(jié)構(gòu)物相互作用的研究主要通過三種手段進(jìn)行：其一是通過現(xiàn)場(chǎng)觀測(cè)研究；其二是用流體力學(xué)或數(shù)理統(tǒng)計(jì)或能量平衡方法，在某種假設(shè)基礎(chǔ)上，把自然界的波浪歸結(jié)為某一模式， 用數(shù)學(xué)分析的方法進(jìn)行研究；其三是模擬實(shí)驗(yàn)的方法。隨著電子計(jì)算機(jī)的發(fā)展和普及，波浪的數(shù)值模擬得到了迅速的發(fā)展，它彌補(bǔ)了實(shí)驗(yàn)室模擬的不足， 而且易于實(shí)現(xiàn)、成本低廉， 同時(shí)也彌補(bǔ)了純數(shù)學(xué)演算的抽象和失真。以數(shù)值模擬的波浪數(shù)據(jù)作

3、為輸入可計(jì)算海上和海岸建筑物或船體等的響應(yīng)，又由于數(shù)值模擬的可控性更強(qiáng)，可通過輸入得到海上和海岸建筑物等長(zhǎng)期（甚至數(shù)百年）響應(yīng)的某些重要特征，如最大響應(yīng)和某些臨界值等。 20 世紀(jì) 80 年代以來，波浪的數(shù)值模擬與物理模擬相結(jié)合，即計(jì)算機(jī)控制下的物理模擬，已成為波浪研究的更有力的手段。隨著社會(huì)經(jīng)濟(jì)的增長(zhǎng)，人類對(duì)海洋的認(rèn)識(shí)不斷提高，利用海洋資源的能力不斷增強(qiáng)，對(duì)海洋空間的探索也不斷擴(kuò)大。越來越多的領(lǐng)域需要對(duì)波浪進(jìn)行模擬。特別在海洋工程領(lǐng)域，波浪的模擬已成為研究波浪特性、波浪作用的一個(gè)重要手段，因此在“海洋平臺(tái)設(shè)計(jì)原理”這門課程中我也嘗試采用莫里森公式計(jì)算多樁腿的波浪和水流作用力力矩。2. 波浪理

4、論2.1 波浪理論概述在海洋工程設(shè)計(jì)中，常采用的波浪理論有如下三種：(1) 微幅波理論；(2) Stokes波浪理論（二階近似、 三階近似、五階近似）；(3) 流函數(shù)理論。對(duì)于微幅波理論和 Stokes 波浪理論，要計(jì)算水質(zhì)點(diǎn)的速度和加速度，須首先知道波長(zhǎng)，而波長(zhǎng)需通過求解波長(zhǎng)方程獲得。微幅波理論和 Stokes二階波浪理論波長(zhǎng)方程相同，均是一元非線性方程。三階和五階 Stokes 波浪理論的波長(zhǎng)方程均由色散關(guān)系式和附加方程組成，它們都是二元非線性方程組，由于該二元方程組表達(dá)式過于復(fù)雜， 需要進(jìn)行數(shù)值分析求解。流函數(shù)理論直接假設(shè)了波面方程和水質(zhì)點(diǎn)速度的形式， 其波剖面參數(shù)和速度參數(shù)需通過優(yōu)化方

5、法獲得。在這里我們主要介紹微幅波理論和Stokes波浪理論（五階）。2.2 微幅波理論微幅波理論（ Airy 理論）是應(yīng)用勢(shì)函數(shù)來研究波浪運(yùn)動(dòng)的一種線性波浪理論，是波浪理論中最基本、最重要的內(nèi)容，也是海洋工程中應(yīng)用的最為廣泛的波浪理論。微幅波理論的波面方程、速度勢(shì)函數(shù)和色散關(guān)系式如下：波面方程：=H2coskx-t速度勢(shì)： =gH2chk(z+d)chkdsin(kx-t)色散關(guān)系：2=gkthkd式中：d為水深（m）； H為波高（m）； T為波浪周期（s）； k為波數(shù)，k=2/L； 為圓頻率，=2/T。2.3 Stokes波浪理論為了更準(zhǔn)確的描述波浪運(yùn)動(dòng)， Stokes 提出了一種有限振幅重

6、力波的高階理論。他的基本假定是，波浪運(yùn)動(dòng)能用小擾動(dòng)級(jí)數(shù)表示，并且認(rèn)為，考慮的量階越高越接近實(shí)際波浪情形。這樣，就得到了計(jì)入不同量階的波浪理論，即所謂的二階、三階和五階 Stokes 波浪理論等。其中二階 Stokes 波浪理論的波面方程、速度勢(shì)函數(shù)和色散關(guān)系式如下：二階 Stokes 波浪理論波面方程： 二階 Stokes 波浪理論速度勢(shì)函數(shù)：五階 Stokes 波理論是目前工程計(jì)算中應(yīng)用廣泛的波浪，與二階、三階Stokes 波浪理論相比，它更能反映波浪的非線性特性。 其波面方程、速度勢(shì)函數(shù)和色散關(guān)系式如下：五階 Stokes 波浪理論波面方程：k=coskx-t+2B22+4B24cos2k

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

8、2dL1+2C1+4C2其中：L0=gT22。已知波高H、波周期T、水深d后，由于系數(shù)B33、B35、B55、C1、C2僅僅是d/L的函數(shù)，聯(lián)立求解非線性方程組即可確定系數(shù)和L，然后便可得出Stokes五階波浪理論中的其他18個(gè)系數(shù)以及波浪特征參數(shù)，由此可以確定該波浪的速度勢(shì)。2.4 波浪力及波浪力矩在海洋工程實(shí)際工程應(yīng)用中，當(dāng)物體的尺度與波長(zhǎng)相比是微小量的情況下，可忽略物體對(duì)波浪運(yùn)動(dòng)的影響，這個(gè)比值一般定為D/L0.2(其中D是物體的特征長(zhǎng)度，如圓柱體則D是直徑，L是波長(zhǎng))。D/L0.2的構(gòu)件，一般稱為小尺度構(gòu)件。對(duì)于小尺度構(gòu)件上的波浪力，通常采用著名的Morison公式計(jì)算。自升式平臺(tái)，無

9、論樁腿是圓柱式還是析架式(可折合成圓柱式計(jì)算)都可看作是小尺度構(gòu)件。Morison方程理論假定，柱體的存在對(duì)波浪運(yùn)動(dòng)無顯著影響，認(rèn)為波浪對(duì)柱體的作用主要是粘滯效應(yīng)和附加質(zhì)量效應(yīng)。取如圖所示的坐標(biāo)系，莫里森公式給出，作用于單個(gè)鋼樁、高dz上的水平波浪力為：dFH=fHdz12=CDDuxuxdz+CMD24uxtdz式中：CD為拖曳力系數(shù)；CM為慣性力系數(shù)。莫里森等認(rèn)為作用于柱體任意高度z處的水平波浪力包括兩個(gè)分量:一是波浪水質(zhì)點(diǎn)運(yùn)動(dòng)的水平速度ux，引起的對(duì)柱體的作業(yè)力一水平拖曳力，另一是水質(zhì)點(diǎn)運(yùn)動(dòng)的水平加速度ax引起的對(duì)柱體的作業(yè)力水平慣性力。又認(rèn)為波浪作用在柱體上的拖曳力的模式與單向定常水流

10、作用在柱體上的拖曳力模式相同，即它與波浪水質(zhì)點(diǎn)的水平速度的平方和單位柱高垂直于波向的投影面積成正比。不同的是波浪水質(zhì)點(diǎn)作周期性的往復(fù)的振蕩運(yùn)動(dòng)，水平速度是時(shí)正時(shí)負(fù)，因而對(duì)柱體的拖曳力也是時(shí)正時(shí)負(fù)，故在式中，取uxux代替ux了以保持拖曳力的正負(fù)性質(zhì)。 整個(gè)鋼樁受到的水平波浪力為：FH=-ddFH=-dfHdz=-d12CDDuxuxdz+-dCMD24uxtdz整個(gè)鋼樁的總水平波力矩（對(duì)海底求矩）為：MH=-dzfHdz=-dzfD+fIdz=-d12CDDuxuxzdz+-dCMD24uxtzdz2.5 水流力取水流為剪切流，水流速度沿深度方向的變化分布由挪威船級(jí)社(DNV)推薦的公式計(jì)算：

11、Vz=V0(d+zd)式中：為速度分布指數(shù)，取為1/7。作用于單個(gè)鋼樁，高dz上的水平波浪力為：dFC=fCdz=12CDDVzVzdzFC=-ddFC=-dfCdz=-d12CDDVzVzdz3. 牛頓迭代法求解非線性方程組牛頓迭代法（Newton's method）是一種在實(shí)數(shù)域和復(fù)數(shù)域上近似求解方程的方法。多數(shù)方程不存在求根公式，因此求精確根非常困難，甚至不可能，從而尋找方程的近似根就顯得特別重要。方法使用函數(shù)f(x)的泰勒級(jí)數(shù)的前面幾項(xiàng)來尋找方程f(x) = 0的根。牛頓迭代法是求方程根的重要方法之一，其最大優(yōu)點(diǎn)是在方程f(x) = 0的單根附近具有平方收斂，而且該法還可以用來

12、求方程的重根、復(fù)根，此時(shí)線性收斂，但是可通過一些方法變成超線性收斂。另外該方法廣泛用于計(jì)算機(jī)編程中。用牛頓迭代法解非線性方程，是把非線性方程f(x)=0 線性化的一種近似方法。把f(x)在點(diǎn)x0的某鄰域內(nèi)展開成泰勒級(jí)數(shù): 取其線性部分（即泰勒展開的前兩項(xiàng)），并令其等于0，即  以此作為非線性方程  f(x)=0 的近似方程，若   ，則其解為這樣，得到牛頓迭代法的一個(gè)迭代關(guān)系式：4. MATLAB計(jì)算實(shí)例4.1 程序流程1.在MATLAB中定義函數(shù)c（x），s（x），以及A11，B35等，以備后續(xù)計(jì)算調(diào)用；2.定義

13、前面所述的波面高度以及色散關(guān)系方程組；3.設(shè)置誤差范圍eps以及最大迭代次數(shù)N，利用newton迭代法數(shù)值求解2中定義的非線性方程組；4.將3中所解未知數(shù)的數(shù)值作為自變量帶入已定義的A11，B35等函數(shù)中計(jì)算出所有stokes五階波中所有參數(shù)的數(shù)值，進(jìn)而得到速度勢(shì)方程；5.通過速度勢(shì)方程得到速度方程；6.利用morison方程和5中所得速度方程計(jì)算波浪力和波浪力矩，用梯形法積分得到整個(gè)鋼柱的波浪力和力矩；7.計(jì)算水流力。4.2 海況等參數(shù)計(jì)算用平臺(tái)參數(shù)及海況條件等參數(shù)設(shè)置如下：參數(shù)數(shù)值水深d/m20波高H/m4.5周期T/s13鋼柱直徑D/m2CD1.1CM1.8V0/m/s1鋼柱1位置X1/

14、m0鋼柱2位置X2/m10計(jì)算時(shí)間/s30密度/kg/m310254.3 計(jì)算結(jié)果1.輸入水深、波高、周期之后用newton迭代法解非線性方程組得到：=0.07930800925395002.計(jì)算30s的波浪力、波浪力矩和水流力，計(jì)算之后繪圖如下：水平波浪力隨時(shí)間的變化圖：波浪力矩隨時(shí)間變化圖：水流力隨時(shí)間變化圖：5. 總結(jié)雖然我們專業(yè)的名稱是船舶與海洋工程，但是專業(yè)必修課全部以船舶為對(duì)象的討論，船舶與海洋工程結(jié)構(gòu)物雖然有很多相似的地方，但也有很多不同之處，所以選修這門“海洋平臺(tái)設(shè)計(jì)原理”對(duì)于拓展知識(shí)面非常有幫助，在這門課程中也了解了很多關(guān)于海洋平臺(tái)的知識(shí)，對(duì)于未來的學(xué)習(xí)、工作都很有幫助。這門

15、課程的大作業(yè)是自己編程計(jì)算海洋平臺(tái)樁腿的波浪力，在編程的過程中遇到了很多問題，比如如何解非線性方程組、如何用梯形法求積分，在查閱了很多資料之后最終一一解決了這些問題完成了大作業(yè)，也體會(huì)到了解決問題的快樂。大學(xué)四年，個(gè)人覺得我們專業(yè)在學(xué)生的實(shí)踐能力培養(yǎng)方面做得工作還不夠，如果更多的課程能有這樣的大作業(yè)，而不是科普式的概述作業(yè)，對(duì)于學(xué)生的成長(zhǎng)肯定大有裨益。最后非常感謝兩位老師和助教的辛勤付出，感謝同學(xué)在我完成大作業(yè)的過程中給予的幫助。6.附錄：MATLAB源代碼將所有.m文件復(fù)制在matlab工作目錄下運(yùn)行platform.m即可開始計(jì)算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