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

    水下溢油數(shù)值模擬研究

    2015-11-22 05:30:50陳海波尤云祥雷方輝趙宇鵬李建偉
    海洋工程 2015年2期
    關鍵詞:羽流單元體油滴

    陳海波,安 偉,楊 勇,尤云祥,雷方輝,趙宇鵬,李建偉

    (1.上海交通大學 海洋工程國家重點實驗室,上海 200240;2.中海石油環(huán)保服務(天津)有限公司,天津300452;3.中海油研究總院,北京 100027)

    我國南海擁有豐富的油氣資源,石油地質儲量約為230 ~300 億噸,占我國油氣總資源量的三分之一,其中70%蘊藏于深海區(qū)域,但深水區(qū)域特殊的自然環(huán)境和復雜的油氣儲藏條件決定了深水油氣勘探開發(fā)具有高投入、高回報、高技術、高風險的特點[1]。油氣開采過程中的溢油事故時有發(fā)生,2010年美國墨西哥灣“深水地平線”溢油事故警示我們,油氣勘探開發(fā)的高風險性應該受到足夠的重視,同時開采深海油氣資源的國家和公司也應該更加重視深海溢油的預防和治理。通過水下溢油數(shù)值模擬研究加深對水下溢油行為和歸宿的認識,為水下溢油應急策略的制定提供技術支持,而一個好的水下溢油數(shù)值模型也有助于準確評估溢油對環(huán)境的影響。20 世紀70年代國際上就有一些學者開始致力于水下溢油數(shù)值模擬的研究,建立了一些溢油模型,例如淺水溢油模型McDougall[2]、Fannelop 和Sjoen[3]、Milgram[4]、Lee 和Cheung[5]、Yapa 等[6],深水溢油模型Johansen[7]的DeepBlow 模型、Zheng 等[8]的CDOG 模型。近幾年,國內(nèi)也有一些學者開展了水下溢油數(shù)值模擬研究,例如汪守東等[9]、王晶等[10]、高清軍等[11]、廖國祥等[12-14],亓俊良等[15]。但與國外研究相比,目前國內(nèi)還沒有同時結合實驗室數(shù)據(jù)和現(xiàn)場觀測數(shù)據(jù)所進行的更充分的模型驗證研究。

    綜合前人研究,將水下溢油過程分為三個階段:噴發(fā)階段、浮力羽流階段和對流擴散階段(見圖1)。其中,噴發(fā)階段在距離噴口幾米范圍內(nèi),由于溢油一般以一定的初始速度噴射而出,噴出的溢油會具有較高的動能(如果有氣體溢出,動能會更高),溢油會因此分裂成大量直徑不等的小油滴,這些小油滴和海水組成的油團在高動能狀態(tài)下會繼續(xù)運動一段距離。受到水環(huán)境的緩沖作用,溢油在運動幾米后其初始動能很快損失殆盡,繼而進入浮力羽流階段。在浮力羽流階段,溢油運動具有一定的“主動”性,其運動狀態(tài)有別于水環(huán)境的流動狀態(tài),水環(huán)境通過浮力作用和卷吸作用間接地對溢油運動施加影響,這使油團持續(xù)抬升,體積不斷擴大,其密度不斷趨近于水環(huán)境密度。在這一過程中油團也會將部分底層水夾帶至較淺的水域。當油團成長到一定程度,油團整體的浮力水平接近中性,其運動的主動性不再明顯,溢油過程進入對流擴散階段。在對流擴散階段,由于油團體積較大,油滴比較分散,大量的油滴在水環(huán)境的流動和浮力的直接作用下被動地漂移和擴散。水下溢油的這三個階段沒有明顯界線,并且前兩個階段也包含對流擴散過程,只是在前兩個階段中噴流/羽流過程占主導地位,而對流擴散過程相對較弱,可以忽略。前兩個階段都是將溢油當作一個整體研究,可以用羽流動力模型模擬,第三個階段需要用對流擴散模型進行模擬。需要指出的是,如果溢油發(fā)生在淺水區(qū),溢油過程可能只需用羽流動力模型模擬。本文基于水下溢油三階段理論框架,通過分析溢油在水下環(huán)境的行為過程和影響因素,建立了一個水下溢油數(shù)值模型,并結合實驗室水槽實驗和水下溢油現(xiàn)場實驗的觀測資料進行模型驗證。

    圖1 水下溢油過程示意[8]Fig.1 Sketch of underwater oil spill process [8]

    1 數(shù)值模型

    該溢油模型由兩個子模型組成,包括羽流動力模型和對流擴散模型,其中羽流動力模型模擬噴發(fā)階段和浮力羽流階段,將一定量的溢油視為一個整體油團,考慮油團與海水的相互作用;對流擴散模型模擬對流擴散階段,將溢油離散為一定數(shù)量的“油粒子”,每個粒子代表一個由大量油滴組成的集合,并具有一定的質量、體積、濃度、油滴直徑等屬性,這些油粒子在海流的作用下漂移擴散。

    1.1 羽流動力模型

    羽流動力模型采用Lagrange 積分法模擬溢油的噴發(fā)階段和浮力羽流階段。如圖2 所示,將溢油持續(xù)時間平均分為若干等份,時間步長為Δt,每一個時間步對應一小份溢油,將每一小份溢油看作一個圓柱形的控制單元體,控制單元體的半徑為b(m),厚度為為控制單元體的移動速度(m/s),控制單元體質量為m = ρπb2h(kg),ρ 為控制單元體的密度(kg/m3)。控制單元體的底面法線平行于移動速度假設相鄰兩控制單元體之間互不影響。某一時刻所有控制單元體的中心連線即為該時刻的溢油軌跡??刂茊卧w的行為滿足以下控制方程:

    (1)質量守恒方程

    控制單元體在運動過程中會因卷吸作用而發(fā)生質量變化,此過程滿足以下質量守恒方程:

    式中:ρa為水環(huán)境的密度(kg/m3),Qe為卷吸體積通量(m3/s)。

    圖2 羽流動力模型控制單元體示意[5]Fig.2 Sketch of control volume of plume dynamics model [5]

    (2)動量守恒方程

    溢油進入水環(huán)境后,會與周圍水體相互作用導致動量變化,根據(jù)動量守恒可得到控制單元體的動量變化滿足以下關系:

    (3)熱量、鹽度和濃度守恒方程

    溢油進入水環(huán)境后,會與周圍水體發(fā)生物質交換,由此而導致的控制單元體的溫度、鹽度和濃度變化可用以下守恒方程描述:

    式中:I 為一個標量,表示熱含量(I = CpT)、鹽度(I = S)或油濃度(I = C)。Cp為比熱容(J/kg·K),T 為溫度(K),水比熱容Cp= 4 216.3 J/kg·K,油的比熱容約為水的一半。

    (4)狀態(tài)方程

    狀態(tài)方程描述了控制單元體的密度與溫度、鹽度和濃度之間的關系。Neumann 和Pierson[16]給出了水密度的計算方法,Bobra 和Chung[17]則給出了多種油的密度計算公式。在本文的溢油模型中,控制單元體中的物質是油和水組成的混合物,控制單元體在成長過程中通過卷吸作用將周圍海水不斷納入其中,導致控制單元體的溫度、鹽度和濃度發(fā)生變化,從而產(chǎn)生密度的改變。本文溢油模型考慮控制單元體的溫度和鹽度變化對密度的影響,并采用Bemporad[18]狀態(tài)方程計算控制單元體的密度:

    式中:ρ0表示參考密度值(kg/m3);βT和βS分別表示熱量膨脹系數(shù)(℃-1)和溶質的膨脹系數(shù)(‰-1),本文分別取值為5 ×10-4℃-1和8 ×10-4‰-1;ΔT 和ΔS 分別表示溫度和鹽度的變化量。

    (5)卷吸

    溢油進入水環(huán)境后,卷吸作用是影響羽流形態(tài)和行為的主要因素。卷吸作用是指羽流表面與水體之間相互作用,使周圍海水進入羽流中。根據(jù)Frick[19]和Lee 和Cheung[5]的研究,卷吸作用可分為兩部分:剪切卷吸和強迫卷吸。其中,剪切卷吸是由溢油運動導致的油團與水環(huán)境之間產(chǎn)生的剪切應力引起的,即使水環(huán)境是靜止的,剪切卷吸依然存在;強迫卷吸是在流動的水環(huán)境中水流迫使海水通過羽流的迎風面進入油團內(nèi)部引起的,強迫卷吸存在的前提為水環(huán)境是流動的。據(jù)此,本文溢油模型將卷吸體積通量Qe分為兩部分計算:

    式中:Qs代表剪切卷吸體積通量(m3/s),Qf代表強迫卷吸體積通量(m3/s)。剪切卷吸體積通量(Qs)采用Fischer[20]提出的公式進行計算:

    式中:α 為卷吸系數(shù)。Schatzmann[21-22]基于大量噴射流的研究,由動能方程和動量方程推導出了卷吸系數(shù)的計算公式,Lee 和Cheung[5]通過將該計算公式與噴射流實驗數(shù)據(jù)進行對比,給出了卷吸系數(shù)的具體計算形式:

    另一方面,如果水環(huán)境是流動的,羽流與水環(huán)境之間還存在強迫卷吸,強迫卷吸體積通量(Qf)可表達為:

    式中:Qfx、Qfy和Qfz分別為強迫卷吸在x、y 和z 方向上的分量。Frick[19]認為強迫卷吸由增長項、圓筒項和彎曲項三部分組成。其中,增長項由溢油軌跡半徑的增大導致;圓筒項由羽流側面投影到水流上的面積導致;彎曲項由溢油軌跡的彎曲導致?;贔rick[19]的研究,Lee 和Cheung[5]推導出同向水流中三維溢油強迫卷吸體積通量的表達式:

    式中:ua、va、wa分別為水流在x、y、z 方向的分量;θ 為控制單元體的速度矢量在水平面上的投影與x 軸的夾角(圖2);為控制單元體在一個時間步內(nèi)的位移大小,Δx、Δy、Δz 分別表示控制單元的體位移矢量在x、y、z 方向的分量。式(10)~式(12)等號右端括號中的三項分別為增長項、圓筒項和彎曲項的具體表達式。

    1.2 對流擴散模型

    油團上浮到一定高度后,由于卷吸了大量海水,其密度與水環(huán)境的密度已經(jīng)非常接近,尤其在密度層化的水環(huán)境中,油團會夾帶著來自于底層密度較大的水體,上浮到一定高度后其密度會等于甚至超過周圍水體的密度,這使油團失去浮力的動力支持而無法繼續(xù)抬升。然而,實際溢油事故并非如此。實際上,當油團發(fā)展到一定程度,油團內(nèi)外的海水性質和湍流強度已經(jīng)相差很小,羽流動力過程不再占主導地位,而水環(huán)境中的對流擴散過程開始主導溢油的運動,因而浮力羽流動力模型不再適用。取而代之,利用對流擴散模型繼續(xù)模擬溢油的行為。對流擴散模型的控制方程(對流擴散方程)為:

    方程(13)可以用多種數(shù)值方法求解[24-25],如有限差分法和有限元法,數(shù)值方法的主要困難在于穩(wěn)定性問題,即數(shù)值離散過程中可能會引進與物理擴散無關的數(shù)值擴散,數(shù)值擴散很大時,會使模擬結果失真,不能描述真實的物理過程。另外,還有一些過程難以用對流擴散方程模擬,如比較典型的非Fick 擴散問題[26]。對于溢油模擬,更多的學者采用的是Lagrange 粒子追蹤法(也叫隨機走動法),該方法是將溢油離散為一定數(shù)量具有拉格朗日性質的“油粒子”,每個粒子具有一定的質量、體積、濃度等屬性,粒子的平流過程主要受控于海浪、海流等環(huán)境動力過程的影響,可用確定性方法模擬,海洋湍流引起的擴散過程具有隨機性,用隨機運動來模擬。通過統(tǒng)計所有粒子的位置,可以確定溢油在海洋環(huán)境中的時空分布。Lagrange 粒子追蹤法,很早就在海洋和大氣研究中有著應用[27],Johansen[28]和Elliott[29]首先將此方法應用于溢油模擬(用于溢油模擬時也稱油粒子法)。該方法不僅避免了數(shù)值離散方法本身帶來的數(shù)值擴散問題,同時還可以準確重現(xiàn)海上油膜的破碎分離現(xiàn)象,準確描述溢油的真實擴散過程。大量研究表明,將該方法用于溢油模擬是行之有效的[30-33]。

    本文溢油模型利用Lagrange 粒子追蹤法模擬溢油的對流擴散過程。具體地,油粒子在三維空間中的運動滿足以下關系:

    式中:R 為在區(qū)間 [-1,1] 內(nèi)均勻分布的隨機數(shù)。

    油滴的浮力速度wb由油滴的直徑和密度以及海水的粘度和密度決定,可以根據(jù)油滴直徑大小分兩種情況計算[32]。首先定義一個臨界直徑

    式中:μ 表示海水的動力粘度,量級為10-3Ns/m2;ρ 和ρa分別為油滴和海水的密度。當油滴直徑d <dc時,根據(jù)Stokes 定律計算浮力速度:

    當油滴直徑d ≥dc時,根據(jù)Reynolds 定律計算浮力速度:

    具體計算中,將油粒子所在的控制單元體在浮力羽流階段終結時的位置作為該油粒子在對流擴散階段的初始位置。對于浮力羽流階段終結的判斷指標,Yapa 等[32]在層化的水環(huán)境中進行溢油模擬時將油團與水環(huán)境的密度差,作為判斷是否終結羽流動力模型的依據(jù)。Dasanayaka 和Yapa[33]對比了多個終結標指標的模擬效果,包括油滴速度指標(油團的上浮速度等于單個油滴的浮力速度)、中性浮力水平指標(油團的密度等于水環(huán)境的密度)和零速度指標(油團上浮速度為零),研究結果表明,油滴速度指標優(yōu)于其他兩個指標,并且基于平均油滴速度指標的溢油量估計更為準確。然而,在本文現(xiàn)場實驗的模擬中,由于考慮了海水的層化情況,仍然采用中性浮力水平指標作為判斷結束浮力羽流階段的依據(jù)。

    2 模型驗證

    為了檢驗本文水下溢油模型的準確性,將溢油模型應用到若干個水下溢油場景的數(shù)值模擬,并將模擬結果與實驗數(shù)據(jù)進行對比。實驗數(shù)據(jù)包括7 組來自實驗室水槽模擬實驗的觀測數(shù)據(jù)和一組在現(xiàn)場實驗中獲取的觀測數(shù)據(jù),其中實驗室數(shù)據(jù)包含了層化的和未層化的水環(huán)境中的噴流軌跡數(shù)據(jù)。

    2.1 實驗室水槽實驗驗證

    首先通過數(shù)值實驗將模擬結果與實驗室水槽噴射實驗的觀測數(shù)據(jù)進行對比,對本文水下溢油模型的準確性進行檢驗。共進行7 個數(shù)值實驗,分別考察了7 個噴射場景,考慮了不同的噴射條件和水環(huán)境層化情況。為區(qū)分不同的噴射場景,首先定義如下三個參數(shù):

    式中:帶有下標0 的變量代表初始值,F(xiàn)r 為Froude 數(shù);g' = g (Δρ/ρ) ,Δρ = ρa-ρ;D 為噴嘴直徑;St 為水環(huán)境的層化數(shù);R 為噴射速度與水環(huán)境橫流速度之比。各噴射實驗的參數(shù)取值見表1。在二維噴射模擬實驗中(實驗編號1 ~實驗編號5),為便于描述且不失一般性,將水環(huán)境的流向定義為x 軸的正方向,即( ua,0,0 )。

    表1 實驗室實驗參數(shù)Tab.1 Laboratory experimental parameters

    (1)未層化水環(huán)境中二維溢油模擬實驗驗證

    數(shù)值實驗1、2 和3 考察了在未層化且存在水平橫流的水環(huán)境中的二維溢油軌跡,三個實驗的噴射方向均垂直向上。圖3 給出了三個實驗中羽流控制體中心點的移動軌跡,并分別與Hirst[34]的三組實驗室水槽觀測數(shù)據(jù)進行了對比,可以看出,本文溢油模型的模擬結果與實驗數(shù)據(jù)具有很好的一致性。圖3 表明,垂直向上噴射的溢油,在水平橫流的作用下會向水流的下游漂移,這使溢油軌跡逐漸向下游傾斜,參數(shù)R0越大,傾斜程度越大。實際上,在初始時刻,油團的動量方向是垂直向上的,隨后,油團通過卷吸作用從水環(huán)境中獲得水平方向的動量,使水平方向的移動速度不斷增加。同時,雖然浮力能夠為油團提供向上的動量,但由于卷吸進來了大量垂向速度為零的水,最終導致油團向上的移動速度不斷減小,這使羽流形狀出現(xiàn)彎曲,當羽流成長到一定程度,溢油在水平方向和垂直方向的漂移速度均不再有明顯變化,此后溢油軌跡趨近于一條斜線。

    (2)層化水環(huán)境中二維溢油模擬實驗驗證

    實驗4 和實驗5 考察了垂直向上噴出的溢油在層化水流中的移動軌跡。圖4 將數(shù)值模擬的結果與Hirst[34]的實驗室水槽觀測數(shù)據(jù)進行了對比,可以看出,模擬值與觀測值的一致性較好,模擬結果非常理想,這說明本文模型對于層化水環(huán)境中的溢油模擬也具有很好的適用性。

    圖3 未層化的水流環(huán)境中的溢油軌跡(數(shù)據(jù)源自Hirst[34])Fig.3 Trajectories for oil spilled in an unstratified flowing environment (data from Hirst[34])

    圖4 層化的水流環(huán)境中的溢油軌跡(數(shù)據(jù)源自Hirst[34])Fig.4 Trajectories for oil spilled in a stratified flowing environment (data from Hirst[34])

    (3)三維溢油模擬實驗驗證

    下面利用Doneker 和Jirka[35]的兩個實驗室水槽溢油實驗(未層化的水環(huán)境)觀測數(shù)據(jù)檢驗本文溢油模型的準確性。在這兩個實驗中,溢油都是沿水平方向噴出,且噴射方向與水流方向垂直。在數(shù)值模擬中,將水流方向定義為x 軸的正方向,溢油的初始噴射方向為y 軸的正方向,見圖5。

    圖5 三維溢油軌跡(數(shù)據(jù)源自Doneker 和Jirka[35])Fig.5 Three-dimensional trajectories for oil spills (data from Doneker and Jirka[35])

    利用本文模型模擬的溢油軌跡和水槽實驗的觀測數(shù)據(jù)如圖5 所示,其中圖5(a)為溢油軌跡在xy 平面上的投影,圖5(b)為溢油軌跡在xz 平面上的投影,圖5(c)為溢油軌跡的三維圖像。通過圖5 模擬結果與實驗數(shù)據(jù)的對比可以看出,總體而言,兩個實驗的模擬結果都比較理想。然而由圖5(b)可以發(fā)現(xiàn),在溢油的垂向抬升高度(僅由浮力作用決定)方面,實驗6 的模擬結果要明顯高于水槽實驗的觀測數(shù)據(jù),這可能是由水槽實驗的邊界效應導致。如果提高水流速度(抑或保持流速不變,減小溢油噴射速度,從而削弱邊界效應的影響),會使模擬結果與觀測數(shù)據(jù)更為吻合,這一點在圖5 中實驗7 的結果對比中得到了證實。

    2.2 現(xiàn)場實驗驗證

    以上利用實驗室水槽中的小型溢油實驗驗證了本文溢油模型中的羽流動力模型的準確性。如要繼續(xù)檢驗另外一個子模型——對流擴散模型的準確性,從而完成對整個水下溢油模型的驗證,還需要大型的現(xiàn)場溢油實驗的支持。IKU、NOFO、Esso Norge 和Norsk Hydro 于1995年在北海的挪威海岸附近進行了大型的水下溢油實驗,Rye 等[36]和Rye 和Brandvik[37]對這次實驗進行了詳細描述。表2 給出了該實驗的參數(shù)設置,實驗海區(qū)溫度、鹽度和密度的垂向結構以及當時的海流實測值分別見圖6 和表3。Rye 等[36]通過分析現(xiàn)場拍攝的水下溢油圖像指出,噴出的溢油首先以羽流的形式運動和變化,在上浮到一定高度后轉變?yōu)橛纱罅坑偷谓M成的油團在浮力和海流的作用下被動地漂移擴散,轉變區(qū)域在50 ~60 m 水深之間。根據(jù)聲吶記錄數(shù)據(jù)得到的溢油水下范圍在東南-西北方向上的投影如圖7 所示,其中不同深度的橫線對應不同的觀測時刻?,F(xiàn)場觀測還表明,溢油最早到達海面的時間為溢油開始10 min 之后。

    表2 現(xiàn)場實驗參數(shù)[36]Tab.2 Field experimental parameters [36]

    圖6 現(xiàn)場實驗海域的溫度、鹽度和密度的垂向結構[36]Fig.6 Vertical profiles of temperature,salinity and density of field experiment [36]

    表3 現(xiàn)場實驗的海流信息[36]Tab.3 Current velocity profile of field experiment [36]

    基于表2 的實驗參數(shù)以及表3 和圖6 中的海洋環(huán)境數(shù)據(jù),利用本文模型進行水下溢油模擬,其中,海流數(shù)據(jù)由表3 中的數(shù)值通過線性插值得到。油滴直徑分布是對流擴散模型的一個重要參量,本文基于Delvigne[38]和Delvigne 和Sweeney[39]的實驗室研究,假定油滴直徑服從數(shù)學期望為250 μm,標準差為75 μm的正態(tài)分布,油滴直徑的最小值取為10 μm,最大值則根據(jù)Rye 等[36]的記錄取為5 000 μm??刂茊卧w通過卷吸作用從水壞境中不斷吸收海水,其密度會隨之增大,但由于海水的密度隨著水深的減小而減小,當控制單元體上升到某一深度后,其密度會等于或大于其所在位置周圍的海水密度。如前所述,在本文模型設置中,當控制單元體密度大于或等于周圍海水密度時,終止羽流動力模型的計算,繼而利用對流擴散模型取而代之。在模擬油粒子的隨機運動時,將水平和垂向擴散系數(shù)在海面以下10 m 的深度內(nèi)分別取為常數(shù)0.05和0.001 m2/s,由10 m 深度至海底,其取值線性減?。?2]。

    利用本文模型模擬的溢油分布在東南-西北方向的投影如圖7 所示,投影的下半部分由羽流動力模型計算得到,上半部分由對流擴散模型計算得到。數(shù)值模擬結果表明,浮力羽流階段的最大高度大約在52.1 m深度,溢油上浮至海面的最短時間為11.4 min,這說明本文模型的模擬結果與現(xiàn)場溢油實驗的現(xiàn)場數(shù)據(jù)具有很好的一致性。盡管如此,由圖7 還可以發(fā)現(xiàn),數(shù)值模擬結果與觀測數(shù)據(jù)之間也存在一定的偏差:(1)在100 m深度的羽流寬度小于觀測值,據(jù)Rye 等[36]記載,在現(xiàn)場實驗過程中,海洋波浪使噴油裝置存在小幅度的上下顛簸,這會增大噴射流初期的卷吸作用,而在本文數(shù)值模擬中并未考慮這一點,從而導致羽流寬度的模擬值較小;(2)在溢油圖像的上半部分(0 ~50 m 深度),溢油主體位置的模擬結果與觀測數(shù)據(jù)之間存在一定偏移,其原因應該是,本文數(shù)值模擬所用的流場只是根據(jù)在4 個深度上測得的海流數(shù)據(jù)插值得到,因而在細節(jié)方面不夠準確,而海流數(shù)據(jù)的精度能直接影響到對流擴散模型的模擬精度,另一方面,在現(xiàn)場觀測時不能保證觀測斷面是嚴格地東南-西北方向的,存在一定的角度偏差,這也是導致本文模擬結果出現(xiàn)位置偏差的一個可能原因。綜合以上兩點考慮,本文模型對現(xiàn)場實驗的模擬結果還是令人滿意的,從而驗證了本文模型的合理性和準確性。

    圖7 溢油水下分布在東南-西北方向的投影Fig.7 Projection of underwater oil spill in southeast-northwest direction

    進一步結果分析表明,在浮力羽流階段終結后,油滴直徑的選取對數(shù)值模擬結果有較大影響。直徑較大的油滴由于具有較快的浮力速度,能更快地到達海面,其運動軌跡更接近于垂直方向,并且由于受水下橫流的影響時間較短,其到達海面的位置會更靠近溢油點;反之,直徑較小的油滴到達海面的時間更長,其運動軌跡更接近于水平方向,并且到達海面的位置離溢油點也更遠。

    3 結 語

    通過分析溢油在水下環(huán)境的行為過程和影響因素,基于Lagrange 積分法和Lagrange 粒子追蹤法建立了一個水下溢油數(shù)值模型。該溢油模型由兩個子模型組成,包括羽流動力模型和對流擴散模型,其中羽流動力模型采用Lagrange 積分法模擬噴發(fā)階段和浮力羽流階段,將一定量的溢油視為一個整體油團,考慮油團與海水間的相互作用;對流擴散模型采用Lagrange 粒子追蹤法模擬對流擴散階段,將溢油離散為一定數(shù)量的“油粒子”,每個粒子代表一個由大量油滴組成的集合,并具有一定的質量、體積、濃度、油滴直徑等屬性,這些油粒子在海水浮力和海流的共同作用下漂移擴散。

    通過數(shù)值實驗,結合實驗室水槽實驗和水下溢油現(xiàn)場實驗的觀測資料進行模型驗證。數(shù)值實驗考察了層化和未層化的水流環(huán)境中的水下溢油過程以及實際海洋中的水下溢油過程。通過對比模擬結果與觀測數(shù)據(jù),驗證了本文水下溢油模型的合理性和準確性。進一步分析表明,在數(shù)值模型中,羽流動力模型的準確模擬能為對流擴散模型提供一個可靠的源,進入對流擴散階段后,海流、海水的垂向密度結構和油滴的直徑分布是影響溢油運動和分布的主要因素。該模型可以為水下溢油應急計劃提供技術參考。

    [1]李清平.我國海洋深水油氣開發(fā)面臨的挑戰(zhàn)[J].中國海上油氣,2006,18(2):130-133.(LI Qingping.The situation and challenges for deepwater oil and gas exploration and exploitation in China[J].China Offshore Oil and Gas,2006,18(2):130-133.(in Chinese))

    [2]MCDOUGALL T J.Bubble plume in stratified environments[J].Journal of Fluid Mechanics,1978,85(4):655-672.

    [3]FANNELOP T K,SJOEN K.Hydrodynamics of underwater blowouts[J].Norwegian Maritime Research,1980,(4):17-33.

    [4]MILGRAM J H.Mean flow in round bubble plumes[J].Journal of Fluid Mechanism,1983,133:345-376.

    [5]LEE J H W,CHEUNG V.Generalized Lagragian model for buoyant jets in current[J].Journal of Environmental Engineering,1990,116(6):1085-1106.

    [6]YAPA P D,ZHENG L.Simulation of oil spills from under water accidents I:Model development[J].Journal of Hydraulic Research,1997,35(5):673-687.

    [7]JOHANSEN ?.DeepBlow-a Lagrangian plume model for deep water blowouts[J].Spill Science & Technology Bulletin,2000,6(2):103-111.

    [8]ZHENG L,YAPA P D,CHEN F H.A model for simulating deepwater oil and gas blowouts Part I:Theory and model formulation[J].Journal of Hydraulic Research,2003,41(4):339-351.

    [9]汪守東,沈永明.海底管線溢油數(shù)學模型研究[J].大連理工大學學報,2006,46(S1):191-197.(WANG Shoudong,SHEN Yongming.Study of mathematical model for oil spills from seabed pipeline[J].Journal of Dalian University of Technology,2006,46(S1):191-197.(in Chinese))

    [10]王晶,李志軍,GONCHAROV V K,等.渤海海底管線溢油污染預測模型[J].海洋環(huán)境科學,2007,26(1):10-13.(WANG Jing,LI Zhijun,GONCHAROV V K,et al.A forecast model on oil pollution of leaks from seabed pipeline in the Bohai Sea[J].Marine Environmental Science,2007,26(1):10-13.(in Chinese))

    [11]高清軍,褚云峰,林建國.海底管線溢油的數(shù)值模擬[J].大連海事大學學報,2007,33(S2):169-171.(GAO Qingjun,ZHU Yunfeng,LIN Jianguo.Mathematical simulation of oil spills from seabed pipeline[J].Journal of Dalian Maritime University,2007,33(S2):169-171.(in Chinese))

    [12]廖國祥,高振會.水下溢油事故污染物輸移擴散的數(shù)值模擬研究[J].海洋環(huán)境科學,2011,30 (4):578-582.(LIAO Guoxiang,GAO Zhenhui.Numerical simulation of pollutant transport and diffusion in underwater oil spill[J].Marine Environmental Science,2011,30(4):578-582.(in Chinese))

    [13]廖國祥,高振會,楊建強.深海溢油輸移擴散模型研究[J].海洋環(huán)境科學,2012,31(5):718-723.(LIAO Guoxiang,GAO Zhenhui,YANG Jianqiang.A model for simulating the transport and behavior of oil and gas released from deepwater[J].Marine Environmental Science,2012,31(5):718-723.(in Chinese))

    [14]廖國祥,楊建強,高振會,等.深海中噴灑分散劑后溢油運動軌跡預測[J].大連海事大學學報,2013,39(2):103-107.(LIAO Guoxiang,YANG Jianqiang,GAO Zhenhui,et al.Numerical prediction of oil spill trajectory after dispersant application in deepwater environment[J].Journal of Dalian Maritime University,2013,39(2):103-107.(in Chinese))

    [15]亓俊良,李建偉,安偉,等.深水區(qū)水下溢油行為及歸宿研究[J].海洋開發(fā)與管理,2013,8:77-84.(QI Junliang,LI Jianwei,AN Wei,et al.Study on behavior and fate of oil spilled from deepwater[J].Ocean Development and Management,2013,8:77-84.(in Chinese))

    [16]NEUMANN G,PIERSON W J Jr.Principles of physical oceanography[M].Prentice-Hall,Inc.,Englewood Cliffs N J,1996.

    [17]BOBRA M A,CHUNG P T.A catalogue of oil properties[R].Report EE-77,Consultchem,Ottawa,Canada,March,1986.

    [18]BEMPORAD G A.Simulation of round buoyant jet in stratified flowing environment[J].Journal of Hydraulic Engineering,1994,120(5):529-543.

    [19]FRICK W E.Non-empirical closure of the plume equations[J].Atmospheric Environment,1984,18(4):653-662.

    [20]FISCHER H B,LIST E J,KOH R C Y,et al.Mixing in inland and coastal waters[M].Academic Press,New York,1979.

    [21]SCHATZMANN M.Calculation of submerged thermal plumes discharged into air and water flows [C]//Proc.18th IAHR Congress,Int.Assoc.for Hydr.Res.,Delft,The Netherlands,1979,4:379-385.

    [22]SCHATZMANN M.Mathematical modelling of submerged discharges into coastal waters[C]//Proc.19th IAHR Congress,Int.Assoc.for Hydr.Res.,Delft,The Netherlands,1981,3:239-246.

    [23]ZHENG L,YAPA P D.Simulation of oil spills from underwater accidents II:model verification[J].Journal of Hydraulic Research,1998,36(1):117-134.

    [24]金梅兵.近岸溢油的全動力預測方法研究[J].海洋環(huán)境科學,1997,16(1):30-36.(JIN Meibing.Study on the method of the dynamical prediction of spill oil inshore[J].Marine Environmental Science,1997,16(1):30-36.(in Chinese))

    [25]GOMEZ S,IVORRA B,RAMOS A M.Optimization of a pumping ship trajectory to clean oil contamination in the open sea[J].Mathematical and Computer Modelling,2011,54:477-489.

    [26]劉偉峰,孫英蘭.海上溢油運動數(shù)值模擬方法的探討與改進[J].華東師范大學學報:自然科學版,2009,3:90-97.(LIU Weifeng,SUN Yinglan.Study and improvement of oil spill simulation methods[J].Journal of East China Normal University:Natural Science,2009,3:90-97.(in Chinese))

    [27]CSANADY G T.Turbulent diffusion in the environment[M].Dordrecht,Boston,1973.

    [28]JOHANSEN O.The Halten Bank experiment-observations and model studies of drift and fate of oil in the marine environment[C]//Proceedings of the 11th Arctic Marine Oil Spill Program (AMOP)Technical Seminar,Environment Canada,1984:18-36.

    [29]ELLIOTT A J.Shear diffusion and the spread of oil in the surface layers of the North Sea[J].Deutsche Hydrographische Zeitschrift,1986,39(3):113-137.

    [30]HUNTER J R.The application of Lagrangian particle-tracking technique to modelling of dispersant in the sea [C]//In Numerical Modelling:Application to Marine Systems,ed.Noye J,North-Holland,1987.

    [31]WANG S D,SHEN Y M,GUO Y K,et al.Three-dimensional numerical simulation for transport of oil spills in seas[J].Ocean Engineering,2008,35:503-510.

    [32]YAPA P D,ZHENG L,NAKATA K.Modeling underwater oil/gas jets and plumes[J].Journal of Hydraulic Engineering,1999,125(5):481-491.

    [33]DASSANAYAKA L K,YAPA P D.Role of plume dynamics phase in a deepwater oil and gas release model[J].Journal of Hydro-environment Research,2009,(2):243-253.

    [34]HIRST E.Buoyant jets with three-dimensional trajectories[J].Journal of the Hydraulics Division,1972,98(11):1999-2014.

    [35]DONEKER R L,JIRKA G H.CORMIXI:an expert system for hydrodynamic mixing zone analysis of conventional and toxic submerged single port discharge[R].Tech.Rep.EPA/600/3-90/012,U.S.Environmental Protection Agency,Athens,Ga,1990.

    [36]RYE H,BRANDVIK P J,REED M.Subsurface oil release field experiment-observations and modelling of subsurface plume behavior[C]//Proc.19th Arctic and Marine Oil Spill Program (AMOP)Tech.Seminar,Vol.2,Environment Canada,Ottawa,1996:1417-1435.

    [37]RYE H,BRANDVIK P J.Verification of subsurface oil spill models[C]//Proc.1997 Int.Oil Spill Conf.,Environment Canada,Ottawa,1997:551-577.

    [38]DELVIGNE G A L.Natural and chemical dispersion of oil[J].Journal of Advanced Marine Technology Conference of Japan,1994,11:23-40.

    [39]DELVIGNE G A L,SWEENEY C E.Natural dispersion of oil[J].Oil and Chemical Pollution,1989,4:281-310.

    猜你喜歡
    羽流單元體油滴
    超高層單元體吊裝技術及安裝施工方法研究
    建筑與裝飾(2024年1期)2024-01-25 08:47:56
    圓形的油滴
    小主人報(2022年18期)2022-11-17 02:19:56
    水下羽流追蹤方法研究進展
    球墨鑄鐵復合仿生耦合單元體結構參數(shù)變化對摩擦應力的影響模擬研究
    某渦軸發(fā)動機單元體設計分析
    密里根油滴實驗的理論分析和測量結果討論
    水下管道向下泄漏的羽/射流特性
    化工學報(2016年12期)2016-12-14 09:28:02
    典型民用航空發(fā)動機單元體劃分淺析
    烷烴油滴在超臨界二氧化碳中溶解的分子動力學模擬
    軸承腔潤滑油沉積特征分析
    五月伊人婷婷丁香| av在线蜜桃| 日本色播在线视频| 舔av片在线| 日韩免费高清中文字幕av| 免费黄网站久久成人精品| 精品久久久久久久人妻蜜臀av| av天堂中文字幕网| 我的女老师完整版在线观看| 最近手机中文字幕大全| 国产日韩欧美在线精品| .国产精品久久| 日韩伦理黄色片| 国国产精品蜜臀av免费| 久久久久网色| 人妻 亚洲 视频| 全区人妻精品视频| 国产成年人精品一区二区| 国产精品久久久久久av不卡| 亚洲av不卡在线观看| 人妻一区二区av| 日韩欧美 国产精品| 免费不卡的大黄色大毛片视频在线观看| 亚洲高清免费不卡视频| 亚洲天堂国产精品一区在线| 日本免费在线观看一区| 国产欧美另类精品又又久久亚洲欧美| 一区二区三区乱码不卡18| 欧美成人午夜免费资源| 美女脱内裤让男人舔精品视频| 亚洲美女视频黄频| 日韩免费高清中文字幕av| av在线app专区| 三级经典国产精品| 日韩大片免费观看网站| 18禁裸乳无遮挡免费网站照片| 别揉我奶头 嗯啊视频| 久久久久国产网址| 一区二区av电影网| 婷婷色综合www| 高清在线视频一区二区三区| 日本-黄色视频高清免费观看| 亚洲av男天堂| 一级毛片我不卡| 尤物成人国产欧美一区二区三区| 夫妻性生交免费视频一级片| 亚洲自偷自拍三级| 狂野欧美白嫩少妇大欣赏| 秋霞伦理黄片| 自拍偷自拍亚洲精品老妇| 舔av片在线| 99热这里只有是精品在线观看| 午夜视频国产福利| 欧美日韩视频精品一区| 色播亚洲综合网| 在线亚洲精品国产二区图片欧美 | 男女啪啪激烈高潮av片| 欧美国产精品一级二级三级 | 王馨瑶露胸无遮挡在线观看| 十八禁网站网址无遮挡 | 欧美xxⅹ黑人| 国产亚洲午夜精品一区二区久久 | 亚洲熟女精品中文字幕| 毛片一级片免费看久久久久| 国产色爽女视频免费观看| 日本免费在线观看一区| 亚洲内射少妇av| 亚洲va在线va天堂va国产| 中文字幕人妻熟人妻熟丝袜美| 美女脱内裤让男人舔精品视频| 一级毛片 在线播放| 日韩欧美精品免费久久| 精品视频人人做人人爽| 亚洲精品乱码久久久久久按摩| 国产成人精品婷婷| 亚洲国产精品成人久久小说| 在线播放无遮挡| 91aial.com中文字幕在线观看| 色视频在线一区二区三区| 超碰97精品在线观看| 在线精品无人区一区二区三 | 亚洲欧美日韩卡通动漫| 国产91av在线免费观看| 熟女电影av网| 久久久精品94久久精品| 国产精品一二三区在线看| 亚洲,欧美,日韩| 日本午夜av视频| 亚洲,一卡二卡三卡| 亚洲人成网站在线播| 久久久久久久久久久免费av| 国产精品秋霞免费鲁丝片| 能在线免费看毛片的网站| 日日啪夜夜爽| 日韩av不卡免费在线播放| 91久久精品国产一区二区三区| 麻豆精品久久久久久蜜桃| 人妻一区二区av| 3wmmmm亚洲av在线观看| 亚洲欧美清纯卡通| 免费人成在线观看视频色| 亚洲精品亚洲一区二区| 久久97久久精品| 亚洲精品乱码久久久v下载方式| 夫妻性生交免费视频一级片| 夫妻午夜视频| 如何舔出高潮| 午夜福利在线观看免费完整高清在| 国产一区亚洲一区在线观看| a级毛色黄片| 免费大片黄手机在线观看| 日韩不卡一区二区三区视频在线| 久久久久九九精品影院| 久久久欧美国产精品| 街头女战士在线观看网站| 最近中文字幕高清免费大全6| 国产91av在线免费观看| 成人毛片a级毛片在线播放| 日韩强制内射视频| 别揉我奶头 嗯啊视频| 欧美三级亚洲精品| 听说在线观看完整版免费高清| 女人被狂操c到高潮| 99久久九九国产精品国产免费| 精品久久久久久电影网| 免费观看的影片在线观看| 亚洲丝袜综合中文字幕| 婷婷色综合大香蕉| 成年人午夜在线观看视频| 小蜜桃在线观看免费完整版高清| 99热国产这里只有精品6| 国语对白做爰xxxⅹ性视频网站| 夜夜爽夜夜爽视频| 亚洲精品久久久久久婷婷小说| 亚洲av一区综合| 成年av动漫网址| 国产精品久久久久久久电影| 国国产精品蜜臀av免费| 国产在线一区二区三区精| 国产免费福利视频在线观看| 色视频在线一区二区三区| 国产乱人偷精品视频| 身体一侧抽搐| 国产高清有码在线观看视频| 99热这里只有是精品50| 日韩av不卡免费在线播放| av.在线天堂| 国产黄a三级三级三级人| 国产乱来视频区| 看免费成人av毛片| 国产精品三级大全| 又爽又黄无遮挡网站| 久久女婷五月综合色啪小说 | 别揉我奶头 嗯啊视频| 麻豆久久精品国产亚洲av| 欧美性猛交╳xxx乱大交人| 在线亚洲精品国产二区图片欧美 | 三级国产精品片| 久久鲁丝午夜福利片| 七月丁香在线播放| 最近2019中文字幕mv第一页| 国产精品一区二区在线观看99| 日韩大片免费观看网站| 亚洲精品456在线播放app| 欧美日韩综合久久久久久| 日韩av免费高清视频| 国产亚洲一区二区精品| 国产免费福利视频在线观看| 天美传媒精品一区二区| 真实男女啪啪啪动态图| 免费观看a级毛片全部| 男插女下体视频免费在线播放| 国产成人精品一,二区| 欧美精品人与动牲交sv欧美| 美女内射精品一级片tv| 一级毛片aaaaaa免费看小| 国产精品成人在线| 99久国产av精品国产电影| 欧美日本视频| 黑人高潮一二区| 麻豆精品久久久久久蜜桃| 一个人看的www免费观看视频| 一级片'在线观看视频| 免费看日本二区| 免费观看在线日韩| 嫩草影院新地址| 国产 一区 欧美 日韩| 日韩av在线免费看完整版不卡| 欧美性感艳星| 青春草视频在线免费观看| 亚洲精品日韩av片在线观看| 毛片女人毛片| 黄色一级大片看看| 综合色av麻豆| av线在线观看网站| 一区二区三区精品91| 欧美成人一区二区免费高清观看| 人妻一区二区av| 97超视频在线观看视频| 18禁动态无遮挡网站| 永久网站在线| 一级毛片电影观看| 男女那种视频在线观看| 看非洲黑人一级黄片| 亚洲精品日本国产第一区| 亚洲在久久综合| 免费不卡的大黄色大毛片视频在线观看| 99热这里只有是精品在线观看| 国产av码专区亚洲av| 激情 狠狠 欧美| 晚上一个人看的免费电影| 身体一侧抽搐| 日韩欧美精品免费久久| 国产一区二区三区综合在线观看 | 黄色视频在线播放观看不卡| 白带黄色成豆腐渣| av在线天堂中文字幕| 97在线人人人人妻| 日产精品乱码卡一卡2卡三| 亚洲va在线va天堂va国产| 亚洲国产精品999| 亚洲一级一片aⅴ在线观看| 热re99久久精品国产66热6| 人体艺术视频欧美日本| 亚洲色图综合在线观看| 91精品一卡2卡3卡4卡| 免费av毛片视频| 建设人人有责人人尽责人人享有的 | 国产大屁股一区二区在线视频| 亚洲真实伦在线观看| 禁无遮挡网站| h日本视频在线播放| 男女国产视频网站| 国产探花在线观看一区二区| 国产精品一及| 婷婷色av中文字幕| 国产精品一区二区性色av| 日日摸夜夜添夜夜爱| 午夜免费观看性视频| 少妇人妻 视频| 亚洲精品成人av观看孕妇| 99精国产麻豆久久婷婷| 国产免费视频播放在线视频| 亚洲无线观看免费| 欧美另类一区| 日韩电影二区| 亚洲av免费在线观看| 久热久热在线精品观看| 午夜日本视频在线| 日韩欧美精品v在线| 日本一二三区视频观看| 久久久久久久久久人人人人人人| 青春草亚洲视频在线观看| 真实男女啪啪啪动态图| 少妇人妻久久综合中文| 热re99久久精品国产66热6| 我的老师免费观看完整版| 免费观看无遮挡的男女| 人人妻人人爽人人添夜夜欢视频 | 亚洲av欧美aⅴ国产| 六月丁香七月| 女人久久www免费人成看片| 国产视频首页在线观看| 可以在线观看毛片的网站| 精品久久久久久电影网| 蜜桃久久精品国产亚洲av| 亚洲精品视频女| 国语对白做爰xxxⅹ性视频网站| 一个人看视频在线观看www免费| 男女边吃奶边做爰视频| 日日啪夜夜爽| 亚洲国产av新网站| 99热这里只有精品一区| 国产黄色免费在线视频| 久久精品夜色国产| 国产乱人视频| 99久久精品国产国产毛片| 国产成人精品福利久久| 一区二区三区四区激情视频| 日韩在线高清观看一区二区三区| 久久人人爽人人爽人人片va| 久久久久久久亚洲中文字幕| 亚洲综合色惰| 国产乱人偷精品视频| 亚洲精品日韩av片在线观看| 女人被狂操c到高潮| 男女啪啪激烈高潮av片| 国产 精品1| 在线天堂最新版资源| 国产成人免费无遮挡视频| 免费av观看视频| xxx大片免费视频| 日日啪夜夜爽| www.色视频.com| 免费看日本二区| 久久久午夜欧美精品| 国产精品嫩草影院av在线观看| 三级经典国产精品| 一个人观看的视频www高清免费观看| 中文资源天堂在线| 综合色av麻豆| 免费看不卡的av| 在线观看三级黄色| 国产日韩欧美在线精品| 五月玫瑰六月丁香| 国产精品麻豆人妻色哟哟久久| 又大又黄又爽视频免费| 国产精品伦人一区二区| 亚洲成人精品中文字幕电影| 国产成人aa在线观看| 日本黄大片高清| 亚洲va在线va天堂va国产| 成人免费观看视频高清| av卡一久久| 伊人久久精品亚洲午夜| 精品少妇久久久久久888优播| 男女无遮挡免费网站观看| 最近最新中文字幕免费大全7| 中文字幕免费在线视频6| 欧美日韩视频精品一区| 免费看日本二区| 只有这里有精品99| 性色av一级| 蜜臀久久99精品久久宅男| 国产 精品1| 我的老师免费观看完整版| 美女cb高潮喷水在线观看| 亚洲av二区三区四区| 日日摸夜夜添夜夜爱| 免费看光身美女| 久久人人爽av亚洲精品天堂 | 只有这里有精品99| 热re99久久精品国产66热6| 汤姆久久久久久久影院中文字幕| a级毛色黄片| 亚洲av国产av综合av卡| 51国产日韩欧美| 边亲边吃奶的免费视频| 熟女人妻精品中文字幕| 亚洲精品乱码久久久久久按摩| 下体分泌物呈黄色| 亚洲激情五月婷婷啪啪| 中文乱码字字幕精品一区二区三区| 精品国产一区二区三区久久久樱花 | 搡老乐熟女国产| 少妇的逼水好多| 亚洲内射少妇av| 大又大粗又爽又黄少妇毛片口| 国产欧美另类精品又又久久亚洲欧美| 夜夜爽夜夜爽视频| 亚洲国产欧美人成| 国产在线男女| 亚洲最大成人av| 赤兔流量卡办理| 亚洲国产欧美人成| 亚洲国产av新网站| 精品99又大又爽又粗少妇毛片| av在线观看视频网站免费| 熟女人妻精品中文字幕| 国产黄色免费在线视频| 久久精品久久久久久久性| freevideosex欧美| 国产视频内射| 日本色播在线视频| 亚洲久久久久久中文字幕| 日韩中字成人| 蜜桃久久精品国产亚洲av| 最新中文字幕久久久久| 丝瓜视频免费看黄片| 国产伦精品一区二区三区视频9| 精品99又大又爽又粗少妇毛片| 26uuu在线亚洲综合色| 特大巨黑吊av在线直播| 国产精品av视频在线免费观看| 成人鲁丝片一二三区免费| 免费高清在线观看视频在线观看| 干丝袜人妻中文字幕| 亚洲激情五月婷婷啪啪| 韩国高清视频一区二区三区| 亚洲成人精品中文字幕电影| 国产精品久久久久久精品古装| 街头女战士在线观看网站| 又爽又黄无遮挡网站| 亚洲天堂av无毛| 美女视频免费永久观看网站| 少妇人妻 视频| 亚洲国产精品国产精品| 日韩三级伦理在线观看| 99久久精品一区二区三区| 亚洲欧洲日产国产| 国产伦理片在线播放av一区| 亚洲色图综合在线观看| 国产高清不卡午夜福利| 色婷婷久久久亚洲欧美| 亚洲精品亚洲一区二区| 男人狂女人下面高潮的视频| 少妇人妻 视频| 97在线人人人人妻| 啦啦啦在线观看免费高清www| 午夜精品一区二区三区免费看| 精品午夜福利在线看| 亚洲三级黄色毛片| 成人特级av手机在线观看| 日韩在线高清观看一区二区三区| 69人妻影院| 大香蕉久久网| 丰满少妇做爰视频| 午夜福利网站1000一区二区三区| 午夜激情福利司机影院| av免费观看日本| 国产精品一二三区在线看| 亚洲成人一二三区av| 亚洲精品成人久久久久久| 免费黄网站久久成人精品| 亚洲欧美成人综合另类久久久| 久久久久久久久久久丰满| 丰满人妻一区二区三区视频av| av免费在线看不卡| 欧美成人午夜免费资源| 亚洲伊人久久精品综合| 日韩三级伦理在线观看| 久久久精品欧美日韩精品| 小蜜桃在线观看免费完整版高清| 在线观看国产h片| 成人综合一区亚洲| 一级毛片电影观看| 免费看光身美女| 国产精品久久久久久av不卡| 天美传媒精品一区二区| av在线亚洲专区| 欧美zozozo另类| 精品国产一区二区三区久久久樱花 | 久久精品国产自在天天线| 国产乱人偷精品视频| 黑人高潮一二区| 午夜免费观看性视频| 一级二级三级毛片免费看| 久久久久国产网址| 午夜老司机福利剧场| 制服丝袜香蕉在线| a级毛色黄片| 日韩伦理黄色片| 国产亚洲av片在线观看秒播厂| 777米奇影视久久| 99热这里只有是精品50| 久久精品综合一区二区三区| 国产成人freesex在线| 欧美成人a在线观看| 国产熟女欧美一区二区| 亚洲成人久久爱视频| 国产精品偷伦视频观看了| 精品99又大又爽又粗少妇毛片| 99久久中文字幕三级久久日本| 精品视频人人做人人爽| 色视频在线一区二区三区| 亚洲成色77777| 男人和女人高潮做爰伦理| 久久精品人妻少妇| 特大巨黑吊av在线直播| 亚洲av免费在线观看| 久久久精品免费免费高清| 国产精品女同一区二区软件| h日本视频在线播放| 一本一本综合久久| 汤姆久久久久久久影院中文字幕| 成年av动漫网址| 国产成人91sexporn| 午夜精品一区二区三区免费看| 黑人高潮一二区| 乱码一卡2卡4卡精品| 高清av免费在线| kizo精华| 亚洲精品影视一区二区三区av| 寂寞人妻少妇视频99o| 99热全是精品| 国产精品av视频在线免费观看| 在线观看人妻少妇| 日韩av免费高清视频| 久久久欧美国产精品| 亚洲av一区综合| 天天一区二区日本电影三级| 男男h啪啪无遮挡| 亚洲成人av在线免费| 日韩免费高清中文字幕av| 国内精品美女久久久久久| 下体分泌物呈黄色| 国产精品福利在线免费观看| 狠狠精品人妻久久久久久综合| 亚洲美女视频黄频| 在线观看一区二区三区激情| 丰满乱子伦码专区| www.av在线官网国产| 一区二区三区免费毛片| 国产精品一区二区性色av| 人妻系列 视频| 少妇人妻一区二区三区视频| 亚洲久久久久久中文字幕| 国产精品一及| 国产探花在线观看一区二区| 在线观看国产h片| 久久久久国产网址| 大香蕉97超碰在线| 久久久亚洲精品成人影院| 热re99久久精品国产66热6| 最新中文字幕久久久久| 日本一本二区三区精品| 嫩草影院精品99| 国产成人freesex在线| 国产精品久久久久久av不卡| 亚洲欧美一区二区三区黑人 | 成年女人看的毛片在线观看| 色婷婷久久久亚洲欧美| 欧美亚洲 丝袜 人妻 在线| 欧美成人精品欧美一级黄| 国产伦精品一区二区三区视频9| 亚洲人成网站高清观看| 国产视频首页在线观看| 国产毛片a区久久久久| 一区二区三区免费毛片| 国产女主播在线喷水免费视频网站| 久久精品夜色国产| 3wmmmm亚洲av在线观看| 精品久久久噜噜| 亚洲精品国产成人久久av| 天天一区二区日本电影三级| 午夜福利视频1000在线观看| 午夜日本视频在线| 如何舔出高潮| 亚洲av福利一区| 97超碰精品成人国产| 亚洲色图综合在线观看| 欧美三级亚洲精品| 日韩一区二区三区影片| 国内揄拍国产精品人妻在线| 国产日韩欧美在线精品| 一级二级三级毛片免费看| 国产白丝娇喘喷水9色精品| 国内精品宾馆在线| 欧美亚洲 丝袜 人妻 在线| 成年av动漫网址| 国产男女内射视频| 国产探花极品一区二区| 欧美日韩国产mv在线观看视频 | 久久精品久久久久久噜噜老黄| 国产 一区 欧美 日韩| 激情 狠狠 欧美| 99视频精品全部免费 在线| 亚洲精品成人av观看孕妇| www.色视频.com| 91精品一卡2卡3卡4卡| 亚洲va在线va天堂va国产| 国产黄频视频在线观看| 真实男女啪啪啪动态图| 天美传媒精品一区二区| 国产女主播在线喷水免费视频网站| 在线亚洲精品国产二区图片欧美 | 久久久久国产网址| 人妻少妇偷人精品九色| 精品国产三级普通话版| 小蜜桃在线观看免费完整版高清| 国模一区二区三区四区视频| 亚洲精品久久久久久婷婷小说| 中文字幕人妻熟人妻熟丝袜美| 美女xxoo啪啪120秒动态图| 尾随美女入室| 亚洲av电影在线观看一区二区三区 | 夜夜看夜夜爽夜夜摸| 日韩免费高清中文字幕av| 亚洲电影在线观看av| 毛片女人毛片| 高清午夜精品一区二区三区| 狂野欧美激情性xxxx在线观看| 少妇的逼水好多| 色视频在线一区二区三区| 黄色配什么色好看| 亚洲欧美日韩卡通动漫| 伊人久久精品亚洲午夜| 欧美性感艳星| 熟女av电影| av在线天堂中文字幕| 免费av不卡在线播放| 久久精品国产亚洲网站| 国产精品一及| 不卡视频在线观看欧美| 国产 精品1| 久久久成人免费电影| 91狼人影院| 亚洲精品日本国产第一区| 成人国产麻豆网| 亚洲精华国产精华液的使用体验| 午夜福利视频1000在线观看| 五月伊人婷婷丁香| 另类亚洲欧美激情| 免费观看的影片在线观看| eeuss影院久久| 亚洲av中文字字幕乱码综合| 18禁裸乳无遮挡免费网站照片| 亚洲精品日韩av片在线观看| av专区在线播放| 日韩欧美精品免费久久| 国产 一区精品| 欧美潮喷喷水| 成人无遮挡网站| 黄片无遮挡物在线观看| 欧美潮喷喷水| 亚洲va在线va天堂va国产| av国产免费在线观看| 美女被艹到高潮喷水动态| 国模一区二区三区四区视频| 国产 一区精品| 亚洲精品乱久久久久久| 亚洲欧美成人精品一区二区| 免费观看的影片在线观看| 最新中文字幕久久久久| 亚洲国产精品999|