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

    有限點(diǎn)法在液艙晃蕩問題中的應(yīng)用*

    2015-02-18 04:59:08盧雨胡安康劉亞沖

    盧雨 胡安康,2? 劉亞沖

    (1.哈爾濱工程大學(xué) 船舶工程學(xué)院, 黑龍江 哈爾濱 150001; 2.中集船舶海洋工程設(shè)計(jì)研究院,上海 201206)

    有限點(diǎn)法在液艙晃蕩問題中的應(yīng)用*

    盧雨1胡安康1,2?劉亞沖1

    (1.哈爾濱工程大學(xué) 船舶工程學(xué)院, 黑龍江 哈爾濱 150001; 2.中集船舶海洋工程設(shè)計(jì)研究院,上海 201206)

    摘要:液艙晃蕩問題是近年來液貨船設(shè)計(jì)中十分關(guān)注的問題之一,它涉及自由液面的大幅運(yùn)動(dòng)等一系列強(qiáng)非線性問題,因此給基于歐拉觀點(diǎn)下空間拓?fù)渚W(wǎng)格的數(shù)值模擬帶來了較大的難度.文中基于拉格朗日觀點(diǎn)下的有限點(diǎn)法(FPM),對(duì)不同充容率的液艙晃蕩問題進(jìn)行了數(shù)值模擬,并將數(shù)值計(jì)算結(jié)果與試驗(yàn)值進(jìn)行對(duì)比驗(yàn)證.結(jié)果表明:當(dāng)激勵(lì)頻率與液艙固有頻率一致時(shí),艙內(nèi)流體的晃蕩較為劇烈;而當(dāng)充容率不同時(shí),艙壁對(duì)液體流動(dòng)的抑制作用使得高充容液艙內(nèi)的流體翻卷變形較為平緩.通過有限點(diǎn)法有效且準(zhǔn)確地模擬了不同充容率液艙在各激勵(lì)頻率下的水面晃蕩情況,與試驗(yàn)現(xiàn)象基本吻合.同時(shí)通過對(duì)艙壁監(jiān)測(cè)點(diǎn)的壓力計(jì)算得知,有限點(diǎn)法的數(shù)值結(jié)果與實(shí)測(cè)壓力曲線基本保持一致.FPM針對(duì)不同工況下的二維液艙晃蕩情況均給出了較為準(zhǔn)確的數(shù)值模擬結(jié)果,有效驗(yàn)證了有限點(diǎn)法的可靠性.

    關(guān)鍵詞:液艙晃蕩;拉格朗日觀點(diǎn);有限點(diǎn)法;移動(dòng)最小二乘法;泊松方程

    近年來,液貨船在船舶交易市場(chǎng)中的需求不斷增長(zhǎng),顯著刺激了船舶運(yùn)輸技術(shù)的不斷提高.液貨船的關(guān)鍵技術(shù)之一就是液艙內(nèi)流體的晃蕩問題.常見的晃蕩現(xiàn)象有:駐波、行進(jìn)波、水躍和組合波,還伴有水躍、翻卷等強(qiáng)非線性現(xiàn)象.而液體晃蕩產(chǎn)生的沖擊壓力會(huì)嚴(yán)重威脅到艙壁以及船體結(jié)構(gòu)的安全性能,并在一定程度上影響航行中船舶的穩(wěn)性,強(qiáng)烈的晃蕩將直接導(dǎo)致嚴(yán)重的海洋結(jié)構(gòu)物的破壞[1-2],因此液艙晃蕩研究已在全世界范圍內(nèi)受到越來越多的關(guān)注.

    液體晃蕩是一種高度非線性的復(fù)雜流動(dòng)現(xiàn)象,這種非線性來源于自由液面的大幅度運(yùn)動(dòng)、濕邊界的變化以及流固耦合等因素,嚴(yán)重阻礙了問題的理論研究.隨著計(jì)算機(jī)技術(shù)的迅猛發(fā)展,數(shù)值模擬技術(shù)成為研究液艙晃蕩問題的重要手段.Armenio[3]采用改進(jìn) MAC (Marker-And-Cell) 方法(即標(biāo)記子與單元方法)將流體壓力和速度作為求解變量,成功地求解了帶有自由面的流體大幅晃蕩問題.Kim等[4-6]基于有限差分法數(shù)值模擬了液艙晃蕩運(yùn)動(dòng),Hu等[7]采用約束插值輪廓法(CIP)研究了劇烈液艙晃蕩現(xiàn)象.祁江濤等[8]采用兩相流 VOF (Volume of Fluid) 法,結(jié)合動(dòng)網(wǎng)格技術(shù),針對(duì)二維及三維液艙晃蕩情況進(jìn)行了數(shù)值計(jì)算.沈猛等[9-10]改進(jìn)了傳統(tǒng)VOF法,基于部分單元參數(shù)概念,引入了一種混合自由表面邊界速度條件的基本原理,解決了具有復(fù)雜邊界的液艙晃蕩問題.方智勇等[11-12]依據(jù)Level-set法隱式跟蹤運(yùn)動(dòng)界面原理,數(shù)值模擬了由箱體做簡(jiǎn)諧運(yùn)動(dòng)引起的液艙晃蕩問題,證明了此法的可行性.然而這些 CFD (Computational Fluid Dynamics) 技術(shù)在處理液艙晃蕩問題時(shí)往往是基于空間網(wǎng)格拓?fù)浣Y(jié)構(gòu),極易造成數(shù)值擴(kuò)散導(dǎo)致捕捉液面模糊等問題,并且很難模擬出液體飛濺、融合等復(fù)雜自由表面流現(xiàn)象,因此這些基于歐拉觀點(diǎn)的數(shù)值求解技術(shù)在處理大變形的流體運(yùn)動(dòng)方面遇到了一些難以克服的問題.

    在處理流體復(fù)雜變形的問題上,相比于傳統(tǒng)歐拉觀點(diǎn)的網(wǎng)格類方法,基于拉格朗日觀點(diǎn)發(fā)展而來的無網(wǎng)格粒子法受到眾多學(xué)者的青睞[13].有限點(diǎn)法Finite Point Method,FPM) 是近些年發(fā)展起來的新型Lagrange法,是一種真正意義上的無網(wǎng)格技術(shù).它首先由西班牙學(xué)者Oate等[14]于1996年提出,用有限空間粒子點(diǎn)來離散表達(dá)流體域,不受固定網(wǎng)格結(jié)構(gòu)形式的束縛,每個(gè)粒子點(diǎn)攜帶特定的物理信息,如速度壓力等,粒子間的相互作用通過移動(dòng)最小二乘法近似得到,這樣就較為容易地構(gòu)造出所期望的順序一致性,并基于配點(diǎn)法離散求解流體運(yùn)動(dòng)控制方程,達(dá)到時(shí)刻監(jiān)測(cè)粒子運(yùn)動(dòng)狀態(tài)并準(zhǔn)確追蹤粒子運(yùn)動(dòng)軌跡的效果.此外有限點(diǎn)法對(duì)于邊界的處理采用較為通用的邊界粒子分布法,易于數(shù)值實(shí)現(xiàn).FPM在數(shù)值求解與液面捕捉方面的靈活性,使其被應(yīng)用到許多流動(dòng)問題的求解.早年Oate等[15]基于該有限點(diǎn)法, 對(duì)不可壓縮流體運(yùn)動(dòng)問題進(jìn)行了數(shù)值研究.隨后Oate等[16]將FPM應(yīng)用于求解結(jié)構(gòu)力學(xué)中的彈性力學(xué)問題,并在線彈性結(jié)構(gòu)分析中取得了一定的進(jìn)展.2005年,Shu 等[17]基于無網(wǎng)格法中的離散原理及局部等參插值理論,提出了一種新型無網(wǎng)格法-等參有限點(diǎn)法(IFPM),并將用于其求解粘性不可壓縮流等問題,得到了滿意的結(jié)果.之后,Pandey 等[18]利用有限點(diǎn)法數(shù)值求解了熱力學(xué)流動(dòng)中的 LLNS (Landau-Lifshitz Navier-Stokes) 方程,并與理論值進(jìn)行有效性驗(yàn)證,結(jié)果吻合較好.2014年,Reséndiz 等[19]利用泰勒級(jí)數(shù)形式的有限點(diǎn)法對(duì)二維熱傳導(dǎo)問題進(jìn)行研究,并與解析解對(duì)比,得到了滿意的精度.由此可見,有限點(diǎn)法的研究已受到各國(guó)學(xué)者的重視.

    文中基于有限點(diǎn)法對(duì)液艙晃蕩問題進(jìn)行了數(shù)值模擬,對(duì)流體大幅運(yùn)動(dòng)特性、艙壁周期性砰擊壓力以及液面大變形進(jìn)行了分析.主要采用投影法對(duì)不可壓縮Navier-Stokes方程進(jìn)行求解,壓力泊松方程中的空間導(dǎo)數(shù)由移動(dòng)最小二乘法近似獲得.文中選用新四次函數(shù)作為計(jì)算的權(quán)函數(shù),是因?yàn)樾滤拇魏瘮?shù)是一條連續(xù)曲線,且二階導(dǎo)數(shù)曲線具有較好的光滑性,可光順連續(xù)地描述大變形自由面.最后通過數(shù)值計(jì)算結(jié)果與試驗(yàn)值的對(duì)比,對(duì)有限點(diǎn)法在求解液艙晃蕩問題中的可靠性進(jìn)行了驗(yàn)證.

    1數(shù)值方法

    1.1 控制方程

    不可壓縮流體的運(yùn)動(dòng)控制方程(質(zhì)量守恒、動(dòng)量守恒)表述為

    (1)

    (2)

    式中,ρ為流體密度, ui為坐標(biāo)系中對(duì)應(yīng)i軸的速度分量, xi為對(duì)應(yīng)i軸的方向分量,p為壓力,μ為流體的動(dòng)力粘性系數(shù),fi為源項(xiàng)(通常取重力加速度gi).由于是從拉格朗日觀點(diǎn)出發(fā),故式(1)和式(2)右端的時(shí)間導(dǎo)數(shù)項(xiàng)是以物質(zhì)導(dǎo)數(shù)的形式給出.

    1.2 投影法

    (3)

    (4)

    (5)

    式(5)中增加了壓力梯度項(xiàng),通過引入不可壓縮流體的連續(xù)方程,可顯式獲得關(guān)于壓力的泊松方程,即為

    (6)

    將式(5)沿邊界Γ的外單位法向量n投影,得到關(guān)于壓力p的Neumann 邊界條件,即

    (7)

    1.3 壓力泊松方程的導(dǎo)數(shù)離散

    在投影法求解壓力泊松方程過程中,產(chǎn)生了壓力項(xiàng)的方向?qū)?shù)及Laplacian算子.本節(jié)利用移動(dòng)最小二乘法的思想,將空間任意點(diǎn)的待求函數(shù)導(dǎo)數(shù)用支持域內(nèi)的周圍粒子信息來近似表達(dá).

    首先令流場(chǎng)中全部粒子點(diǎn)(粒子總數(shù)為n)的函數(shù)向量為U=[u1(x),u2(x),…,un(x)],u(x)為任意粒子點(diǎn)x處的待求函數(shù),則其近似函數(shù)u(x)可表示為

    (8)

    式中:b(x)T=[b1(x)b2(x)…bk(x)],bi(x)為基函數(shù),k為基函數(shù)的個(gè)數(shù);a(x)T=[a1(x)a2(x)…ak(x)],ai(x)為待定系數(shù).文中計(jì)算取二維空間單項(xiàng)式基函數(shù):

    b(x)T=[1xyx2xyy2],k=6

    (9)

    考慮其計(jì)算通用性,將基函數(shù)bi(x)寫入全局矩陣中,即

    (10)

    根據(jù)最小二乘法求取極值原理,可得如下關(guān)系:

    S(x)a(x)=r(x)

    (11)

    其中S(x)=VW(x)VT,r(x)=VW(x)U,W(x)為權(quán)函數(shù)矩陣,即

    (12)

    分別對(duì)式(11)兩邊求一階及二階導(dǎo)數(shù),可得

    a′=S-1(r′-S′S-1r)

    (13)

    a″=S-1((r″-S″S-1r)-2S′S-1(r′-S′S-1r))

    (14)

    根據(jù)式(13)可得近似函數(shù)的一階導(dǎo)數(shù):

    (∏u)′=(b′)Ta+bTa′=

    (b′)TS-1r+bTS-1(r′-S′S-1r)

    (15)

    聯(lián)立式(13)與(14)可得其近似函數(shù)二階導(dǎo)數(shù):

    (∏u)″=(b″)Ta+2(b′)Ta′+bTa″=

    (b″)TS-1r+2(b′)TS-1(r′-S′S-1r)+

    bTS-1((r″-S″S-1r)-2S′S-1(r′-S′S-1r))

    (16)

    借助式(16)及(17),即可離散壓力泊松方程,最后形成關(guān)于壓力的大型稀疏矩陣,采用較為通用的雙共軛梯度穩(wěn)定法(Bi-ConjugateGradientsStabilizedMethod)[20]來迭代求解,并設(shè)定迭代前后時(shí)間步內(nèi)的壓力誤差限作為計(jì)算收斂條件.同時(shí)為滿足Dirichlet條件,在計(jì)算開始時(shí)需指定每個(gè)粒子點(diǎn)的初始?jí)毫?

    文中計(jì)算權(quán)函數(shù)w(x)選取新4次權(quán)重函數(shù),因?yàn)樵摵瘮?shù)本身是1條連續(xù)曲線,且2階導(dǎo)數(shù)曲線具有較好的光滑性.具體取為

    (17)

    式中,d=x-xi,h為支持半徑.

    1.4 自由面判斷

    FPM 對(duì)自由液面的判斷采用粒子數(shù)密度的概念來區(qū)分自由液面與內(nèi)部流體粒子,當(dāng)粒子數(shù)密度〈ni〉滿足:

    〈ni〉<βn0

    (18)

    即可判定為自由面粒子,其中n0為初始粒子數(shù)密度.在求解壓力 Poisson 方程時(shí),自由面粒子的壓力賦值為 0.β是一參數(shù),對(duì)自由面的判斷有一定的影響, 一般取值為β=0.80~0.99.

    2數(shù)值算例

    為驗(yàn)證文中提出的有限點(diǎn)法的有效性,依據(jù)文獻(xiàn)[21]中的試驗(yàn),選取4個(gè)不同工況的液艙晃蕩模型(見表1),對(duì)其進(jìn)行數(shù)值模擬.同時(shí)在計(jì)算中設(shè)置壓力監(jiān)測(cè)點(diǎn)O,以便與試驗(yàn)中測(cè)得的艙壁壓力進(jìn)行對(duì)比,定量分析數(shù)值計(jì)算結(jié)果的準(zhǔn)確性.

    表1 液艙晃蕩計(jì)算工況Table 1 Calculation conditions of tank sloshing

    液艙晃蕩計(jì)算模型與試驗(yàn)[21]保持一致,取二維矩形艙,寬度為0.6 m,高度為0.3 m.液體的晃蕩由液艙做規(guī)則橫蕩運(yùn)動(dòng)激發(fā)產(chǎn)生,其表達(dá)式為

    (19)

    式中:A為運(yùn)動(dòng)幅值;Amax為振幅,取值0.05 m.具體計(jì)算模型見圖1.

    圖1 液艙模型及模擬工況示意圖(單位:m)Fig.1 Schematic sketch of the liquid tank models and simulation conditions(Unit:m)

    計(jì)算中物理參數(shù)的設(shè)置為:水的密度為 1 000 kg/m3,運(yùn)動(dòng)粘性系數(shù)取 1.0×10-6m2/s,重力加速度為 9.81 m/s2,粒子支持域半徑為 0.007 m,時(shí)間步長(zhǎng)取為0.001 s,計(jì)算中分別數(shù)值模擬了10個(gè)周期的晃蕩運(yùn)動(dòng).

    2.1 晃蕩模型A

    在晃蕩模型A中,液艙充容率40%為中等充容,激勵(lì)頻率小于共振頻率.圖2為晃蕩模型A在一個(gè)周期內(nèi)的不同時(shí)刻(t=0.1T,0.2T,0.3T和0.4T)運(yùn)動(dòng)流場(chǎng)的計(jì)算結(jié)果,以壓力分布的形式給出.由圖2可看到,在t=0.1T時(shí), 液艙內(nèi)流體在激勵(lì)力的作用下向右側(cè)壁運(yùn)動(dòng),并伴有行進(jìn)波.當(dāng)達(dá)到 0.2T時(shí),撞壁流體在慣性的作用下開始沿壁面上升.當(dāng)t=0.3T時(shí), 沿右側(cè)壁面運(yùn)動(dòng)的流體再次撞擊到液艙頂棚,同時(shí)底部的液體由于越積越多,并在液艙的激勵(lì)運(yùn)動(dòng)下愈要形成反向運(yùn)動(dòng)的“水包”.當(dāng)接近 0.4T時(shí),抨擊右側(cè)壁面的流體以凸起“水包”狀開始向艙壁左側(cè)運(yùn)動(dòng).在艙壁激勵(lì)力的作用下,流體發(fā)生了不同形態(tài)的晃蕩,出現(xiàn)了行進(jìn)波、水躍、翻卷等強(qiáng)非線性現(xiàn)象.這些復(fù)雜變形的水面通過FPM得以準(zhǔn)確表達(dá),相比于傳統(tǒng)基于空間網(wǎng)格的數(shù)值模擬方法,有限點(diǎn)法不但消除了對(duì)流項(xiàng)的數(shù)值耗散問題,還具有較強(qiáng)的自由面捕捉能力,這對(duì)自由面大幅運(yùn)動(dòng)的類似求解方法具有一定的指導(dǎo)意義.

    圖2 模型 A 在不同時(shí)刻的運(yùn)動(dòng)流場(chǎng)示意Fig.2 Numerical results of liquid sloshing case A at different time

    為定量分析FPM在求解液艙晃蕩問題時(shí)的準(zhǔn)確性,圖3給出了壓力監(jiān)測(cè)點(diǎn)O在不同時(shí)刻的試驗(yàn)測(cè)量值[21]與數(shù)值結(jié)果.由圖3可知,數(shù)值計(jì)算的壓力峰值產(chǎn)生了一定幅度的非物理震蕩,前4個(gè)周期內(nèi)壓力峰值較試驗(yàn)值偏大.對(duì)第1個(gè)周期內(nèi)的兩條曲線進(jìn)行對(duì)比分析可知,試驗(yàn)測(cè)量峰值(1 300 N/m2)約為數(shù)值計(jì)算峰值(2 400 N/m2)的1/2.在之后的周期內(nèi),數(shù)值計(jì)算結(jié)果與試驗(yàn)測(cè)量值基本保持一致.對(duì)比兩條曲線可知,每個(gè)周期內(nèi)均有兩個(gè)壓力峰值,且第2峰值的數(shù)值結(jié)果與試驗(yàn)值基本相同,壓力的整體變化趨勢(shì)基本吻合,從定量角度驗(yàn)證了FPM對(duì)解決液艙晃蕩問題具有很好的可靠性.

    圖3 模型A監(jiān)測(cè)點(diǎn)O壓力計(jì)算結(jié)果與試驗(yàn)值對(duì)比

    Fig.3 Pressure comparison between calculation and experiment for measuring pointOof case A

    2.2 晃蕩模型B

    晃蕩模型B為中等充容,同樣液艙充容率為40%,激勵(lì)周期為1.3 s.圖4列出了試驗(yàn)[21]與數(shù)值場(chǎng)結(jié)果在4個(gè)時(shí)刻(t=0.1T,0.2T,0.3T和0.4T)的流場(chǎng)對(duì)比,其中數(shù)值計(jì)算的瞬時(shí)流場(chǎng)同樣以壓力形式給出.從圖4可看到,由于激勵(lì)頻率與共振頻率相同,液艙的橫蕩運(yùn)動(dòng)激發(fā)液體發(fā)生劇烈的晃蕩,自由液面凹凸明顯;尤其當(dāng)t=0.2T時(shí), 靠近液艙右側(cè)壁的流體壓強(qiáng)激增,形成較大的壓力梯度.當(dāng)達(dá)到0.4T時(shí),由于頂部與左側(cè)艙壁與液體發(fā)生強(qiáng)烈的拍擊,促使流體翻卷下落,并伴有迸濺的水花.通過將數(shù)值解與試驗(yàn)圖像定性對(duì)比可知,這些復(fù)雜變形通過有限點(diǎn)法得到了準(zhǔn)確的捕捉,且與試驗(yàn)結(jié)果吻合較好,驗(yàn)證了FPM在解決液艙晃蕩問題中的可靠性.

    圖4 模型 B 在不同時(shí)刻(t=0.1T,0.2T,0.3T和0.4T)的運(yùn)動(dòng)流場(chǎng)試驗(yàn)對(duì)比Fig.4 Qualitative comparison between FPM and experiment att=0.1T, 0.2T, 0.3Tand 0.4T, Case B-sloshing in medium-filling tank

    2.3 晃蕩模型C

    在前面的晃蕩模型A、B中,討論的是液艙中等充容情況,下面給出高充容率的液艙晃蕩計(jì)算結(jié)果.晃蕩模型C中的充容率高達(dá)83%,圖5為試驗(yàn)[21]與數(shù)值結(jié)果在4個(gè)時(shí)刻(t=0.02T,0.17T,0.35T和0.50T)的流場(chǎng)對(duì)比,其中數(shù)值結(jié)果同樣以壓力場(chǎng)形式給出.從圖5可看到,在t=0.02T時(shí),液艙的激勵(lì)作用使得流體開始做晃蕩運(yùn)動(dòng),水面呈凹陷狀.在接近t=0.17T時(shí),液艙內(nèi)運(yùn)動(dòng)流體逐漸靠近右側(cè)艙壁,并伴有飛濺的液滴.直到t=0.35T時(shí),由于右側(cè)壁面的阻擋, 在慣性的驅(qū)使下,使得撞擊到右側(cè)艙壁的液體沿壁面上升,并再次與上壁面發(fā)生拍擊.當(dāng)達(dá)到t=0.50T時(shí),沿上壁面運(yùn)動(dòng)的流體由于受到重性力的作用,開始下落, 并與從底部向上運(yùn)動(dòng)的水流匯集,融合的水面發(fā)生極復(fù)雜的變形.觀察整個(gè)運(yùn)動(dòng)過程,由于液艙的充容率很高,其橫蕩運(yùn)動(dòng)使得流體與左右艙壁發(fā)生了強(qiáng)烈的拍擊,液面變化更加明顯,其劇烈程度明顯高于中等充容模型.通過與試驗(yàn)圖片的對(duì)比可知,數(shù)值模擬結(jié)果與試驗(yàn)值基本吻合,說明有限點(diǎn)法對(duì)解決高充容率的液艙晃蕩具有很好的數(shù)值模擬能力.

    圖5 模型 C 在不同時(shí)刻(t=0.02T,0.17T,0.35T和0.50T)的運(yùn)動(dòng)流場(chǎng)試驗(yàn)對(duì)比Fig.5 Qualitative comparison between FPM and experiment att=0.02T, 0.17T, 0.35Tand 0.50T, Case C-slo-shing in high-filling tank

    2.4 晃蕩模型D

    晃蕩模型D與晃蕩模型C的充容率相同,但激勵(lì)周期為1.0 s.本晃蕩模型的計(jì)算是為了從定量角度出發(fā),來驗(yàn)證有限點(diǎn)法在數(shù)值計(jì)算高充容液艙晃蕩問題時(shí)的準(zhǔn)確性.為與試驗(yàn)[21]保持一致,計(jì)算中的壓力監(jiān)測(cè)點(diǎn)O高度為0.235 m.下面給出晃蕩模型D中壓力監(jiān)測(cè)點(diǎn)在不同時(shí)刻的試驗(yàn)值[21]與數(shù)值結(jié)果,詳見圖6.

    圖6 模型D監(jiān)測(cè)點(diǎn)O壓力計(jì)算結(jié)果與試驗(yàn)值對(duì)比Fig.6 Measuring pointOof case D pressure comparison in between calculation and experiment

    由圖6可知,壓力數(shù)值計(jì)算峰值較實(shí)測(cè)壓力值略有偏大.分析其原因主要是由于數(shù)值模擬時(shí)采用的是單相流,而在真實(shí)試驗(yàn)中當(dāng)晃蕩流體撞擊到壁面時(shí),由于下落、翻卷、融合的過程中融入一定組分的空氣,促使流體的拍擊壓力起到了緩沖作用,造成計(jì)算結(jié)果大于試驗(yàn)值.從壓力的相位角度觀察,圖中前幾個(gè)周期內(nèi)有一定的相位差存在,隨著時(shí)間的推移,相位差基本消除,峰值也漸趨保持一致.總體來說,數(shù)值結(jié)果與試驗(yàn)基本吻合,說明FPM對(duì)于液艙晃蕩問題的求解具有較好的準(zhǔn)確性.

    3結(jié)語(yǔ)

    文中采用有限點(diǎn)法對(duì)不同充容率的二維液艙在不同頻率下的受迫運(yùn)動(dòng)進(jìn)行了數(shù)值模擬研究.結(jié)果表明:當(dāng)激勵(lì)頻率與液艙固有頻率一致時(shí),艙內(nèi)流體的晃蕩較為劇烈;而當(dāng)充容率不同時(shí),由于艙壁對(duì)液體流動(dòng)的抑制作用,使得高充容液艙內(nèi)的流體翻卷變形較為平緩.通過有限點(diǎn)法準(zhǔn)確有效地捕捉到液艙在不同激勵(lì)頻率下的復(fù)雜水面晃蕩情況,如水躍、翻卷以及融合等強(qiáng)非線性現(xiàn)象,F(xiàn)PM均給予較為準(zhǔn)確的模擬,與試驗(yàn)現(xiàn)象基本吻合.與此同時(shí),在晃蕩模型B和D中,通過分析監(jiān)測(cè)點(diǎn)的壓力曲線,進(jìn)一步驗(yàn)證有限點(diǎn)法在定量分析上也有滿意的精度.由于模型D為高充容狀態(tài),相對(duì)模型B而言,空氣組分影響較小,所以模型D中壓力曲線更為準(zhǔn)確.為提高數(shù)值模擬精度,今后可在FPM法中采用多相流模型.縱觀定性對(duì)比與定量分析,有限點(diǎn)法在解決不同工況的二維液艙晃蕩問題時(shí)均得到較為準(zhǔn)確的模擬結(jié)果,說明其在求解液艙晃蕩問題中是可靠、有效的.

    參考文獻(xiàn):

    [1]Pitblado Robin M,Woodward John L.Highlights of LNG risk technology [J].Journal of Loss Prevention in the Process Industries,2011,24(6):827-836.

    [2]Pistani Fabrizio,Thiagarajan Krish.Experimental mea-surements and data analysis of the impact pressures in a sloshing experiment [J].Ocean Engineering,2012,52(1):60-74.

    [3]Armenio V.An improved MAC method (SIMAC) for unsteady high-reynolds free surface flows [J].International Journal for Numerical Methods in Fluids,1998,24(2):185-214.

    [4]Kim Y.Numerical simulation of sloshing flows with impact load [J].Applied Ocean Research,2001,23(1):53-62.

    [5]Kim Y,Shin Y S,Lee K H.Numerical study on slosh-induced impact pressure on three-dimensional prismatic tanks [J].Applied Ocean Research,2004,26(5):213-226.

    [6]Kim Y,Nam B W,Kim D W,et al.Study on coupling effects of ship motion and sloshing [J].Ocean Engineering,2007,34(16):2176-2187.

    [7]Hu C H,Yang K K,Kim Y H.3-D numerical simulations of violent sloshing by CIP-based method [J].Journal of Hydrodynamics:Ser B,2010,22(5):253-258.

    [8]祁江濤,顧民,吳乘勝.液艙晃蕩的數(shù)值模擬 [J].船舶力學(xué),2008,12(4):574-581.

    Qi Jiang-tao,Gu Min,Wu Cheng-sheng.Numerical simulation of sloshing in liquid tank [J].Journal of Ship Mechanics,2008,12(4):574-581.

    [9]沈猛.基于改進(jìn)VOF法的棱形液艙晃蕩分析及應(yīng)用 [D].上海:上海交通大學(xué)船舶海洋與建筑工程學(xué)院,2008.

    [10]沈猛,王剛,唐文勇.基于改進(jìn)VOF法夫人棱形液艙液體晃蕩分析 [J].中國(guó)造船,2009,50(1):1-9.

    Shen Meng,Wang Gang,Tang Wen-yong.Liguid sloshing analysis on prismatic tanks based on improved VOF method [J].Shipbuilding of China,2009,50(1):1-9.

    [11]方智勇,楊松林,朱仁慶.Level-set法在液體晃蕩研究中的應(yīng)用 [J].水動(dòng)力研究與進(jìn)展,2007,22(2):150-156.

    Fang Zhi-yong,Yang Song-lin,Zhu Ren-qing.Application of Level-set method in the research on liquid sloshing [J].Journal of Hydrodynamics,2007,22(2):150-156.

    [12]方智勇,端木玉,朱仁慶. 基于Level-set 法的液艙液體晃蕩數(shù)值模擬 [J].船舶力學(xué),2007,11(1):62-67.

    Fang Zhi-yong,Duan Mu-yu,Zhu Ren-qing.Numerical simulation of liquid sloshing in a liquid tank based on Level-set method [J].Journal of Ship Mechanics,2007,11(1):62-67.

    [13]趙相坤,李鳳霞,戰(zhàn)守義.基于GPU的面向SPH流體模擬的鄰居查找算法 [J].華南理工大學(xué)學(xué)報(bào):自然科學(xué)版,2011,39(7):150-155.

    Zhao Xiang-kun,Li Feng-xia,Zhan Shou-yi.CPU-based neighbor search algorithm for SPH fluid simulations [J].Journal of South China University of Technology:Natural Science Edition,2011,39(7):150-155.

    [17]Shu C,Ding H,Chen H Q,et al.An upwind local RBF-DQ method for simulation of inviscid compressible flows [J].Computer Methods in Applied Mechanics and Engineering,2005,194(18-20):2001-2017.

    [18]Pandey A,Klar A,Tiwari S.Meshfree method for fluctuating hydrodynamics [J].Mathematics and Computers in Simulation,2012,82(11):2157-2166.

    [19]Reséndiz-Flores E O,García-Calvillo I D.Application of the finite pointset method to non-stationary heat conduction problems [J].International Journal of Heat and Mass Transfer,2014,71:720-723.

    [20]Gu T X,Zuo X Y,Liu X P,et al.An improved parallel hybrid bi-conjugate gradient method suitable for distributed parallel computing [J].Journal of Computational and Applied Mathematics,2009,226(1):55-65.

    [21]Kishev Z R,Hu C,Kashiwagi M.Numerical simulation of violent sloshing by a CIP-based method [J].Journal of marine science and technology,2006,11(2):111-122.

    Application of Finite Point Method to Liquid Sloshing in Tanks

    LuYu1HuAn-kang1,2LiuYa-chong1

    (1.College of Shipbuilding Engineering, Harbin Engineering University, Harbin 150001, Heilongjiang, China;

    2.CIMC Ocean Engineering Design & Research Institute Co., Ltd., Shanghai 201206, China)

    Abstract:In recent years, the structural design of liquid cargo ships has been focused on the liquid sloshing in tanks, which involves the high non-linearity problem related to the large amplitude motion of free liquid surface and therefore causes some difficulties in Euler numerical simulation depending upon spatial topological meshes. In this paper, the liquid sloshing of a tank at different filling ratios is numerically simulated by means of the Lagrangian view-based finite point method (FPM), and the numerical results are compared with the experimental ones. The results show that when the excitation frequency is consistent with the natural frequency of the tank, the fluid sloshing is severer; and that, at different filling ratios, the fluid in the tank with a high filling ratio will flow with mode-rate deformations due to the inhibition of bulkheads. Moreover, the FPM can effectively and accurately simulate the violent changes of the liquid sloshing of the tank at different filling ratios under the multi-frequency excitation, and the simulation results accord well with the experimental phenomena. Meanwhile, an impact pressure on the bulkheads is observed and the calculated values obtained through the FPM are consistent with the measured pressure curves. It is, therefore, found that the FPM performs accurate simulations on the two-dimensional liquid sloshing in the tank under different conditions, which effectively verifies the reliability of the FPM.

    Key words:liquid sloshing in tanks; Lagrangian view; finite point method; moving least square; Poisson equation

    中圖分類號(hào):U661.1

    doi:10.3969/j.issn.1000-565X.2015.08.021

    文章編號(hào):1000-565X(2015)08-0144-07

    作者簡(jiǎn)介:盧雨(1988-),男,博士生,主要從事船舶流體力學(xué)研究.E-mail: luyu90627@126.com?通信作者: 胡安康(1956-),女,博士,教授,主要從事船舶總體設(shè)計(jì)理論與方法研究.E-mail: ankang.hu@cimc.com

    *基金項(xiàng)目:國(guó)家自然科學(xué)基金資助項(xiàng)目(51379040)

    收稿日期:2015-01-12

    Foundation item: Supported by the National Natural Science Foundation of China(51379040)

    97在线视频观看| 99久久人妻综合| 亚洲国产毛片av蜜桃av| 岛国毛片在线播放| 精品人妻偷拍中文字幕| 中国三级夫妇交换| 9191精品国产免费久久| 男人舔女人的私密视频| 人人妻人人爽人人添夜夜欢视频| 日韩av在线免费看完整版不卡| 精品午夜福利在线看| 日韩成人av中文字幕在线观看| 亚洲内射少妇av| 99久久中文字幕三级久久日本| 亚洲欧洲日产国产| 久久97久久精品| 免费日韩欧美在线观看| 一区二区三区乱码不卡18| 一区二区三区乱码不卡18| 亚洲成人一二三区av| 97精品久久久久久久久久精品| 18禁动态无遮挡网站| 免费黄色在线免费观看| 国产精品久久久久成人av| 亚洲美女视频黄频| 中国三级夫妇交换| 国产男人的电影天堂91| 精品久久国产蜜桃| 国产高清国产精品国产三级| 久久午夜综合久久蜜桃| 久久久久人妻精品一区果冻| av卡一久久| av在线播放精品| 亚洲 欧美一区二区三区| av在线播放精品| 亚洲,一卡二卡三卡| 中文字幕精品免费在线观看视频 | 高清在线视频一区二区三区| 欧美最新免费一区二区三区| 国产一区二区在线观看av| 在线观看国产h片| 最近的中文字幕免费完整| 成年女人在线观看亚洲视频| 中文天堂在线官网| 日韩一本色道免费dvd| 观看av在线不卡| 男女下面插进去视频免费观看 | av国产久精品久网站免费入址| 性色avwww在线观看| 午夜福利影视在线免费观看| 国产精品99久久99久久久不卡 | 精品国产一区二区久久| 七月丁香在线播放| 黑人猛操日本美女一级片| 好男人视频免费观看在线| 99热6这里只有精品| 欧美日韩亚洲高清精品| 成年女人在线观看亚洲视频| 午夜免费男女啪啪视频观看| 捣出白浆h1v1| 高清黄色对白视频在线免费看| 亚洲熟女精品中文字幕| 大码成人一级视频| 久久人人爽av亚洲精品天堂| av卡一久久| 人妻人人澡人人爽人人| 一本大道久久a久久精品| 色婷婷久久久亚洲欧美| 久久久久国产精品人妻一区二区| 国产精品一二三区在线看| 伊人久久国产一区二区| 韩国精品一区二区三区 | 国产永久视频网站| 一级毛片黄色毛片免费观看视频| 欧美xxⅹ黑人| 这个男人来自地球电影免费观看 | 国产色爽女视频免费观看| 亚洲欧美色中文字幕在线| 国产成人精品一,二区| 日韩欧美一区视频在线观看| 伦精品一区二区三区| 亚洲中文av在线| 亚洲国产成人一精品久久久| 国产色爽女视频免费观看| 欧美+日韩+精品| 免费人妻精品一区二区三区视频| 日本爱情动作片www.在线观看| 999精品在线视频| 午夜福利,免费看| 在线观看免费日韩欧美大片| 久久精品久久精品一区二区三区| 久久精品国产a三级三级三级| 大陆偷拍与自拍| 色吧在线观看| 大香蕉97超碰在线| 考比视频在线观看| 亚洲精品日韩在线中文字幕| 2018国产大陆天天弄谢| a级毛片黄视频| 亚洲国产看品久久| 蜜臀久久99精品久久宅男| 精品第一国产精品| 另类精品久久| 久久久久久久久久久免费av| 国产欧美另类精品又又久久亚洲欧美| a 毛片基地| 精品人妻在线不人妻| 这个男人来自地球电影免费观看 | 国产欧美另类精品又又久久亚洲欧美| 一二三四在线观看免费中文在 | 精品福利永久在线观看| 精品国产一区二区三区久久久樱花| 国产精品熟女久久久久浪| 精品国产国语对白av| 亚洲情色 制服丝袜| 91精品伊人久久大香线蕉| 久久国内精品自在自线图片| 人妻系列 视频| 国产色婷婷99| 欧美激情国产日韩精品一区| 国产免费一级a男人的天堂| 2022亚洲国产成人精品| 极品人妻少妇av视频| 免费久久久久久久精品成人欧美视频 | av播播在线观看一区| 国产精品国产三级国产av玫瑰| 中国国产av一级| 性色av一级| 亚洲国产日韩一区二区| 十分钟在线观看高清视频www| 欧美国产精品va在线观看不卡| 精品视频人人做人人爽| 如日韩欧美国产精品一区二区三区| 久久 成人 亚洲| 亚洲人与动物交配视频| 日韩一区二区三区影片| 搡女人真爽免费视频火全软件| 极品人妻少妇av视频| 日韩大片免费观看网站| 搡老乐熟女国产| 91在线精品国自产拍蜜月| 中文字幕精品免费在线观看视频 | 国产av精品麻豆| 满18在线观看网站| 黄片无遮挡物在线观看| 国产精品人妻久久久影院| 久久热在线av| 亚洲精品日本国产第一区| 成人二区视频| 久久婷婷青草| 人人妻人人爽人人添夜夜欢视频| 国产精品 国内视频| 极品人妻少妇av视频| 亚洲欧洲日产国产| 久久国产亚洲av麻豆专区| 国产精品嫩草影院av在线观看| 中文字幕免费在线视频6| 深夜精品福利| 免费女性裸体啪啪无遮挡网站| 国产黄色免费在线视频| 一级毛片电影观看| 日韩免费高清中文字幕av| 免费观看性生交大片5| av女优亚洲男人天堂| 久久久久久久大尺度免费视频| 赤兔流量卡办理| 日韩熟女老妇一区二区性免费视频| 2022亚洲国产成人精品| 在线观看三级黄色| 一区二区三区四区激情视频| a级毛片黄视频| 国产69精品久久久久777片| 精品99又大又爽又粗少妇毛片| 久久久精品免费免费高清| 日韩一本色道免费dvd| 国国产精品蜜臀av免费| 最后的刺客免费高清国语| 久热久热在线精品观看| 久久av网站| 在线看a的网站| 国产精品成人在线| 精品国产一区二区久久| 久久热在线av| 色婷婷久久久亚洲欧美| 如何舔出高潮| 日韩一区二区视频免费看| 亚洲成色77777| 亚洲精品乱久久久久久| 亚洲 欧美一区二区三区| 一边摸一边做爽爽视频免费| av在线app专区| 国产成人精品福利久久| 最新的欧美精品一区二区| 国产亚洲欧美精品永久| 亚洲欧美日韩另类电影网站| 夫妻性生交免费视频一级片| 精品第一国产精品| 国产免费又黄又爽又色| 亚洲色图 男人天堂 中文字幕 | 青春草视频在线免费观看| 免费人妻精品一区二区三区视频| 国语对白做爰xxxⅹ性视频网站| 成年女人在线观看亚洲视频| 国产麻豆69| 熟女人妻精品中文字幕| 久久久久久久久久人人人人人人| 久久人人爽av亚洲精品天堂| 青春草视频在线免费观看| av免费观看日本| 在现免费观看毛片| 黄色配什么色好看| 97在线视频观看| 韩国精品一区二区三区 | 久热久热在线精品观看| 国产永久视频网站| 在线观看免费高清a一片| 午夜福利乱码中文字幕| 国产麻豆69| 亚洲欧美日韩另类电影网站| 人体艺术视频欧美日本| 人人妻人人添人人爽欧美一区卜| 亚洲国产看品久久| 亚洲国产av影院在线观看| 国产成人aa在线观看| 亚洲成色77777| 男女高潮啪啪啪动态图| 18禁观看日本| freevideosex欧美| 我要看黄色一级片免费的| 免费在线观看完整版高清| 亚洲一级一片aⅴ在线观看| 啦啦啦啦在线视频资源| 欧美日韩综合久久久久久| 国产欧美日韩综合在线一区二区| 国产精品嫩草影院av在线观看| 捣出白浆h1v1| 搡老乐熟女国产| 在线观看人妻少妇| 69精品国产乱码久久久| a 毛片基地| 免费av不卡在线播放| 精品一区在线观看国产| 国产精品麻豆人妻色哟哟久久| 91国产中文字幕| 日本欧美国产在线视频| 亚洲人成77777在线视频| 成人漫画全彩无遮挡| 精品亚洲成国产av| 国产欧美日韩综合在线一区二区| 久久精品久久久久久久性| 国产成人精品久久久久久| 啦啦啦中文免费视频观看日本| 丁香六月天网| 亚洲av成人精品一二三区| 99热全是精品| 男的添女的下面高潮视频| 91精品伊人久久大香线蕉| kizo精华| 国产激情久久老熟女| 亚洲伊人色综图| 丝袜人妻中文字幕| av女优亚洲男人天堂| 国产免费现黄频在线看| 看免费av毛片| 亚洲性久久影院| 国产成人精品婷婷| 国产亚洲精品第一综合不卡 | 日日啪夜夜爽| 国产一区二区三区av在线| 美女国产高潮福利片在线看| 久久久久国产网址| 亚洲熟女精品中文字幕| 久久综合国产亚洲精品| 色婷婷久久久亚洲欧美| 欧美日韩亚洲高清精品| 久久久欧美国产精品| 男女啪啪激烈高潮av片| 黑人猛操日本美女一级片| 下体分泌物呈黄色| 国产av一区二区精品久久| 国产有黄有色有爽视频| 色婷婷av一区二区三区视频| 亚洲婷婷狠狠爱综合网| 亚洲av综合色区一区| av有码第一页| 免费播放大片免费观看视频在线观看| 国精品久久久久久国模美| 欧美 日韩 精品 国产| 成人免费观看视频高清| 综合色丁香网| 国产午夜精品一二区理论片| 一级片'在线观看视频| 日日啪夜夜爽| www.色视频.com| 欧美成人精品欧美一级黄| 一级,二级,三级黄色视频| 我要看黄色一级片免费的| 国产白丝娇喘喷水9色精品| 免费大片18禁| 久久av网站| 精品久久国产蜜桃| 国产淫语在线视频| 午夜激情久久久久久久| 十八禁高潮呻吟视频| 欧美人与性动交α欧美精品济南到 | 黑人猛操日本美女一级片| 少妇精品久久久久久久| 成人国产麻豆网| 超碰97精品在线观看| 男女啪啪激烈高潮av片| 水蜜桃什么品种好| 欧美精品人与动牲交sv欧美| 美国免费a级毛片| 岛国毛片在线播放| 日本爱情动作片www.在线观看| 热99久久久久精品小说推荐| 精品一区在线观看国产| 午夜免费男女啪啪视频观看| av一本久久久久| 夜夜爽夜夜爽视频| 亚洲国产精品专区欧美| 久久久久久人妻| 日本色播在线视频| 国产无遮挡羞羞视频在线观看| 1024视频免费在线观看| 在线免费观看不下载黄p国产| 制服丝袜香蕉在线| 老司机亚洲免费影院| 日韩中文字幕视频在线看片| videos熟女内射| 26uuu在线亚洲综合色| 韩国av在线不卡| 久久久久久久久久久久大奶| 九色亚洲精品在线播放| 国产精品 国内视频| 国产片特级美女逼逼视频| 免费女性裸体啪啪无遮挡网站| 9191精品国产免费久久| 女性生殖器流出的白浆| 久久国产亚洲av麻豆专区| 日本黄色日本黄色录像| 日本爱情动作片www.在线观看| 国产伦理片在线播放av一区| 黑人欧美特级aaaaaa片| 亚洲,欧美精品.| 狂野欧美激情性xxxx在线观看| 亚洲精品久久成人aⅴ小说| 大陆偷拍与自拍| 亚洲色图综合在线观看| 黄色配什么色好看| 91精品伊人久久大香线蕉| 国产一区二区三区综合在线观看 | 国产免费现黄频在线看| 毛片一级片免费看久久久久| 少妇高潮的动态图| 亚洲婷婷狠狠爱综合网| 国产av精品麻豆| 欧美日韩国产mv在线观看视频| 国产精品.久久久| 亚洲精品久久午夜乱码| 久久久久精品人妻al黑| 一二三四中文在线观看免费高清| 精品酒店卫生间| 免费观看性生交大片5| 精品人妻熟女毛片av久久网站| 黄片播放在线免费| 亚洲成av片中文字幕在线观看 | 国产成人a∨麻豆精品| 看免费成人av毛片| 精品少妇内射三级| 插逼视频在线观看| 久久精品国产综合久久久 | 老司机影院成人| 亚洲国产精品一区二区三区在线| 精品一区二区三区视频在线| 欧美精品一区二区大全| 一二三四在线观看免费中文在 | 国产在视频线精品| 国产成人精品久久久久久| 亚洲av免费高清在线观看| 你懂的网址亚洲精品在线观看| 久久精品国产综合久久久 | 亚洲欧美一区二区三区国产| 久久久亚洲精品成人影院| 亚洲av免费高清在线观看| 大话2 男鬼变身卡| 水蜜桃什么品种好| 免费人妻精品一区二区三区视频| 国产精品久久久久久久久免| 免费在线观看黄色视频的| 中文字幕免费在线视频6| 日韩精品有码人妻一区| 国产在线一区二区三区精| av片东京热男人的天堂| 亚洲少妇的诱惑av| 各种免费的搞黄视频| 天天影视国产精品| 黑人猛操日本美女一级片| 热99久久久久精品小说推荐| 在线观看美女被高潮喷水网站| 91精品三级在线观看| av在线播放精品| 五月开心婷婷网| 91精品国产国语对白视频| 黄色配什么色好看| 久久国产亚洲av麻豆专区| 久久女婷五月综合色啪小说| 一本久久精品| 中文天堂在线官网| 国产精品一二三区在线看| 午夜av观看不卡| 国产精品蜜桃在线观看| www.色视频.com| 午夜福利视频在线观看免费| 有码 亚洲区| 国产精品国产av在线观看| 伦理电影免费视频| 精品人妻熟女毛片av久久网站| 中文字幕av电影在线播放| 人人妻人人澡人人爽人人夜夜| 男男h啪啪无遮挡| 久久女婷五月综合色啪小说| 亚洲av免费高清在线观看| 成人毛片a级毛片在线播放| 亚洲欧美中文字幕日韩二区| 丝袜人妻中文字幕| 女人被躁到高潮嗷嗷叫费观| 少妇人妻久久综合中文| 亚洲国产看品久久| 亚洲av电影在线观看一区二区三区| 亚洲成av片中文字幕在线观看 | 午夜免费鲁丝| 成人免费观看视频高清| 一本久久精品| 免费观看av网站的网址| 久久久久国产精品人妻一区二区| 制服诱惑二区| 伊人久久国产一区二区| 一边亲一边摸免费视频| 少妇的逼水好多| 在线观看国产h片| 亚洲精品美女久久av网站| 大片免费播放器 马上看| 九九在线视频观看精品| 成年人午夜在线观看视频| 男女边摸边吃奶| 国产无遮挡羞羞视频在线观看| 国产极品粉嫩免费观看在线| 国产精品久久久av美女十八| 夫妻午夜视频| 高清视频免费观看一区二区| 精品卡一卡二卡四卡免费| 午夜福利在线观看免费完整高清在| 亚洲精品久久午夜乱码| 高清欧美精品videossex| 国产一区二区激情短视频 | 18禁国产床啪视频网站| 国产欧美另类精品又又久久亚洲欧美| 午夜精品国产一区二区电影| 22中文网久久字幕| 免费观看性生交大片5| 午夜福利乱码中文字幕| 国产一区二区在线观看av| 亚洲av综合色区一区| 边亲边吃奶的免费视频| 一本久久精品| 亚洲中文av在线| 黑人欧美特级aaaaaa片| 80岁老熟妇乱子伦牲交| 亚洲成国产人片在线观看| 高清在线视频一区二区三区| 亚洲色图 男人天堂 中文字幕 | 精品人妻熟女毛片av久久网站| 久久99热这里只频精品6学生| 国产精品一区二区在线观看99| 一本大道久久a久久精品| 久久久亚洲精品成人影院| 国产一级毛片在线| 久久久久久伊人网av| 亚洲国产看品久久| 人人妻人人添人人爽欧美一区卜| 99国产综合亚洲精品| 男人爽女人下面视频在线观看| 免费高清在线观看视频在线观看| 高清毛片免费看| 啦啦啦在线观看免费高清www| 国产又色又爽无遮挡免| 亚洲精品av麻豆狂野| 免费女性裸体啪啪无遮挡网站| 亚洲精品第二区| 日韩中字成人| 久久国产精品男人的天堂亚洲 | 免费av不卡在线播放| 国产片特级美女逼逼视频| 久久99蜜桃精品久久| 久久精品国产鲁丝片午夜精品| 免费久久久久久久精品成人欧美视频 | av网站免费在线观看视频| 丰满饥渴人妻一区二区三| 亚洲欧洲国产日韩| 在线观看免费视频网站a站| 校园人妻丝袜中文字幕| 精品熟女少妇av免费看| 欧美丝袜亚洲另类| 美女国产视频在线观看| 考比视频在线观看| 毛片一级片免费看久久久久| 国产永久视频网站| 18禁观看日本| 午夜久久久在线观看| 捣出白浆h1v1| 老司机亚洲免费影院| 91aial.com中文字幕在线观看| 国产在线一区二区三区精| 国产精品三级大全| 国产精品.久久久| av女优亚洲男人天堂| 国产熟女午夜一区二区三区| 国产精品久久久av美女十八| 青青草视频在线视频观看| 香蕉丝袜av| 亚洲精品第二区| 久久这里只有精品19| 精品人妻在线不人妻| 黄色毛片三级朝国网站| 国产福利在线免费观看视频| 成人二区视频| 在线观看三级黄色| 欧美日本中文国产一区发布| 老司机影院成人| 成人黄色视频免费在线看| 亚洲国产毛片av蜜桃av| 人人妻人人添人人爽欧美一区卜| 欧美人与善性xxx| 一本大道久久a久久精品| 大香蕉久久成人网| 一个人免费看片子| 在线观看美女被高潮喷水网站| 一级毛片黄色毛片免费观看视频| 久久精品久久久久久久性| 十分钟在线观看高清视频www| 亚洲精品一二三| 丰满少妇做爰视频| 亚洲精品,欧美精品| 乱人伦中国视频| 激情五月婷婷亚洲| 欧美最新免费一区二区三区| 少妇 在线观看| 精品一区二区免费观看| a级片在线免费高清观看视频| 亚洲婷婷狠狠爱综合网| 丝袜人妻中文字幕| 精品亚洲成国产av| 欧美日韩视频高清一区二区三区二| 你懂的网址亚洲精品在线观看| 在线观看www视频免费| 99久久中文字幕三级久久日本| 午夜福利,免费看| 国产亚洲精品久久久com| 欧美精品一区二区免费开放| 免费少妇av软件| 免费观看性生交大片5| 天天操日日干夜夜撸| 国产欧美日韩一区二区三区在线| 午夜激情久久久久久久| 国产精品.久久久| 少妇 在线观看| 桃花免费在线播放| 国产毛片在线视频| 老司机影院成人| 免费在线观看黄色视频的| 国产av一区二区精品久久| 国产免费一区二区三区四区乱码| 免费高清在线观看视频在线观看| 国产色爽女视频免费观看| 精品少妇久久久久久888优播| 国产成人精品久久久久久| 亚洲国产成人一精品久久久| 免费av中文字幕在线| 久久女婷五月综合色啪小说| 国产日韩欧美视频二区| 日韩三级伦理在线观看| 久久久久国产网址| 一二三四在线观看免费中文在 | 如何舔出高潮| 亚洲欧洲日产国产| 黑丝袜美女国产一区| 欧美丝袜亚洲另类| 久久人人爽人人爽人人片va| videosex国产| 国产精品国产三级国产专区5o| 黑人高潮一二区| 黄色配什么色好看| 亚洲精品第二区| 日韩三级伦理在线观看| 国产在线免费精品| 亚洲精品成人av观看孕妇| 桃花免费在线播放| 日本黄色日本黄色录像| 日韩 亚洲 欧美在线| 亚洲av.av天堂| 国产精品人妻久久久久久| 一本大道久久a久久精品| 少妇被粗大猛烈的视频| 国产不卡av网站在线观看| 2021少妇久久久久久久久久久| 26uuu在线亚洲综合色| 亚洲精华国产精华液的使用体验| 一本大道久久a久久精品| 久久精品国产自在天天线| 亚洲精华国产精华液的使用体验| 日本欧美国产在线视频| 人人妻人人澡人人看| 人妻一区二区av| 亚洲国产精品成人久久小说|