• <tr id="yyy80"></tr>
  • <sup id="yyy80"></sup>
  • <tfoot id="yyy80"><noscript id="yyy80"></noscript></tfoot>
  • 99热精品在线国产_美女午夜性视频免费_国产精品国产高清国产av_av欧美777_自拍偷自拍亚洲精品老妇_亚洲熟女精品中文字幕_www日本黄色视频网_国产精品野战在线观看 ?

    船舶在波浪中運(yùn)動(dòng)的強(qiáng)非線性時(shí)域模擬

    2017-04-19 09:51:05何廣華陳麗敏王佳東
    關(guān)鍵詞:船模入射波船體

    何廣華, 陳麗敏, 王佳東

    (哈爾濱工業(yè)大學(xué)(威海) 船舶與海洋工程學(xué)院,山東 威海 264209)

    船舶在波浪中運(yùn)動(dòng)的強(qiáng)非線性時(shí)域模擬

    何廣華, 陳麗敏, 王佳東

    (哈爾濱工業(yè)大學(xué)(威海) 船舶與海洋工程學(xué)院,山東 威海 264209)

    為準(zhǔn)確預(yù)報(bào)在非線性波浪中航行船舶的大幅運(yùn)動(dòng),基于CIP方法建立了波浪中船舶強(qiáng)非線性運(yùn)動(dòng)響應(yīng)的時(shí)域分析模型.該模型基于粘性流理論,對(duì)N-S方程進(jìn)行求解,可以考慮船舶的運(yùn)動(dòng)響應(yīng)和自由面變形過(guò)程中的粘性效應(yīng)和強(qiáng)非線性效應(yīng).首先,通過(guò)模擬船舶的強(qiáng)制垂蕩運(yùn)動(dòng),對(duì)該數(shù)值模型的收斂性展開(kāi)了研究;然后,基于該模型計(jì)算了不同入射波幅和不同入射波長(zhǎng)對(duì)S175船模的垂蕩和縱搖運(yùn)動(dòng)響應(yīng)的影響;最后,模擬了S175船模在大幅波浪中的強(qiáng)非線性運(yùn)動(dòng)響應(yīng)現(xiàn)象.模擬結(jié)果表明:基于自編程建立的船舶強(qiáng)非線性運(yùn)動(dòng)響應(yīng)模型具有較好的收斂性和精確性,與其他數(shù)值計(jì)算結(jié)果吻合較好;尤其擅長(zhǎng)模擬水花飛濺、甲板上浪、底部砰擊等強(qiáng)非線性問(wèn)題.

    CIP;強(qiáng)非線性;運(yùn)動(dòng)響應(yīng);大幅度運(yùn)動(dòng);時(shí)域模擬

    船舶耐波性的準(zhǔn)確預(yù)報(bào)是船舶設(shè)計(jì)的關(guān)鍵要素之一.船舶在波浪中航行時(shí),其產(chǎn)生的垂蕩和縱搖運(yùn)動(dòng)響應(yīng)會(huì)明顯增加船舶阻力;船舶的大幅運(yùn)動(dòng)會(huì)造成螺旋槳出水,推進(jìn)效率降低;同時(shí)還會(huì)影響船上設(shè)備的正常運(yùn)行,嚴(yán)重時(shí)導(dǎo)致設(shè)備的損壞,甚至船舶會(huì)因此短暫或徹底失去控制,發(fā)生傾覆.關(guān)于船舶耐波性的研究,可分為物理模型試驗(yàn)和數(shù)值模擬研究.物理模型試驗(yàn)最大的優(yōu)點(diǎn)是可靠性高,但其具有成本高、周期長(zhǎng)、難以進(jìn)行系統(tǒng)性研究等缺點(diǎn).此外,對(duì)于甲板上浪、砰擊、水花飛濺等強(qiáng)非線性現(xiàn)象,實(shí)驗(yàn)條件和測(cè)量精度仍不成熟.因此目前關(guān)于強(qiáng)非線性船舶運(yùn)動(dòng)響應(yīng)物理模型試驗(yàn)的報(bào)道仍較為罕見(jiàn);數(shù)值模擬具有經(jīng)濟(jì)性、快捷性、靈活性等優(yōu)點(diǎn),傳統(tǒng)的船舶耐波性計(jì)算方法有切片理論法[1-2],EUT (enhanced unified theory)法[3],Rankine源法[4-10]等.這些計(jì)算工具大都是基于勢(shì)流理論求解的,忽略了流體的粘性作用,無(wú)法考慮船體大幅運(yùn)動(dòng)、波浪破碎、水花飛濺、甲板上浪等強(qiáng)非線性現(xiàn)象,而這些強(qiáng)非線性效應(yīng)對(duì)船舶的運(yùn)動(dòng)響應(yīng)有著顯著的影響.相對(duì)于勢(shì)流方法,基于粘性流理論的CFD (computational fluid dynamics),可充分考慮船舶航行過(guò)程中的粘性效應(yīng)和強(qiáng)非線性效應(yīng).近年來(lái),得到了國(guó)內(nèi)外廣大學(xué)者的關(guān)注,并取得了若干研究成果.如:沈志榮等[11]基于開(kāi)源代碼OpenFOAM開(kāi)發(fā)的數(shù)值模型計(jì)算分析了Wigley Ⅲ 船模在迎浪中的運(yùn)動(dòng)響應(yīng).石博文等[12]基于粘性流理論建立了三維數(shù)值水池,分析了DTMB5512船模在不規(guī)則波浪中頂浪航行時(shí)的垂蕩與縱搖運(yùn)動(dòng)時(shí)間歷程曲線.

    CIP (constrained interpolation profile)方法由Takewaki等[13]于1985年首先提出;Yabe等[14-15]將該方法推廣到物理學(xué)、電磁學(xué)和激光學(xué)等學(xué)科.Hu等[16]于2004年首次將CIP方法應(yīng)用到船舶與海洋工程領(lǐng)域,用來(lái)計(jì)算極端非線性自由液面的流動(dòng).隨后,He等[17-20]將其應(yīng)用到處理潰壩問(wèn)題、波浪與建筑物相互作用、船體耐波性等海洋工程和船舶工程中的強(qiáng)非線性問(wèn)題.CIP法利用三次多項(xiàng)式進(jìn)行插值,并同時(shí)對(duì)物理量和其空間導(dǎo)數(shù)進(jìn)行求解,保證了數(shù)值求解在時(shí)間、空間上的三階精度.其具有精度高、穩(wěn)定性好、數(shù)值耗散小等優(yōu)點(diǎn).本文采用自編程方式,基于CIP法建立了船舶耐波性強(qiáng)非線性分析模型;并基于該數(shù)值模型對(duì)波浪中航行的S175船模的運(yùn)動(dòng)響應(yīng)展開(kāi)了研究,分析了不同入射波高和不同入射波長(zhǎng)對(duì)S175船模的垂蕩和縱搖運(yùn)動(dòng)的影響,計(jì)算結(jié)果與其他數(shù)值方法進(jìn)行了對(duì)比驗(yàn)證.最后,計(jì)算了S175船模在大浪中的大幅度運(yùn)動(dòng)響應(yīng),結(jié)果表明,該計(jì)算模型能較好地模擬船體大幅運(yùn)動(dòng)、甲板上浪、底部砰擊、水花飛濺等強(qiáng)非線性現(xiàn)象.

    1 理論與數(shù)值模型

    本研究采用自編程技術(shù),基于CIP方法建立船舶耐波性分析的粘性流數(shù)值模型.文獻(xiàn)[17]采用CIP數(shù)值模型對(duì)Wigley船的強(qiáng)制垂蕩運(yùn)動(dòng)和縱搖運(yùn)動(dòng)進(jìn)行了線性研究,計(jì)算結(jié)果與實(shí)驗(yàn)值吻合較好.限于篇幅,本文不再詳細(xì)、全面地?cái)⑹鲈揅IP模型的基礎(chǔ)理論和數(shù)值技術(shù);關(guān)于該部分,二維模型具體可參考文獻(xiàn)[16,20],三維模型可參考文獻(xiàn)[17].

    1.1 控制方程

    模型以連續(xù)性方程、N-S方程為基本控制方程.不考慮溫度變化,假設(shè)流體是不可壓縮的,可以得到以下方程:

    (1)

    (2)

    式中:ρ為密度;ui(i=1, 2, 3)為速度;μ為粘性系數(shù);Sij=?ui/?xj+?uj/?xi;fi為質(zhì)量力,如:重力等. 其中式(2)可用分裂法,將其可分成:

    (3)

    (4)

    (5)

    對(duì)流項(xiàng),式(3)采用CIP法計(jì)算;擴(kuò)散項(xiàng),式(4)采用中心差分法;壓力-速度耦合項(xiàng),式(5)通過(guò)引入連續(xù)性方程,可得以下泊松方程:

    (6)

    式(6)為壓力泊松方程,同時(shí)適用于固體、液體、氣體, 可采用SOR(超松弛迭代法)來(lái)迭代求解.

    1.2CIP法

    CIP法由文獻(xiàn)[13]首先提出;文獻(xiàn)[16]將其用于船舶與海洋工程領(lǐng)域的強(qiáng)非線性問(wèn)題的求解.隨后文獻(xiàn)[17-20]將其在海洋工程和船舶工程中進(jìn)一步推廣.

    圖1給出了一維CIP方法的原理圖.為降低數(shù)值耗散,CIP法通過(guò)補(bǔ)償網(wǎng)格間的內(nèi)部信息來(lái)提高數(shù)值精度,基本思想為:對(duì)于對(duì)流項(xiàng)變量f,不僅要計(jì)算輸運(yùn)方程的值,而且要計(jì)算輸運(yùn)方程的空間梯度(圖1(d)),即gi=?f/?xi,其中f可為ρ、ui和p.f的輸運(yùn)方程可以寫(xiě)成:

    (7)

    將式 (7)對(duì)空間坐標(biāo)進(jìn)行離散,可以得到gi的輸運(yùn)方程:

    (8)

    式(8)的計(jì)算可以分為對(duì)流項(xiàng)和非對(duì)流項(xiàng)兩步,對(duì)于式(7)和式(8)的對(duì)流項(xiàng)計(jì)算,可以用以下半拉格朗日方法求得(如圖2所示):

    圖1 一維CIP法原理示意

    式中,多項(xiàng)式的系數(shù)ai,bi,ci,di,可通過(guò)強(qiáng)加在xi-1和xi兩個(gè)網(wǎng)格點(diǎn)處的物理量fn的值以及其梯度值gn共4個(gè)量來(lái)確定.

    圖2 CIP法的半拉格朗日時(shí)間步進(jìn)法

    1.3 界面捕捉技術(shù)

    為實(shí)現(xiàn)自由面的精確重構(gòu),本數(shù)值模型選擇高精度THINC法[21]來(lái)重構(gòu)自由面.其基本思路是將流體的密度函數(shù)近似為雙曲正切函數(shù),然后通過(guò)半拉格朗日法積分得到流量,詳細(xì)論述請(qǐng)參考文獻(xiàn)[21].流體的密度函數(shù)的插值函數(shù)如下:

    1)當(dāng)ui+1/2≥0時(shí),

    2)當(dāng)ui+1/2≤0時(shí),

    式中:F(x) 為插值函數(shù);u為流速;g為流量; 未知參數(shù)α、β、γ、δ可通過(guò)等量關(guān)系和經(jīng)驗(yàn)公式得到.圖3給出了THINC法的示意圖(ui+1/2≤0).

    圖3 THINC法插值函數(shù)示意

    2 數(shù)值模擬與結(jié)果分析

    在計(jì)算波浪中船舶的運(yùn)動(dòng)響應(yīng)時(shí),船體一般可看作6自由度運(yùn)動(dòng)的剛體,且需要定義兩個(gè)坐標(biāo)系統(tǒng):一個(gè)是大地坐標(biāo)系,另一個(gè)是隨體坐標(biāo)系.本模型中采用的計(jì)算域隨著船體航速移動(dòng),船體的運(yùn)動(dòng)響應(yīng)也在該隨體坐標(biāo)系下求解.該坐標(biāo)系采用右手笛卡爾坐標(biāo),如圖4所示. 原點(diǎn)置于在船中與靜水面的交點(diǎn),x-方向以船尾指向船首為正,z-方向以豎直向上為正.

    2.1 船體模型

    本文采用ITTC-2010[22](international towing tank conference)算例對(duì)比中的標(biāo)準(zhǔn)船型——S175.S175是一艘長(zhǎng)為175 m的集裝箱船,其具體主尺度參數(shù)見(jiàn)表1, 船體表面的重構(gòu)采用虛擬粒子法,其網(wǎng)格(粒子)模型如圖5所示. 計(jì)算域?yàn)榉叫?,如圖6所示.計(jì)算域的大小x、y、z方向分別設(shè)置為:-2.7L~5.9L, -1.1L~1.1L, -0.8L~1.3L;其中L為船長(zhǎng),已在表1中列出.采用非均勻網(wǎng)格,如:在x、y、z方向的網(wǎng)格數(shù)分別為221×81×81.為提高計(jì)算精度,在靠近船首附近采用精細(xì)網(wǎng)格,最小的網(wǎng)格尺度達(dá)Δx=0.010L;在靠近船體和自由水面處也采用較精細(xì)網(wǎng)格,最小的網(wǎng)格尺度為Δy=Δz=0.005L.

    圖4 計(jì)算模型

    表1 S175集裝箱船的主尺度參數(shù)

    2.2 收斂性驗(yàn)證

    通過(guò)模擬S175船模的強(qiáng)制垂蕩運(yùn)動(dòng),對(duì)建立的數(shù)值模型展開(kāi)了收斂性研究. 計(jì)算域大小x、y、z方向分別設(shè)置為:-2.7L~5.9L, -1.1L~1.1L, -0.8L~1.3L. 為保證數(shù)值計(jì)算的穩(wěn)定性, 船舶從靜止?fàn)顟B(tài)逐漸加速到設(shè)定航速U. 鑒于篇幅,收斂性研究只展示了對(duì)數(shù)值模擬影響較大的關(guān)鍵參數(shù)的計(jì)算結(jié)果, 如:網(wǎng)格數(shù)量和時(shí)間步.

    船舶強(qiáng)制垂蕩運(yùn)動(dòng)可以表示為

    式中:j= 3表示垂蕩運(yùn)動(dòng)模態(tài);為垂蕩幅值;ωe為垂蕩運(yùn)動(dòng)圓頻率. 將縱蕩力Fx和垂蕩力Fz以ρωe進(jìn)行無(wú)因次化, 其中,ρ為流體密度,為排水體積; 時(shí)間的無(wú)因次化因子為Te= 2π/ωe(Te為船舶在波浪中的遭遇周期或強(qiáng)制運(yùn)動(dòng)周期). 圖7~8為S175船模在航速Fn= 0.2處以頻率為kL=14、振幅為kα=0.28進(jìn)行強(qiáng)制垂蕩運(yùn)動(dòng)的數(shù)值模擬結(jié)果.

    圖5 船體模型(粒子法)

    圖6 計(jì)算域與網(wǎng)絡(luò)劃分

    圖7給出了3種網(wǎng)格下無(wú)因次化垂蕩力的計(jì)算結(jié)果.3種網(wǎng)格數(shù)量分別約為Mesh-1 (1218181)、Mesh-2 (1818181) 和Mesh-3 (2218181). 相比于Mesh-1,Mesh-2與Mesh-3的計(jì)算結(jié)果更加接近, 為保證計(jì)算精度,以下的數(shù)值模擬采用Mesh-3 (2218181).

    圖8給出了3種不同時(shí)間步長(zhǎng)對(duì)無(wú)因次化垂蕩力的影響,時(shí)間步長(zhǎng)分別取為Te/1 000,Te/2 000和Te/4 000, 計(jì)算結(jié)果并沒(méi)有較大差異,Te/2 000更接近于Te/4 000.綜合考慮計(jì)算的時(shí)間和精度, 在以下算例中時(shí)間步長(zhǎng)采用Te/2 000.

    圖7 計(jì)算網(wǎng)格對(duì)垂蕩力的影響 (Fn =0.2,kL=14,kα=0.28)

    Fig.7 Sensitivity of hydrodynamic force to domain grids for the S175 model

    2.3 波浪中航行船舶的運(yùn)動(dòng)響應(yīng)

    經(jīng)過(guò)收斂性研究, 以下數(shù)值模擬均采用Mesh-3網(wǎng)格和Te/2 000的時(shí)間步長(zhǎng). 本文對(duì)S175船模在波浪中以航速Fn=0.2迎浪航行的工況展開(kāi)了系統(tǒng)的研究,分析了不同入射波幅和不同入射波長(zhǎng)對(duì)船舶運(yùn)動(dòng)響應(yīng)的影響.鑒于目前關(guān)于強(qiáng)非線性船舶運(yùn)動(dòng)響應(yīng)物理模型試驗(yàn)仍較為罕見(jiàn),ITTC于2010年以S175為標(biāo)準(zhǔn)船模進(jìn)行了多家數(shù)值模型計(jì)算結(jié)果的綜合比對(duì),供后繼程序的驗(yàn)證與校驗(yàn)之用.在本文中每個(gè)波況均計(jì)算了12個(gè)運(yùn)動(dòng)周期Te, 選取穩(wěn)定后的最后2個(gè)周期的垂蕩和縱搖運(yùn)動(dòng)響應(yīng)時(shí)間歷程線,與ITTC-2010[22]中的兩種數(shù)值模型的結(jié)果展開(kāi)了對(duì)比.計(jì)算采用主頻為3.40 GHz的Intel I3-3240 CPU, 單核計(jì)算占用內(nèi)存為680 M, 以Te/2 000的時(shí)間步長(zhǎng)計(jì)算12個(gè)運(yùn)動(dòng)周期Te所需計(jì)算時(shí)間約為11 h(隨輸出數(shù)據(jù)類型及數(shù)據(jù)量的大小有所變化,如:Tecplot流場(chǎng)數(shù)據(jù)等).

    圖8 時(shí)間步長(zhǎng)對(duì)垂蕩力的影響 (Fn =0.2,kL=14,kα=0.28)

    Fig.8 Sensitivity of hydrodynamic force to time step for the S175 model

    在圖9~11中, 帶實(shí)心方點(diǎn)的曲線為本數(shù)值模型的計(jì)算結(jié)果, 三角點(diǎn)為NMRI (national maritime research institute)的數(shù)值模擬結(jié)果; ‘×’ 為NTUA (national technical university of athens) 的數(shù)值模擬結(jié)果. NMRI所開(kāi)發(fā)的數(shù)值模型是基于弱非線性理論,采用二維邊界元方法(2D-BEM)進(jìn)行求解.而NTUA開(kāi)發(fā)的數(shù)值模型是基于三維邊界元方法,可以考慮線性和弱非線性問(wèn)題.圖9~11給出了S175船在入射波為λ/L=1.0中,以Fn=0.2航行時(shí)的垂蕩和縱搖運(yùn)動(dòng)的響應(yīng)時(shí)間歷程線; 圖9~11(a) 和 (b) 分別為無(wú)因次化的垂蕩運(yùn)動(dòng)和縱搖運(yùn)動(dòng)的計(jì)算結(jié)果,kA為無(wú)因次化入射波波幅.

    從圖9~11可以看出:1) 3種數(shù)值方法計(jì)算的垂蕩和縱搖運(yùn)動(dòng)響應(yīng)結(jié)果總體吻合較好.2) 從圖9(a)~11(a)可以看出:入射波幅對(duì)船體垂蕩運(yùn)動(dòng)的影響.當(dāng)入射波幅較小時(shí) (kA=0.04, 圖9 (a)),船體運(yùn)動(dòng)的非線性現(xiàn)象不是太明顯,仍可以按線性理論進(jìn)行考慮; NTUA的線性結(jié)果與NMRI的結(jié)果吻合較好.隨著入射波幅增加到kA=0.08和kA=0.12,船體運(yùn)動(dòng)和自由面變形的非線性現(xiàn)象更加突出,此時(shí)必須采用弱非線性或強(qiáng)非線性理論才能更好地反映和模擬真實(shí)的結(jié)果.從圖10 (a) 和圖11 (a) 可以看出:在kA=0.08和kA=0.12時(shí),本文的強(qiáng)非線性模型的計(jì)算結(jié)果與非線性理論的NMRI的結(jié)果更趨于一致.而NTUA的結(jié)果與其他兩非線性理論(NMRI和本數(shù)值方法)的計(jì)算結(jié)果相差偏大.3) 無(wú)論是垂蕩還是縱搖運(yùn)動(dòng),總的來(lái)說(shuō):本計(jì)算方法所得的運(yùn)動(dòng)響應(yīng)幅值要略小于其他兩種數(shù)值方法.其原因?yàn)?,即使在kA=0.04時(shí),本模型的計(jì)算結(jié)果也能看到波浪破碎等非線性現(xiàn)象,因此能量也會(huì)因?yàn)椴ɡ似扑榈榷兴鶞p小,從而造成船體運(yùn)動(dòng)幅值的減小.而基于邊界元方法 (NMRI和NTUA) 是無(wú)法模擬波浪破碎等強(qiáng)非線性現(xiàn)象,也不存在因?yàn)閺?qiáng)非線性因素造成的能量損失.因此船舶運(yùn)動(dòng)響應(yīng)的計(jì)算結(jié)果會(huì)較實(shí)際值偏大.

    圖9 船體的垂蕩和縱搖運(yùn)動(dòng) (Fn =0.2,λ/L=1.0, kA=0.04)

    圖10 船體的垂蕩和縱搖運(yùn)動(dòng) (Fn=0.2,λ/L=1.0,kA=0.08)

    圖11 船體的垂蕩和縱搖運(yùn)動(dòng) (Fn=0.2,λ/L=1.0,kA=0.12)

    圖12給出了S175船模在不同入射波長(zhǎng)下,以航速Fn=0.2航行時(shí)的運(yùn)動(dòng)響應(yīng)時(shí)間歷程線.眾所周知,當(dāng)船體運(yùn)動(dòng)的固有周期(如:縱搖、橫搖等周期)與波浪周期產(chǎn)生共振時(shí),船體將會(huì)激起較大的運(yùn)動(dòng)幅值.不同的運(yùn)動(dòng)模態(tài)(如:垂蕩、縱搖、橫搖),其各自的固有周期是不同的.從圖12可以看出:隨著入射波長(zhǎng)的增長(zhǎng),船體垂蕩運(yùn)動(dòng)在λ/L=1.2時(shí)產(chǎn)生較大的運(yùn)動(dòng)響應(yīng)振幅,而縱搖運(yùn)動(dòng)則在λ/L=1.4產(chǎn)生較大的運(yùn)動(dòng)響應(yīng).

    圖12 3種不同入射波長(zhǎng)下,船體的垂蕩和縱搖運(yùn)動(dòng) (Fn=0.2,kA=0.04)

    2.4 船舶強(qiáng)非線性運(yùn)動(dòng)與波形分析

    圖13給出了S175船在λ/L=1.2,kA=0.16的波浪中,以航速Fn=0.2迎浪航行時(shí)的4個(gè)不同時(shí)間的運(yùn)動(dòng)截圖.圖14給出了無(wú)因次化的船體運(yùn)動(dòng)及其所受波浪載荷隨時(shí)間的變化歷程;其無(wú)因次化形式與上述結(jié)果一致,即:縱蕩運(yùn)動(dòng)為Surge/A, 垂蕩運(yùn)動(dòng) Heave/A, 縱搖運(yùn)動(dòng)Pitch/A;縱蕩力為Fx/ρAωe2, 垂蕩力Fz/ρAωe2, 縱搖力矩My/ρALωe2. 圖13 (a) 為船體處于接近平衡位置;圖13 (b) 為船首正俯沖下水,激起較大浪花,船首部受到較大的沖擊載荷;圖13 (c) 為船體發(fā)生甲板上浪,甲板受到上浪拍擊的時(shí)刻;圖13 (d) 為船首出水后繼續(xù)前行,水花飛濺和波浪破碎在船首附近發(fā)生.從圖13 (b)~ (d) 可以看出, 該船體運(yùn)動(dòng)過(guò)程是一個(gè)伴隨著艏部沖擊、甲板上浪和水花飛濺等復(fù)雜現(xiàn)象的強(qiáng)非線性過(guò)程.因此,相應(yīng)的船體運(yùn)動(dòng)響應(yīng)和所受的波浪載荷也將是一個(gè)強(qiáng)非線性過(guò)程,這一點(diǎn)可以從圖14中較易看出,尤其是圖14 (b) 中,由于沖擊載荷的存在,波浪力曲線的局部呈現(xiàn)為鋸齒形.

    圖13 船舶強(qiáng)非線性運(yùn)動(dòng)一個(gè)周期內(nèi)的4個(gè)不同時(shí)刻的快照(Fn =0.2,λ/L=1.2, kA=0.16)

    3 結(jié) 論

    1)本文建立的波浪中航行船舶的強(qiáng)非線性運(yùn)動(dòng)響應(yīng)的數(shù)值模型,具有較好的收斂性和精確性.

    2)本模型給出的船模垂蕩和縱搖運(yùn)動(dòng)響應(yīng),與其他數(shù)值計(jì)算結(jié)果吻合較好;由于本模型考慮了粘性和強(qiáng)非線性效應(yīng),隨著非線性現(xiàn)象的增強(qiáng),本模型的垂蕩運(yùn)動(dòng)計(jì)算結(jié)果與非線性理論的NMRI的結(jié)果更趨于一致.

    3)對(duì)S175船在大幅波浪中航行時(shí)的強(qiáng)非線性運(yùn)動(dòng)響應(yīng)和自由液面大變形(或破波)的成功模擬,表明本模型具有較好地模擬水花飛濺、甲板上浪、底部砰擊等強(qiáng)非線性問(wèn)題.

    圖14 船舶強(qiáng)非線性運(yùn)動(dòng)響應(yīng)和船舶所受波浪力的時(shí)間歷程線 (Fn=0.2,λ/L=1.2,kA=0.16)

    [1] 戴仰山,沈進(jìn)威,宋競(jìng)正. 船舶波浪載荷[M]. 北京: 國(guó)防工業(yè)出版社,2007.

    [2] 陳良偉,肖桃云. 艦船波浪中非線性運(yùn)動(dòng)響應(yīng)與砰擊研究[J]. 船海工程, 2011, 40(6): 31-34. DOI: 10.3963/J.ISSN.1671-7953.2011.06.009. CHEN Liangwei, XIAO Taoyun. Study on nonlinear motion response and slam of ship in waves[J]. Ship & Ocean Engineering, 2011, 40(6): 31-34. DOI: 10.3963/J.ISSN.1671-7953.2011.06.009.

    [3] KASHIWAGI M. Prediction of surge and its effect on added resistance by means of the enhanced unified theory [J]. Transactions of the West-Japan Society of Naval Architects, 1995, 89: 77-89.

    [4] HE Guanghua, KASHIWAGI M. A time-domain higher-order boundary element method for 3D forward-speed radiation and diffraction problems [J]. Journal of Marine Science and Technology, 2014, 19(2): 228-244. DOI: 10.1007/s00773-013-0242-1.

    [5] 陳京普, 朱德祥. 船舶在波浪中運(yùn)動(dòng)的非線性時(shí)域數(shù)值模擬[J]. 水動(dòng)力學(xué)研究與進(jìn)展(A輯), 2010, 25(6): 830-836.DOI: 10.3969/j.issn.1000-4874.2010.06.015. CHEN Jingpu, ZHU Dexiang. Numerical simulations of nonlinear ship motions in waves by a Rankine panel method[J]. Journal of Hydrodynamics, 2010, 25(6): 830-836.DOI: 10.3969/j.issn.1000-4874.2010.06.015.

    [6] LIN W M, YUE D K P. Numerical solution for large amplitude ship motions in the time domain [C]// Proceedings of 18th Symposium on Naval Hydrodynamics. Washington DC: National Academy Press, 1991: 41-66.

    [7] NAKOS D, SCLAVOUNOS P. Ship motions by a three-dimensional Rankine panel method [C]// Proceedings of 18th Symposium on Naval Hydrodynamics. Washington DC: National Academy Press, 1991: 21-39.

    [8] KRING D , SCLAVOUNOS P. Numerical stability analysis for time-domain ship motion simulations [J]. Journal of Ship Research, 1995, 39(4): 313-320.

    [9] HUANG Yifeng. Nonlinear ship motions by a Rankine panel method [D]. Massachusetts: Massachusetts Institute of Technology, 1997.

    [10] KIM Y, KIM K H, KIM J H, et al. Time domain analysis of nonlinear motion responses and structural loads on ships and offshore structures: development of wish programs [J]. International Journal of Naval Architecture and Ocean Engineering, 2011, 3(1): 37-52.DOI: 10.2478/IJNAOE-2013-0044.

    [11] 沈志榮, 葉海軒, 萬(wàn)德成. 船舶在迎浪中運(yùn)動(dòng)響應(yīng)和波浪增阻的RANS數(shù)值模擬[J]. 水動(dòng)力學(xué)研究與進(jìn)展(A輯), 2012, 27(6): 621-633.DOI: 10.3969/j.issn1000-4874.2012.06.001. SHEN Zhirong, YE Haixuan, WAN Decheng. Motion response and added resistance of ship in head waves based on RANS simulations[J]. Journal of Hydrodynamics, 2012, 27(6): 621-633.DOI: 10.3969/j.issn1000-4874.2012.06.001.

    [12] 石博文, 劉正江, 吳明. 船模不規(guī)則波中頂浪運(yùn)動(dòng)數(shù)值模擬研究[J]. 船舶力學(xué), 2014, 18(8): 906-915. DOI: 10.3969/j.issn.1007-7294.2014.08.005. SHI Bowen, LIU Zhengjiang, WU Ming. Numerical simulation of ship motions in irregular head waves[J]. Journal of Ship Mechanics, 2014, 18(8): 906-915. DOI: 10.3969/j.issn.1007-7294.2014.08.005.

    [13] TAKEWAKI H, NISHIGUTI A, YABE T. Cubic interpolated pseudo-particle method (CIP) for solving hyperbolic-type equations [J]. Journal of Computational Physics, 1985, 61(2): 261-268. DOI:10.1016/0021-9991(85)90085-3.

    [14] YABE T, AOKI T. A universal solver for hyperbolic equations by cubic-polynomial interpolation I. One-dimensional solver [J]. Computer Physics Communications, 1991, 66(2/3): 219-232. DOI: 10.1016/0010-4655(91)90071-R.

    [15] YABE T, ISHIKAWA T, WANG Peiyuan, et al. A universal solver for hyperbolic equations by cubic polynomial interpolation II. Two- and three-dimensional solvers [J]. Computer Physics Communications, 1991, 66(2/3): 233-242. DOI:10.1016/0010-4655(91)90072-S.

    [16] HU Changhong, KASHIWAGI M. A CIP-based method for numerical simulations of violent free-surface flows [J]. Journal of Marine Science and Technology, 2004, 9(4): 143-157. DOI: 10.1007/s00773-004-0180-z.

    [17] HE Guanghua, ISSHIKI T, KASHIWAGI M, et al. Prediction of Radiation Forces by Means of a CIP-based Cartesian Grid Method [C]// Proceedings of 21st International Offshore and Polar Engineering Conference. Maui, Hawaii: International Society of Offshore and Polar Engineers, 2011, 2: 633-638.

    [18] ZHAO Xizeng, HU Changhong. Numerical and experimental study on a 2-D floating body under extreme wave conditions [J]. Applied Ocean Research, 2012, 35: 1-13. DOI: 10.1016/j.apor.2012.01.001.

    [19] LIAO Kangping, HU Changhong. A coupled FDMFEM method for free surface flow interaction with thin elastic plate [J]. Journal of Marine Science and Technology, 2013, 18(1): 1-11. DOI: 10.1007/s00773-012-0191-0.

    [20] HE Guanghua. A new adaptive Cartesian-grid CIP method for computation of violent free-surface flows [J]. Applied Ocean Research, 2013, 43: 234-243. DOI: 10.1016/j.apor.2013.10.004.

    [21] XIAO Feng, HONMA Y, KONO T. A simple algebraic interface capturing scheme using hyperbolic tangent function [J]. International Journal for Numerical Methods in Fluids, 2005, 48(9): 1023-1040. DOI: 10.1002/fld.975.

    [22] KIM Y. Comparative study on linear and nonlinear ship motion and loads [C]// Proceedings of ITTC Workshop on Seakeeping. Seoul: Seoul National University, 2010: 283-310.

    (編輯 張 紅)

    Strongly-nonlinear simulation of ship motions in head waves

    HE Guanghua, CHEN Limin, WANG Jiadong

    (School of Naval Architecture and Ocean Engineering, Harbin Institute of Technology in Weihai, Weihai 264209, Shandong, China)

    To simulate the ship advancing in large waves with large-amplitude motion, a CIP-based strongly-nonlinear seakeeping analysis tool is developed. The effects of viscosity and strongly-nonlinearity are considered in the numerical model, which is based on the viscous flow theory. First, convergence studies are carried out by calculating the radiation problems of S175 hull with forward speed. Then, the ship response of S175 hull including the heave and pitch motions in waves with amplitude and length of waves varied is calculated. Finally, strongly-nonlinear ship responses in large head waves are calculated and captured successfully. It is confirmed that the developed seakeeping analysis tool has numerical stability and accuracy, and has the capability of simulating the strongly-nonlinear seakeeping problem, such as, the phenomena of wave breaking, slamming, and green water on deck.

    CIP; strongly-nonlinear; ship motion response; large amplitude motion; time domain simulation

    10.11918/j.issn.0367-6234.201508063

    2015-08-21

    國(guó)家自然科學(xué)基金 (51579058); 山東省自然科學(xué)基金 (ZR2014EEQ016)

    何廣華(1980—),男,教授,博士生導(dǎo)師

    何廣華,ghhe@hitwh.edu.cn

    U661.31

    A

    0367-6234(2017)04-0142-07

    猜你喜歡
    船模入射波船體
    基于模糊PID的船模航向控制研究
    船體行駛過(guò)程中的壓力監(jiān)測(cè)方法
    SHPB入射波相似律與整形技術(shù)的試驗(yàn)與數(shù)值研究
    瞬態(tài)激勵(lì)狀態(tài)下樁身速度以及樁身內(nèi)力計(jì)算
    船模靜水橫搖試驗(yàn)的不確定度分析
    焊接殘余應(yīng)力對(duì)船體結(jié)構(gòu)疲勞強(qiáng)度的影響分析
    焊接(2015年9期)2015-07-18 11:03:51
    對(duì)機(jī)械波半波損失現(xiàn)象的物理解釋
    電子科技(2015年11期)2015-03-06 01:32:24
    赴美軍“仁慈”號(hào)醫(yī)院船駐船體會(huì)
    西洋船模王——童鑑良
    航海(2014年6期)2014-12-12 10:36:03
    水下爆炸氣泡作用下船體總縱強(qiáng)度估算方法
    十八禁网站免费在线| 国产伦人伦偷精品视频| 动漫黄色视频在线观看| 美女福利国产在线| 午夜激情久久久久久久| 国产一区二区三区综合在线观看| 国产精品麻豆人妻色哟哟久久| 成年av动漫网址| 亚洲 欧美一区二区三区| 两人在一起打扑克的视频| 90打野战视频偷拍视频| 欧美日韩中文字幕国产精品一区二区三区 | 亚洲精品美女久久久久99蜜臀| 18禁观看日本| 精品少妇一区二区三区视频日本电影| 国产成人av激情在线播放| 亚洲综合色网址| 亚洲天堂av无毛| 免费在线观看黄色视频的| 午夜福利影视在线免费观看| 国产男女内射视频| 亚洲人成电影免费在线| 久久天躁狠狠躁夜夜2o2o| 国产男女超爽视频在线观看| 宅男免费午夜| 精品亚洲成a人片在线观看| 久久久久国产精品人妻一区二区| 亚洲欧美日韩高清在线视频 | 欧美日韩中文字幕国产精品一区二区三区 | 中文字幕精品免费在线观看视频| 久久99热这里只频精品6学生| 日韩一卡2卡3卡4卡2021年| 久久久国产精品麻豆| 两个人免费观看高清视频| 国产高清视频在线播放一区 | 中文字幕另类日韩欧美亚洲嫩草| 午夜福利一区二区在线看| 新久久久久国产一级毛片| 侵犯人妻中文字幕一二三四区| 亚洲自偷自拍图片 自拍| 欧美大码av| 十分钟在线观看高清视频www| www.自偷自拍.com| 亚洲精品国产区一区二| 国产亚洲精品一区二区www | 亚洲精品国产av蜜桃| 午夜视频精品福利| 黑人巨大精品欧美一区二区mp4| 99久久99久久久精品蜜桃| 国产成人av激情在线播放| 国产亚洲欧美精品永久| e午夜精品久久久久久久| 在线观看免费午夜福利视频| 亚洲成国产人片在线观看| 国产精品亚洲av一区麻豆| 国产精品成人在线| 国产精品 国内视频| 欧美亚洲 丝袜 人妻 在线| 亚洲激情五月婷婷啪啪| 国产成人精品久久二区二区91| 久久精品人人爽人人爽视色| 久久精品aⅴ一区二区三区四区| a 毛片基地| 老司机影院成人| 色播在线永久视频| 精品高清国产在线一区| 久久国产精品人妻蜜桃| 亚洲精品av麻豆狂野| 男女国产视频网站| 国产亚洲午夜精品一区二区久久| 一级毛片女人18水好多| 50天的宝宝边吃奶边哭怎么回事| 国产男女内射视频| 99久久精品国产亚洲精品| 亚洲成国产人片在线观看| 日韩中文字幕视频在线看片| 18禁裸乳无遮挡动漫免费视频| 大码成人一级视频| 一级片免费观看大全| 肉色欧美久久久久久久蜜桃| 亚洲av日韩精品久久久久久密| 美女高潮到喷水免费观看| 色94色欧美一区二区| 亚洲少妇的诱惑av| 一区二区日韩欧美中文字幕| 亚洲熟女精品中文字幕| 人成视频在线观看免费观看| 大香蕉久久网| 亚洲成人国产一区在线观看| 久久久国产欧美日韩av| 在线观看www视频免费| 老熟女久久久| 91字幕亚洲| 纯流量卡能插随身wifi吗| 中文字幕制服av| 一级a爱视频在线免费观看| 不卡av一区二区三区| 精品高清国产在线一区| 无限看片的www在线观看| 亚洲专区字幕在线| 亚洲一区中文字幕在线| 日韩一区二区三区影片| 久久久国产欧美日韩av| 久久精品熟女亚洲av麻豆精品| 欧美在线黄色| 国产免费av片在线观看野外av| 一二三四社区在线视频社区8| kizo精华| 欧美黄色片欧美黄色片| 国产一卡二卡三卡精品| 老汉色av国产亚洲站长工具| 大型av网站在线播放| a在线观看视频网站| 精品久久久久久电影网| 亚洲国产中文字幕在线视频| 91精品国产国语对白视频| av不卡在线播放| 天堂俺去俺来也www色官网| 国产91精品成人一区二区三区 | 亚洲国产av新网站| 色综合欧美亚洲国产小说| 多毛熟女@视频| 黑人猛操日本美女一级片| 人人澡人人妻人| 男男h啪啪无遮挡| 天天躁日日躁夜夜躁夜夜| 欧美日韩国产mv在线观看视频| 免费日韩欧美在线观看| 久久久久久免费高清国产稀缺| 妹子高潮喷水视频| 久久精品亚洲熟妇少妇任你| 777久久人妻少妇嫩草av网站| 国产无遮挡羞羞视频在线观看| 另类精品久久| 中文字幕制服av| 好男人电影高清在线观看| 伊人久久大香线蕉亚洲五| 国产又爽黄色视频| 一个人免费看片子| videos熟女内射| 国产视频一区二区在线看| 又紧又爽又黄一区二区| 欧美成人午夜精品| 蜜桃国产av成人99| 嫩草影视91久久| 精品亚洲成国产av| 欧美变态另类bdsm刘玥| 19禁男女啪啪无遮挡网站| 在线精品无人区一区二区三| 久久国产精品大桥未久av| 亚洲激情五月婷婷啪啪| 精品熟女少妇八av免费久了| 久久久国产一区二区| 色播在线永久视频| netflix在线观看网站| 热99国产精品久久久久久7| 啦啦啦啦在线视频资源| 999精品在线视频| 日韩欧美国产一区二区入口| 老司机靠b影院| tube8黄色片| 一本一本久久a久久精品综合妖精| 视频区欧美日本亚洲| 亚洲一码二码三码区别大吗| 中文欧美无线码| 午夜福利在线免费观看网站| 国产成人精品在线电影| 永久免费av网站大全| 午夜福利在线观看吧| 夜夜夜夜夜久久久久| 午夜日韩欧美国产| 午夜福利影视在线免费观看| 麻豆乱淫一区二区| 精品人妻1区二区| 国产又色又爽无遮挡免| 国产精品.久久久| 中亚洲国语对白在线视频| av福利片在线| 91精品三级在线观看| 久久精品aⅴ一区二区三区四区| 在线精品无人区一区二区三| a级毛片在线看网站| 三级毛片av免费| 精品一区在线观看国产| 王馨瑶露胸无遮挡在线观看| a 毛片基地| 老熟女久久久| 人人妻人人爽人人添夜夜欢视频| 国产欧美亚洲国产| 日本a在线网址| 亚洲国产精品一区二区三区在线| 精品一区二区三区四区五区乱码| 亚洲精品国产一区二区精华液| 国产精品 欧美亚洲| h视频一区二区三区| 一二三四在线观看免费中文在| 蜜桃国产av成人99| 欧美av亚洲av综合av国产av| 亚洲精品久久成人aⅴ小说| 女人久久www免费人成看片| 永久免费av网站大全| 国产av又大| 女人精品久久久久毛片| 国产高清国产精品国产三级| 一区二区三区四区激情视频| 国产亚洲精品一区二区www | 午夜免费观看性视频| 欧美日本中文国产一区发布| 天天操日日干夜夜撸| 男人操女人黄网站| 国产男人的电影天堂91| 免费在线观看日本一区| 国产精品久久久久成人av| 超碰97精品在线观看| 国产亚洲精品第一综合不卡| 男人添女人高潮全过程视频| 在线精品无人区一区二区三| 肉色欧美久久久久久久蜜桃| 后天国语完整版免费观看| 亚洲国产av新网站| 国精品久久久久久国模美| 欧美乱码精品一区二区三区| 丰满饥渴人妻一区二区三| 日韩 亚洲 欧美在线| 可以免费在线观看a视频的电影网站| av又黄又爽大尺度在线免费看| 免费一级毛片在线播放高清视频 | 国产三级黄色录像| 亚洲激情五月婷婷啪啪| 日韩中文字幕视频在线看片| 啦啦啦在线免费观看视频4| 一级a爱视频在线免费观看| 丝袜喷水一区| 亚洲国产精品999| 黄频高清免费视频| 又大又爽又粗| 午夜久久久在线观看| 色播在线永久视频| 老汉色∧v一级毛片| 国产精品麻豆人妻色哟哟久久| 亚洲成av片中文字幕在线观看| 久久99一区二区三区| 国产国语露脸激情在线看| 男女国产视频网站| 老司机亚洲免费影院| 亚洲国产看品久久| 亚洲av电影在线进入| 精品一区二区三区四区五区乱码| 亚洲自偷自拍图片 自拍| 国产精品麻豆人妻色哟哟久久| 亚洲av日韩精品久久久久久密| 免费人妻精品一区二区三区视频| 精品第一国产精品| 亚洲专区字幕在线| 伊人久久大香线蕉亚洲五| 他把我摸到了高潮在线观看 | 人妻 亚洲 视频| 人成视频在线观看免费观看| 性少妇av在线| 人人妻人人添人人爽欧美一区卜| 欧美性长视频在线观看| 午夜视频精品福利| 日韩三级视频一区二区三区| 亚洲国产欧美一区二区综合| 成人国语在线视频| 免费在线观看日本一区| 蜜桃国产av成人99| www.av在线官网国产| 免费高清在线观看视频在线观看| 国产高清videossex| 色播在线永久视频| 97在线人人人人妻| 岛国毛片在线播放| 亚洲美女黄色视频免费看| svipshipincom国产片| 国产精品久久久久成人av| 91国产中文字幕| 别揉我奶头~嗯~啊~动态视频 | 成年av动漫网址| 亚洲男人天堂网一区| 最近最新免费中文字幕在线| 可以免费在线观看a视频的电影网站| 亚洲国产欧美网| 建设人人有责人人尽责人人享有的| 精品乱码久久久久久99久播| 黄频高清免费视频| 99热国产这里只有精品6| svipshipincom国产片| 国产成+人综合+亚洲专区| 精品国产超薄肉色丝袜足j| 精品人妻在线不人妻| 香蕉丝袜av| 99久久国产精品久久久| 午夜免费鲁丝| 亚洲成av片中文字幕在线观看| 色婷婷久久久亚洲欧美| 十八禁高潮呻吟视频| 在线精品无人区一区二区三| 日本av免费视频播放| 啦啦啦在线免费观看视频4| 激情视频va一区二区三区| 精品乱码久久久久久99久播| 一区二区三区精品91| 国产一区二区 视频在线| 老司机亚洲免费影院| 精品人妻1区二区| 欧美日本中文国产一区发布| 嫁个100分男人电影在线观看| 老司机福利观看| netflix在线观看网站| 欧美日韩黄片免| 考比视频在线观看| 精品国产乱码久久久久久男人| 欧美另类一区| 老熟妇仑乱视频hdxx| 天堂中文最新版在线下载| 精品一区二区三区av网在线观看 | 18禁裸乳无遮挡动漫免费视频| 久久精品国产综合久久久| av在线老鸭窝| 91av网站免费观看| 国产精品免费大片| 999精品在线视频| 亚洲性夜色夜夜综合| 久久久久国产一级毛片高清牌| 在线 av 中文字幕| 两个人免费观看高清视频| 亚洲av国产av综合av卡| 久久综合国产亚洲精品| 激情视频va一区二区三区| 精品人妻熟女毛片av久久网站| bbb黄色大片| a级毛片黄视频| 国产在线视频一区二区| 亚洲七黄色美女视频| 桃花免费在线播放| 免费在线观看日本一区| 在线十欧美十亚洲十日本专区| 亚洲 国产 在线| 男女无遮挡免费网站观看| 亚洲精品美女久久av网站| 亚洲av电影在线进入| 国产真人三级小视频在线观看| 亚洲av欧美aⅴ国产| 麻豆国产av国片精品| 久久久久久久精品精品| 女性被躁到高潮视频| 亚洲少妇的诱惑av| 亚洲国产日韩一区二区| 69av精品久久久久久 | 亚洲精品久久久久久婷婷小说| 亚洲第一欧美日韩一区二区三区 | 69av精品久久久久久 | 超色免费av| 午夜激情久久久久久久| 黄色毛片三级朝国网站| 亚洲一区中文字幕在线| 97在线人人人人妻| 十八禁网站网址无遮挡| 伊人亚洲综合成人网| 天堂8中文在线网| 国产深夜福利视频在线观看| 亚洲九九香蕉| 国产97色在线日韩免费| 男人爽女人下面视频在线观看| 亚洲视频免费观看视频| 中文字幕人妻熟女乱码| 一本大道久久a久久精品| 国产av一区二区精品久久| 欧美日韩福利视频一区二区| 丰满少妇做爰视频| 午夜影院在线不卡| 汤姆久久久久久久影院中文字幕| 成人18禁高潮啪啪吃奶动态图| 悠悠久久av| 久久精品aⅴ一区二区三区四区| 男女下面插进去视频免费观看| 亚洲国产av影院在线观看| 亚洲色图综合在线观看| 在线观看人妻少妇| 十八禁人妻一区二区| 亚洲精品av麻豆狂野| 国产又爽黄色视频| 久久香蕉激情| a 毛片基地| 免费高清在线观看日韩| www.熟女人妻精品国产| 黄色 视频免费看| 精品国产乱码久久久久久男人| 91大片在线观看| 国产日韩欧美亚洲二区| a在线观看视频网站| 我要看黄色一级片免费的| 大香蕉久久网| 欧美另类一区| 精品亚洲乱码少妇综合久久| 一级,二级,三级黄色视频| 国产精品成人在线| 久久精品国产综合久久久| 久久毛片免费看一区二区三区| 热99re8久久精品国产| 97精品久久久久久久久久精品| 亚洲人成77777在线视频| www.熟女人妻精品国产| 国产在线一区二区三区精| 又大又爽又粗| 热99国产精品久久久久久7| 久久人妻熟女aⅴ| 久久狼人影院| 亚洲免费av在线视频| 国产亚洲av高清不卡| 久久国产精品大桥未久av| www.熟女人妻精品国产| 亚洲av成人不卡在线观看播放网 | 亚洲精品久久午夜乱码| 天天躁狠狠躁夜夜躁狠狠躁| 久久性视频一级片| 国产成人a∨麻豆精品| 精品国产一区二区三区四区第35| 亚洲一码二码三码区别大吗| 中文字幕最新亚洲高清| 在线观看人妻少妇| 色精品久久人妻99蜜桃| 久久久久久久精品精品| 无遮挡黄片免费观看| av超薄肉色丝袜交足视频| 黄片播放在线免费| 韩国精品一区二区三区| 欧美xxⅹ黑人| 精品一区二区三区av网在线观看 | 精品一区二区三区av网在线观看 | 超碰成人久久| 久久久精品免费免费高清| 精品国产乱码久久久久久小说| 日韩 亚洲 欧美在线| 国产精品.久久久| 午夜久久久在线观看| 亚洲精品在线美女| 亚洲中文av在线| 91麻豆av在线| 午夜精品国产一区二区电影| 麻豆av在线久日| 日日爽夜夜爽网站| 巨乳人妻的诱惑在线观看| 黄片大片在线免费观看| 99九九在线精品视频| 一本一本久久a久久精品综合妖精| 一本—道久久a久久精品蜜桃钙片| 欧美日韩av久久| 黄色视频,在线免费观看| 法律面前人人平等表现在哪些方面 | 一区福利在线观看| 人妻久久中文字幕网| 一区二区三区激情视频| 肉色欧美久久久久久久蜜桃| 国产精品一区二区免费欧美 | 国产av精品麻豆| 国产男人的电影天堂91| 性色av一级| 久久精品国产亚洲av高清一级| 一区二区三区四区激情视频| 日韩视频在线欧美| 国产极品粉嫩免费观看在线| 韩国高清视频一区二区三区| 精品国产一区二区三区久久久樱花| 亚洲精品国产一区二区精华液| 成年av动漫网址| 丝袜人妻中文字幕| 精品国产国语对白av| 老汉色av国产亚洲站长工具| 老司机午夜福利在线观看视频 | 亚洲欧洲精品一区二区精品久久久| av天堂在线播放| 亚洲黑人精品在线| 久久精品国产综合久久久| 99re6热这里在线精品视频| 中文字幕最新亚洲高清| 国产免费av片在线观看野外av| 国产麻豆69| 国产视频一区二区在线看| 纵有疾风起免费观看全集完整版| 丁香六月欧美| 91麻豆精品激情在线观看国产 | 亚洲精品国产一区二区精华液| 亚洲精品一区蜜桃| 欧美亚洲 丝袜 人妻 在线| 我要看黄色一级片免费的| 老司机影院毛片| 黄色怎么调成土黄色| 青春草视频在线免费观看| bbb黄色大片| 国产精品国产三级国产专区5o| 交换朋友夫妻互换小说| 中文字幕高清在线视频| 人人妻人人爽人人添夜夜欢视频| 午夜久久久在线观看| 国产深夜福利视频在线观看| 欧美亚洲 丝袜 人妻 在线| 男人爽女人下面视频在线观看| 久久女婷五月综合色啪小说| 黄色怎么调成土黄色| 在线观看免费视频网站a站| 久久久国产成人免费| 亚洲午夜精品一区,二区,三区| 狠狠婷婷综合久久久久久88av| 日韩免费高清中文字幕av| 久久久国产精品麻豆| 久久综合国产亚洲精品| 国产又色又爽无遮挡免| svipshipincom国产片| 丁香六月天网| 一进一出抽搐动态| 999久久久精品免费观看国产| 十八禁人妻一区二区| 九色亚洲精品在线播放| 国产伦人伦偷精品视频| 国精品久久久久久国模美| 亚洲第一青青草原| 亚洲avbb在线观看| 午夜影院在线不卡| 99精品久久久久人妻精品| 国产日韩欧美视频二区| 99国产精品一区二区蜜桃av | 久久人妻福利社区极品人妻图片| 免费一级毛片在线播放高清视频 | av电影中文网址| 国产一区二区激情短视频 | 午夜免费观看性视频| 一级片'在线观看视频| 欧美av亚洲av综合av国产av| 国产成人av教育| 狠狠精品人妻久久久久久综合| 欧美人与性动交α欧美精品济南到| 欧美 日韩 精品 国产| 国产国语露脸激情在线看| 久久久久久久国产电影| 男人操女人黄网站| 另类精品久久| 国产色视频综合| 天堂俺去俺来也www色官网| 亚洲精品久久久久久婷婷小说| 麻豆乱淫一区二区| 亚洲精华国产精华精| 两性午夜刺激爽爽歪歪视频在线观看 | 久久这里只有精品19| 久久人人97超碰香蕉20202| 后天国语完整版免费观看| 丰满饥渴人妻一区二区三| 中文欧美无线码| 侵犯人妻中文字幕一二三四区| 免费高清在线观看视频在线观看| 下体分泌物呈黄色| 亚洲av欧美aⅴ国产| 777久久人妻少妇嫩草av网站| 欧美xxⅹ黑人| 两性午夜刺激爽爽歪歪视频在线观看 | 脱女人内裤的视频| 亚洲国产日韩一区二区| 免费在线观看日本一区| 99精品欧美一区二区三区四区| 正在播放国产对白刺激| 午夜福利,免费看| 久久人人爽av亚洲精品天堂| 国产免费av片在线观看野外av| 高清在线国产一区| 久久精品国产亚洲av香蕉五月 | 国产淫语在线视频| 精品久久蜜臀av无| 少妇精品久久久久久久| 婷婷丁香在线五月| 日本欧美视频一区| 中文字幕色久视频| 精品一区二区三卡| av网站在线播放免费| 欧美在线黄色| 美女扒开内裤让男人捅视频| 久久精品国产亚洲av香蕉五月 | 人妻 亚洲 视频| 国产激情久久老熟女| 国产精品国产av在线观看| 亚洲精品国产一区二区精华液| 国产一卡二卡三卡精品| 黄片播放在线免费| avwww免费| 亚洲国产av新网站| 国产在线一区二区三区精| 国产精品偷伦视频观看了| av不卡在线播放| 91av网站免费观看| 欧美激情久久久久久爽电影 | 久热爱精品视频在线9| 中亚洲国语对白在线视频| 亚洲成人免费电影在线观看| 亚洲一卡2卡3卡4卡5卡精品中文| 嫩草影视91久久| 中文字幕色久视频| 国产主播在线观看一区二区| 高清在线国产一区| 国产精品二区激情视频| 国产99久久九九免费精品| 五月开心婷婷网| 在线观看免费高清a一片| 国产一区二区三区在线臀色熟女 | 少妇人妻久久综合中文| 亚洲专区中文字幕在线| 免费在线观看黄色视频的| 三上悠亚av全集在线观看| 欧美日韩一级在线毛片| 国产欧美日韩一区二区三区在线| 别揉我奶头~嗯~啊~动态视频 | 欧美在线一区亚洲| 五月天丁香电影| 亚洲一码二码三码区别大吗| 欧美激情久久久久久爽电影 | 97在线人人人人妻|