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

    蓬萊19-3溢油運移與歸宿數(shù)值模擬和分析

    2016-11-23 06:19:50匡翠萍謝華浪
    關(guān)鍵詞:蓬萊溢油擴散系數(shù)

    匡翠萍, 謝華浪, 蘇 平, 陳 括

    (1. 同濟大學(xué) 土木工程學(xué)院, 上海 200092; 2. 廣西壯族自治區(qū)交通規(guī)劃勘察設(shè)計研究院, 廣西 南寧 530029)

    ?

    蓬萊19-3溢油運移與歸宿數(shù)值模擬和分析

    匡翠萍1, 謝華浪1, 蘇 平2, 陳 括1

    (1. 同濟大學(xué) 土木工程學(xué)院, 上海 200092; 2. 廣西壯族自治區(qū)交通規(guī)劃勘察設(shè)計研究院, 廣西 南寧 530029)

    通過耦合水動力模型建立渤海海上溢油運移與歸宿數(shù)學(xué)模型,包含漂移、擴散、蒸發(fā)、乳化等過程,通過拉格朗日方法模擬漂移過程、蒙特卡羅方法模擬油粒子擴散的隨機性,采用變化的紊動擴散系數(shù)使溢油所受紊動擴散作用隨時間變化.通過率定得到蓬萊19-3溢油模擬中風(fēng)拖拽系數(shù)和擴散系數(shù)校正參數(shù)b的合理取值范圍分別為1.4%~2.3%和0.38~0.42.對蓬萊19-3溢油事故進行了模擬,通過溢油附岸位置與時間驗證了模型的準(zhǔn)確性.模擬結(jié)果表明: 由于風(fēng)速較大、方向變化較小,風(fēng)對漂浮于水面上的溢油產(chǎn)生較大作用,且溢油漂移方向與主風(fēng)向一致;渤海海域內(nèi)的潮流以往復(fù)流、旋轉(zhuǎn)流為主且余流小使其對溢油運移的作用弱于風(fēng);溢油初期的蒸發(fā)速度較快,后期因揮發(fā)組分的減少以及乳化作用使溢油含水率增大使得蒸發(fā)速度明顯減?。徽舭l(fā)和附岸是溢油漂浮組分減少的主要原因,也是溢油的主要歸宿.

    渤海; 海上溢油; 運移; 數(shù)值模擬

    目前,海上石油污染已成為一種全球性海洋污染[1].海上溢油以其突發(fā)性、污染的嚴(yán)重性以及應(yīng)急處理的艱巨性使其在過去的幾十年里成為海洋生態(tài)環(huán)境最大的破壞因素,且呈現(xiàn)越來越嚴(yán)重的趨勢[2].尤其當(dāng)溢油事故發(fā)生于近岸或河口水體環(huán)境時,被污染的水體及近岸環(huán)境將受到慘重的破壞[3].2011年渤海蓬萊19-3溢油事故造成了極其惡劣的環(huán)境破壞,海水污染面積高達4 900 km2,造成的海洋生態(tài)損害價值總計16.83億元[4].根據(jù)康菲公司公布數(shù)據(jù),本次事故溢油量為500 t,然而多數(shù)研究學(xué)者認(rèn)為實際溢油量甚至高達該數(shù)據(jù)的十倍以上[5].

    海上溢油屬于突發(fā)性的海洋污染事故,且傳播速度快、影響范圍大.但由于海洋環(huán)境復(fù)雜多變,溢油事故發(fā)生于不同的時間和地點都可能產(chǎn)生不同的漂移軌跡和污染范圍[6].因此,對海上溢油進行數(shù)值模擬,有利于預(yù)測溢油變化趨勢,制定合理的治理方案,及時采取應(yīng)對措施,最大程度地減少溢油事故帶來的破壞和損失.徐青等[7]通過GNOME導(dǎo)入潮流和風(fēng)建立了渤海溢油模型,通過將溢油油膜概化為由多數(shù)圓形斑點組成,并采用拉格朗日方法模擬了渤海蓬萊19-3溢油事故6月11—19日中部分油膜的漂移過程.鄧增安等[3-6]考慮了風(fēng)、潮流、波浪等因素,研發(fā)了嵌套于中國數(shù)字海洋原型系統(tǒng)(CDOPS)中的BSOSSM溢油模擬服務(wù)組件,并對蓬萊19-3溢油事故進行了模擬.汪守東[8-9]與王金華等[10]將海上溢油分為表層浮油和懸浮油滴,通過離散顆粒的拉格朗日算法進行計算,考慮了溢油的漂移、擴散、蒸發(fā)和乳化等過程,建立了渤海海域溢油模型.畢海普等[11]建立了三峽水庫溢油模型,通過在模型中設(shè)置一個衰減參數(shù)來模擬與分析三峽水庫溢油事故發(fā)生后各溢油組分在機械回收和投放分散劑、乳化劑等緊急措施影響下的變化情況.李燕等[12]建立了三維的拉格朗日隨機步溢油模型,并研究了在不同方案下油粒子垂直擴散的變化情況.總體來說,目前各種溢油模型所采用方法與考慮因素各有不同.其中,對溢油的漂移過程模擬均采用拉格朗日方法,計算中主要將風(fēng)和潮流的作用作為溢油漂移的主要影響因素.半封閉的渤海海域中風(fēng)對溢油漂移的影響是最重要的,一般認(rèn)為風(fēng)拖拽系數(shù)取值范圍在1%~3%之間,但對渤海區(qū)域風(fēng)拖拽系數(shù)的取值尚未有詳細(xì)的討論和研究.另外,多數(shù)溢油模型由于對模擬粒子數(shù)限制較大,以至于對大規(guī)模的溢油模擬未能形成一個較為直觀的溢油平面密度分布場.

    利用Delft3D通過對水動力模型進行耦合建立了渤海海上溢油模型,模型計算包括漂移、擴散、蒸發(fā)、乳化等變化過程.本模型將溢油作為離散化的油粒子,通過拉格朗日方法模擬漂移過程、蒙特卡羅方法模擬湍流擴散過程,且采用隨時間變化的紊動擴散系數(shù)使溢油所受湍流作用亦隨時間變化,使得大規(guī)模的渦流和環(huán)流對溢油的擴散作用得以體現(xiàn)出來,并采用大量的油粒子進行模擬使溢油分布密度具有較高分辨率,對于判斷不同區(qū)域的污染程度差異更為清晰直觀.通過模型率定得到了風(fēng)拖拽系數(shù)和擴散系數(shù)校正參數(shù)b的合理取值范圍,在此基礎(chǔ)上對2011年蓬萊19-3溢油事故進行了驗證模擬,研究了渤海海域風(fēng)和潮流對溢油運移的影響,并分析了溢油歸宿的變化特征.

    1 數(shù)學(xué)模型

    1.1 Delft3D-FLOW水動力模型

    FLOW模型的基本方程為二維或三維的非線性淺水方程,采用淺水假定和布西奈斯克(Boussinesq)近似由不可壓縮流體的三維Navier-Stokes方程推求得出.由于垂向加速度遠遠小于重力加速度,故在垂直方向上的動量方程中忽略了垂直加速度的影響,并由此推導(dǎo)出在靜壓假定下的水流方程[13].平面二維水動力模型的基本方程為

    (1) 連續(xù)方程

    (1)

    式中:ζ為水位;d為基準(zhǔn)面下的深度,總水深H=d+ζ;u和ν分別表示水平方向上ζ和η向的流速分量;Gξξ和Gηη分別為正交曲線坐標(biāo)系和Cartesian坐標(biāo)系之間的轉(zhuǎn)換系數(shù);Q為單位面積的源匯項.

    (2) 動量方程

    ξ方向的動量方程

    (2)

    η方向的動量方程

    (3)式中:ξ和η為平面曲線坐標(biāo);f=2Ωsinφ,為科氏力參數(shù),由地球自轉(zhuǎn)角速度和緯度值決定;ρ0為水的密度;Pξ和Pη為壓力梯度;Fξ和Fη為水平雷諾應(yīng)力的不平衡性;Mξ和Mη為外來源匯項引起的附加動量.

    1.2 Delft3D-PART模塊

    Delft3D-PART模型是一個包含隨機運動的粒子追蹤模型,它將溶解或懸浮于水中的物質(zhì)通過一定數(shù)量離散的油粒子來表示,每個油粒子代表一定質(zhì)量的溢油.油粒子的運動包括漂移與擴散.漂移主要由流和風(fēng)的作用引起;擴散則主要由湍流引起[14].本溢油模型計算還包括蒸發(fā)和乳化等變化過程.

    1.2.1 漂移

    漂移包括平流與風(fēng)的拖拽產(chǎn)生的移動.每一個油粒子的平流速度通過所在計算網(wǎng)格節(jié)點的速度(由潮流模型獲得)線性內(nèi)插得到.風(fēng)的拖拽計算中,風(fēng)的影響僅僅作用于溢油,而對流場不再產(chǎn)生影響.風(fēng)的拖曳作用,即將溢油因風(fēng)產(chǎn)生的速度增量表達為風(fēng)速的百分比進行計算,其關(guān)系見式(4).漂移過程通過拉格朗日方法進行模擬,其表達式如式(5)所示.

    (4)

    (5)

    1.2.2 擴散

    在釋放粒子后的初始階段,粒子團相對較小,這時粒子的擴散作用只是由小規(guī)模的湍流作用引起.當(dāng)“粒子云團”充分?jǐn)U散開來,大規(guī)模的渦流將對粒子的擴散產(chǎn)生明顯的影響.故本模型采用隨時間變化紊動擴散系數(shù)D,使得模擬溢油所受湍流過程亦隨時間變化,以更好地模擬大規(guī)模渦流和環(huán)流對后期溢油的擴散作用.紊動擴散系數(shù)D的冪函數(shù)表達式為

    D=atb

    (6)

    式中:t為時間,以粒子釋放時間為初始時刻(t=0);a和b為紊動擴散系數(shù)的校正參數(shù),且0

    湍流擴散是在水動力速度場的基礎(chǔ)上使油粒子產(chǎn)生一個附加位移.這個附加位移的方向是隨機選取的,通過蒙特卡羅方法來模擬.蒙特卡羅方法即利用計算機通過構(gòu)造符合一定規(guī)則的隨機數(shù)以模擬隨機現(xiàn)象進行近似計算的方法.在紊動擴散的過程中,其隨機性的模擬基于一個白噪聲過程,其隨機數(shù)的取值范圍為-1~+1,且均值為0.每個油粒子在一個時間步長Δt內(nèi)的擴散位移為

    (7)

    式中:ΔS為單位時間步長的擴散位移;R為隨機數(shù),取值范圍為[-1,1].

    1.2.3 蒸發(fā)

    PART溢油模型將溢油分為揮發(fā)組分與不揮發(fā)組分,不揮發(fā)組分不發(fā)生蒸發(fā)過程,揮發(fā)組分將產(chǎn)生一階衰減的蒸發(fā)過程.表達如下:

    (8)

    式中:Fvol為揮發(fā)組分分?jǐn)?shù);Fv為已經(jīng)蒸發(fā)掉的組分分?jǐn)?shù);k為蒸發(fā)率.

    1.2.4 乳化

    乳化過程本身是相對快速的,在實驗室中乳化過程從發(fā)生到乳化完成大約需要0.1~3 h[15].但形成油包水乳化物的最長時間一般在10~100 h之間,這主要是因為溢油發(fā)生后往往要經(jīng)過一段時間乳化作用才逐漸發(fā)生.乳化幾乎是一個不可逆轉(zhuǎn)的過程,它使溢油變?yōu)槊芏容^大、具有高粘度的半固體物質(zhì).乳化過程根據(jù)Mackay等人[16]所提出的算法進行計算,其吸水率為

    (9)

    乳化過程不僅起到改變粘度的作用,它還將影響到溢油的蒸發(fā)過程和水下分散過程.根據(jù)Fingas[17]的研究,乳化過程能夠使溢油粘度增大2~3個數(shù)量級,溢油擴散到水體中的過程減慢,而蒸發(fā)過程將近停止.為了使乳化過程和蒸發(fā)過程之間形成關(guān)聯(lián),計算中乳化物的含水率的增大將減小蒸發(fā)率,模型通過減小溢油中蒸發(fā)組分的大小來實現(xiàn),相應(yīng)的計算公式如下:

    (10)

    式中:Fvol為乳化發(fā)生前溢油的蒸發(fā)組分;Few為乳化發(fā)生后溢油的蒸發(fā)組分.當(dāng)乳化物達到最大含水率時,F(xiàn)ew取值為0,蒸發(fā)過程終止.

    1.2.5 附岸

    附岸過程即溢油粘附于岸邊并不再進入水體的過程.當(dāng)溢油進入近岸海域后將發(fā)生附岸過程,并對海岸生態(tài)環(huán)境造成一定的破壞.本模型通過設(shè)定一個粘附概率進行模擬.當(dāng)油粒子進入岸邊或沉至床底時會獲得一個介于0~1之間的一個隨機數(shù),當(dāng)隨機數(shù)小于粘附概率時則油粒子粘附于岸上并不再參與計算.

    2 數(shù)學(xué)模型建立與驗證

    本模型計算范圍為以大連—煙臺為開邊界的整個渤海海域,采用正交曲線網(wǎng)格,網(wǎng)格尺度約為1 000 m.計算網(wǎng)格與邊界如圖1所示.

    圖1 計算區(qū)域與網(wǎng)格

    2.1 水動力模型及驗證

    水動力模型通過環(huán)渤海多個測站的潮位潮流資料進行了驗證.以下僅列出2011年5月20—21日秦皇島與山海關(guān)的潮位驗證如圖2和圖3所示,以及2011年5月20日11點~21日11點V01和V10測站的潮流驗證(圖4).

    圖2 秦皇島潮位驗證(2011-05-20—21)

    圖3 山海關(guān)潮位驗證(2011-05-20—21)

    圖4 V01和V10潮流流速流向驗證 (2011-05-20 11:00—5.21 11:00)

    采用Willmott[18]的統(tǒng)計學(xué)方法對水動力模型計算結(jié)果進行效率評價,計算公式如下:

    (11)

    圖5與圖6分別為2011年6月4日漲急時刻與落急時刻渤海流場圖.渤海灣、遼東灣和萊州灣潮流以往復(fù)流為主,渤海灣漲落潮流以W—E向為主,遼東灣漲落潮流以NE—SW向為主,萊州灣漲落潮流以SW—NE向為主,近岸海域由于岸線影響復(fù)雜多變.渤海中北部以旋轉(zhuǎn)流為主,潮流流向呈順時針旋轉(zhuǎn),漲落潮流向以NW—SE向為主.秦皇島海域的潮流具有明顯的往復(fù)流特征,漲、落潮方向與岸線平行,漲潮方向為西南方向,落潮方向為東北方向,漲、落潮平均流速均較小,并且差別不大.空間分布上,不同海域漲落急時刻并不同步,在三大海灣內(nèi)漲落潮流速變幅較大,而渤海中部區(qū)域?qū)掗?,流速變幅較小,秦皇島海域處于無潮點附近,潮流流速變幅更小.總體上一個潮周期內(nèi)渤海的余流較小.

    圖5 2011年6月4日漲急時刻流場

    圖6 2011年6月4日落急時刻流場

    2.2 溢油模型

    溢油模型將溢油概化為大量的油粒子,每個油粒子代表一定質(zhì)量的溢油, 在拉格朗日方法的基礎(chǔ)上增加了蒙特卡羅法計算溢油擴散隨機性.根據(jù)國家海洋局調(diào)查報告[4],蓬萊19-3(圖1)溢油事故發(fā)生于2011年6月4日,6月21日基本控制溢油,之后較長時間內(nèi)仍有少量溢油溢出.2011年7月中下旬在東戴河岸灘發(fā)現(xiàn)長約4 km,寬約0.5 m呈不均勻帶狀分布的油污帶;在昌黎黃金海岸岸灘高潮線附近發(fā)現(xiàn)長約1.2 km零星分布的油污帶;在河北唐山淺水灣岸灘高潮線附近發(fā)現(xiàn)長約500 m,寬約1~1.5 m呈帶狀分布的油污帶,低潮線附近有長約300 m,寬約1.5~2 m的油污帶.即蓬萊19-3溢油最終在7月中下旬進入東戴河至淺水灣沿岸一帶.根據(jù)溢油最終附岸的時間與位置對模擬結(jié)果進行驗證.

    2.2.1 基本參數(shù)設(shè)置

    模擬對象為2011年蓬萊19-3溢油事故,作為連續(xù)溢油進行考慮,溢油釋放時間為6月4日—7月31日,溢油總量為500 t,總粒子數(shù)為50 000個.溢油密度參數(shù)采用了蓬萊19-3油田中未熟-低熟油密度值(935 kg·m-3)[18].具體參數(shù)設(shè)置如表1所示,其中蒸發(fā)率為每天的蒸發(fā)部分所占比例.

    表1 模型參數(shù)設(shè)置

    2.2.2 風(fēng)拖拽系數(shù)、擴散系數(shù)校正參數(shù)b率定

    溢油運動包括漂移與擴散兩部分:漂移主要由風(fēng)應(yīng)力和潮流引起,其中風(fēng)應(yīng)力起著非常關(guān)鍵的作用;擴散主要為紊動擴散,在初始階段由小規(guī)模的湍流作用引起,當(dāng)粒子“云團”充分?jǐn)U散開來后,大規(guī)模的渦流將對粒子的擴散產(chǎn)生明顯的影響.紊動擴散所引起的附加位移方向是隨機的,可用隨機擴散作用來模擬.在溢油數(shù)值模擬中,風(fēng)拖拽系數(shù)和擴散系數(shù)對溢油的運移與擴散起著關(guān)鍵性的作用.不同的石油具有不同的組分和性質(zhì),同一環(huán)境條件下性質(zhì)不同的溢油在漂移、擴散等變化過程中也會呈現(xiàn)出不同的結(jié)果.對于不同性質(zhì)的溢油,其風(fēng)拖拽系數(shù)和擴散系數(shù)是時空變化的,但是在一定的范圍內(nèi)變化需要通過率定確定.

    由于風(fēng)在溢油的漂移當(dāng)中起著最重要的作用,風(fēng)拖拽系數(shù)的大小對于溢油模擬的準(zhǔn)確性也非常關(guān)鍵.在以往經(jīng)驗當(dāng)中,風(fēng)拖拽系數(shù)一般在1~3%之間進行取值,根據(jù)不同的溢油性質(zhì)進行選擇和校準(zhǔn)[14].為了探索渤海風(fēng)拖拽系數(shù)的合理范圍,本次研究設(shè)置了不同的風(fēng)拖拽系數(shù),根據(jù)溢油的附岸時間與實際情況進行比較來確定.從圖7可以看出,在6種風(fēng)拖拽系數(shù)下的模擬結(jié)果均與實際情況較為接近,當(dāng)Cωd=1.4%~2.3%時,溢油在7月中下旬進入東戴河至淺水灣沿岸一帶,與實際情況較為相符;而當(dāng)Cωd=1.1%、Cωd=2.6%時,溢油略微偏早或偏晚進入近岸海域.在適當(dāng)范圍內(nèi)對風(fēng)拖拽系數(shù)進行取值對于溢油的漂移路徑?jīng)]有太大影響;而由于較大的拖拽系數(shù)下溢油的漂移速度較快,在漂移相同的距離下向兩邊擴散的范圍較小,故溢油的分布形狀較為狹長.根據(jù)模擬結(jié)果與實際情況的對比,對密度大、粘度高的蓬萊19-3溢油進行溢油模擬時,風(fēng)拖拽系數(shù)應(yīng)取1.4%~2.3%.以下蓬萊19-3溢油模擬計算中風(fēng)拖拽系數(shù)取值為Cωd=2%.

    b Cωd=1.4%

    c Cωd=1.7%

    d Cωd=2.0%

    e Cωd=2.3%

    f Cωd=2.6%

    紊動擴散系數(shù)式(6)中含兩個可變參數(shù)a和b.根據(jù)式(6)可知:參數(shù)a的影響較小,可取默認(rèn)值1.0;對于時間跨度較長的溢油模擬,參數(shù)b起主要作用.本文主要針對參數(shù)b的影響進行探討,通過設(shè)置不同的參數(shù)b進行模擬與對比,根據(jù)溢油的附岸范圍和實際情況進行比較以確定合適的參數(shù)范圍.模擬結(jié)果如圖8所示,結(jié)果表明:參數(shù)b的改變對溢油的漂移速度和路徑影響不大,對溢油的擴散范圍影響較大.當(dāng)b=0.38~0.42時,溢油于東戴河至淺水灣沿岸一帶附岸,與實際情況相符;當(dāng)b=0.3,0.35時,淺水灣沿岸未出現(xiàn)溢油附岸,溢油擴散范圍偏小;當(dāng)b=0.45時,溢油附岸范圍兩端分別略微超出了東戴河與淺水灣,擴散范圍偏大.因此,在利用本模型對密度大、粘度高的蓬19-3溢油進行模擬時,參數(shù)b合理的取值范圍為0.38~0.42.以下蓬萊19-3溢油模擬取b=0.4.經(jīng)計算,時間步長設(shè)為1 h,則1 h后紊動擴散系數(shù)D=26.5,則最大擴散位移ΔS=756 m;24 h后D=94.3,則ΔS=1 427 m;10 d后D=236.9,ΔS=2 262 m;50 d后D=451.0,ΔS=3 122 m.模擬期間紊動擴散系數(shù)取值大概在26.5~500范圍內(nèi).由于擴散位移的方向是隨機選取的,故擴散位移主要使擴散范圍逐漸增大,對溢油質(zhì)心位移的影響較小.可見,溢油后期的紊動擴散系數(shù)增大,油粒子的擴散位移也明顯增大.這樣大渦促使溢油后期迅速擴散的作用便得以體現(xiàn)出來.

    a b=0.3

    b b=0.35

    c b=0.38

    d b=0.4

    e b=0.42

    f b=0.45

    3 模擬結(jié)果與分析

    模擬采用風(fēng)場由美國NOAA提供的9個渤海及臨近區(qū)域站點的風(fēng)場數(shù)據(jù),其中渤海中部站點(120°E,39.047°N)的風(fēng)速風(fēng)向如圖9所示,平均風(fēng)速為4.10 m·s-1,最大風(fēng)速為12.04 m·s-1;常風(fēng)向為南向,南向風(fēng)占總數(shù)的25.4%,平均風(fēng)速為5.02 m·s-1;南偏東、南、南偏西風(fēng)向占總數(shù)的57.3%,平均風(fēng)速為4.61 m·s-1.

    在前面參數(shù)率定(即風(fēng)拖拽系數(shù)Cωd=2%、擴散系數(shù)校正參數(shù)b=0.4)的基礎(chǔ)上,對蓬萊19-3溢油事故進行模擬,結(jié)果如圖10.溢油發(fā)生后擴散范圍逐漸增大,前端擴散面積較大,平面分布密度較小,約為0.000 02 kg·m-2.溢油平臺周圍及其北偏西方向溢油分布密度較大,大約為0.000 05~0.000 5 kg·m-2.結(jié)合圖9與圖10進行分析可知:6月22日以前以南風(fēng)為主,溢油運移方向與風(fēng)向大致相同為正北方向;6月22日—7月4日,風(fēng)向以東南偏東向和南向為主,溢油逐漸轉(zhuǎn)為向西北方向運移,并靠近西北沿岸一帶;7月4日—7月16日,風(fēng)向以南向為主,溢油往北運移進入東戴河至黃金海岸沿岸一帶;7月16日—7月22日,風(fēng)向以南偏東向為主,溢油往西北方向運移,附岸范圍向南擴展至淺水灣沿岸;隨后風(fēng)向以南向為主,溢油擴散范圍略微往北偏移.總體上,溢油油污于7月中下旬進入淺水灣至東戴河沿岸一帶,其中東戴河與黃金海岸沿岸一帶受到的溢油污染最為嚴(yán)重,淺水灣受到的污染較輕,模擬結(jié)果與實際情況相符.

    渤海海域夏季以南風(fēng)和東南風(fēng)為主,風(fēng)速較大且風(fēng)向較為穩(wěn)定,溢油的漂移方向與主風(fēng)向大致相同,表明風(fēng)對漂浮于水面上的溢油的運移起著很明顯的作用.然而,渤海海域的潮流流速不大,尤其是秦皇島海域處于無潮點區(qū)域,潮流流速更小.渤海潮流以往復(fù)流和旋轉(zhuǎn)流為主,漲潮潮流與落潮潮流流向大致相反,且流速差異不明顯.一個潮周期內(nèi)溢油在漲落潮流的作用下會來回擺動,余流小、凈位移小.因此,在渤海海域內(nèi)潮流對溢油的漂移不及風(fēng)的作用.

    圖9 風(fēng)速風(fēng)向過程 (2011-06-04—7.31)

    a 6月14日

    b 6月24日

    c 7月4日

    d 7月14日

    e 7月24日

    f 7月31日

    溢油模型計算過程中,溢油主要分為漂浮組分、蒸發(fā)組分和附岸組分.圖11為整個溢油模擬過程中各組分的質(zhì)量變化圖.其中溢油總量設(shè)定為線性增長,16 d(6月19日)基本控制溢油,溢油總量達290 t;后期依然持續(xù)有少量溢油溢出.在溢油前期,溢油漂浮組分和蒸發(fā)組分占溢油總量的絕大部分.漂浮組分在初始階段隨著溢油的發(fā)生呈線性增長,36 d(7月10日)達到最大值204 t,約占溢油總量的50.2%.之后隨著漂浮于水面的溢油逐漸發(fā)生蒸發(fā)與附岸使其逐漸下降,第56 d(7月31日),漂浮組分下降至約155.0 t,占溢油總量30.6%.溢油的蒸發(fā)主要發(fā)生于溢油初期,由于前25 d溢油量較大,故蒸發(fā)量增長較快,于第25 d便達到約160 t,約占溢油總量48.2%;之后隨著溢油揮發(fā)組分減少,且乳化過程使溢油含水率增大,使得蒸發(fā)在溢油的后期階段速度減緩,第57 d(7月31日)蒸發(fā)量達到最大值約為223 t,占溢油總量的44.0%.其中,第22 d(6月26日)出現(xiàn)了12.04 m·s-1的最大風(fēng)速,溢油的蒸發(fā)出現(xiàn)了短暫的陡增,漂浮組分也隨之出現(xiàn)驟減,可見風(fēng)對促進溢油蒸發(fā)也有一定作用.至第35 d(7月9日)之后,溢油進入沿岸海域并逐漸發(fā)生附岸,于7月31日達126 t左右,占溢油總量的24.9%.

    圖11 溢油組分質(zhì)量變化圖

    4 結(jié)論

    通過Delft3D中的FLOW模塊建立了渤海水動力模型,在水動力模型的基礎(chǔ)上利用PART模塊耦合建立了渤海海上溢油模型.率定了該溢油模型中風(fēng)拖拽系數(shù)、擴散系數(shù)校正參數(shù)b的合理范圍,對2011年蓬萊19-3溢油事故進行了模擬與驗證,并對溢油運移與歸宿進行了相關(guān)探討.得出以下結(jié)論:

    (1) 渤海海域潮流以往復(fù)流和旋轉(zhuǎn)流為主,其中渤海中北部以旋轉(zhuǎn)流為主,呈順時針旋轉(zhuǎn),漲落潮流向以NW-SE向為主.渤海海域潮流流速較小,且漲落潮流速相差小,變幅也小.其中,秦皇島海域處于無潮點附近,潮流流速更小.

    (2) 對于密度大粘度高的蓬萊19-3溢油風(fēng)拖拽系數(shù)合理范圍為1.4~2.3%;擴散系數(shù)公式中校正參數(shù)b對溢油的擴散范圍影響較大,其合理取值范圍為0.38~0.42.

    (3) 渤海夏季風(fēng)速較大、風(fēng)向較為穩(wěn)定的南風(fēng)和東南風(fēng)對漂浮于水面上的溢油作用較大,且溢油漂移方向與主風(fēng)向大致相同.而渤海潮流以往復(fù)流和旋轉(zhuǎn)流為主,溢油在漲落潮流的作用下會出現(xiàn)來回擺動,凈位移小,潮流對溢油的漂移不及風(fēng)的作用.

    (4) 溢油初期的蒸發(fā)速度較快;由于揮發(fā)組分的減少以及乳化作用使溢油含水率增大,溢油后期的蒸發(fā)速度明顯減小.進入沿岸海域后溢油容易發(fā)生附岸,附岸和蒸發(fā)使得溢油漂浮組分逐漸減少.最終,漂浮組分、蒸發(fā)組分和附岸組分分別占溢油總量的30.6%、44.0%和24.9%,蒸發(fā)和附岸是溢油的主要歸宿.

    [1] 張麗萍,王輝.不同海面狀況海洋石油污染處理方法優(yōu)化配置[J].海洋技術(shù),2006,25(3):1.

    ZHANG Liping, WANG Hui. Optimization design of countermeasures to cleanup the marine oil pollution in different sea surface conditions.[J].Ocean Technology,2006,25(3):1.

    [2] 徐艷東.海上溢油風(fēng)化過程及其預(yù)測模型研究[D].青島:中國海洋大學(xué),2006.

    XU Yandong. A study on theweathering and predicting model of spilled oils at sea[D]. Qingdao: Ocean University of China,2006.

    [3] DENG Zenggan, YU Ting, JIANG Xiaoyi,etal. Bohai Sea oil spill model: a numerical case study[J]. Marine Geophysical Research.2013,34(2):115.

    [4] 國家海洋局.蓬萊19-3油田溢油事故聯(lián)合調(diào)查組關(guān)于事故調(diào)查處理報告[R].北京:國家海洋局,2011[2015-02-03]. http://www.soa.gov.cn/xw/hyyw_90/201211/t20121109_884.html.

    State Oceanic Administration People’s Republic of China. The joint investigation group report on the investigation and handling of the Penglai 19-3 oil spill accident[R].Beijing: State Oceanic Administration People’s Republic of China, 2011[2015-02-03]. http://www.soa.gov.cn/xw/hyyw_90/201211/t20121109_884.html.

    [5] KUANG Cuiping, PAN Yi, HU Ying,etal. Primary analysis of the 2011ConocoPhillips oil spill in Bohai Bay, China [C]// Advances in Biomedical Engineering——2012 Asia Pacific Conference on Environmental Science and Technology (APEST 2012). Kuala Lumpur: Information Engineering Research Institute,2012:1-5.

    [6] DENG Zenggan, YU Ting, SHI Suixiang,etal. Numerical study of the oil spill trajectory in Bohai Sea, China[J]. Marine Geodesy,2013,36(4): 351.

    [7] XU Qing, LI Xiaofeng, WEI Yongliang,etal. Satellite observations and modeling of oil spill trajectories in the Bohai Sea[J]. Marine Pollution Bulletin,2013,71(1/2):107.

    [8] WANG Shoudong, SHENYongming, ZHENG Y.H. Two-dimensional numerical simulation for transportand fate of oil spills in seas[J]. Ocean Engineering,2005,32(13):1556.

    [9] WANG Shoudong, SHEN Yongming, GUO Yakun,etal. Three-dimensional numerical simulation for transport of oil spills in seas[J]. Ocean Engineering,2008,35(5): 503.

    [10] WANG Jinhua, SHEN Yongming. Development of an integrated model system to simulate transport and fate of oil spills in seas[J].Sciences China-Technological Sciences,2010,53(9): 2423.

    [11] BI Haipu, SI Hu. Numerical simulation of oil spill for the Three Gorges Reservoir in China[J]. Water and Environment Journal,2014,28(2):183.

    [12] LI Yan, ZHU Jiang, WANG Hui. The impact of different vertical diffusion schemes in a three-dimensional oil spill model in the BohaiSea[J]. Advances in Atmospheric Sciences,2013, 30(6): 1569.

    [13] WL|DelftHydraulics. User manual Delft3D-FLOW[M].Delft: WL|Delft Hydraulics, 2013.

    [14] WL|DelftHydraulics. User manual D-WaqPART[M].Delft: WL|Delft Hydraulics, 2013.

    [15] Fingas M, Filedhouse B, Mullin J. Water-in-oil emulsions results of formation studies and applicability to oil spill modeling[J].Spill Science & Technology,1999,5(1): 81.

    [16] Mackay D, Paterson S, Trudel K. A mathematical model of oil spill behavior on water with natural and chemical dispersion. [R]. Ottawa: Environment Canada, 1977.

    [17] Fingas M F. Chemistry of oil and modelling of spills[J].Journal of Advances in Marine Technology Conference,1994, 11: 41.

    [18] Willmott C J. On the validation of models[J].Physical Geography,1981, 2(2):184.

    [19] 王廣源,周心懷,王昕,等.蓬萊19-3/25-6油田未熟-低熟油特征與成因[J].石油實驗地質(zhì), 2014,36(2):231.

    WANG Guangyuan, ZHOU Xinhuai, WANG Xin,etal. Characteristics and genesis of immature and low-matire oils in Penglai 19-3/25-6 oilfields[J]. Petroleum Geology & Experiment, 2014, 36(2):231.

    Numerical Simulation and Analysis of Transport and Fate of Penglai 19-3 Oil Spill

    KUANG Cuiping1, XIE Hualang1, SU Ping2, CHEN KUO1

    (1. College of Civil Engineering, Tongji University, Shanghai 200092, China; 2. Guangxi Communications Planning Surveying and Designing Institute, Nanning 530029, China)

    The numerical model of oil spill transport and fate in the Bohai Sea is set up by coupling with the hydrodynamic model. The change processes of spilled oil included in the model are drift, dispersion, evaporation, emulsification and so on. The drift process is simulated by the Lagrange method, while the turbulent dispersion process is considered by using the ‘Monte Carlo method’, turbulent dispersion of spilled oil changes over time by using the time-varying turbulent dispersion coefficient. The reasonable ranges of wind drag coefficient and the dispersion correction parameterbfor Penglai 19-3 oil spill model obtained by calibration are 1.4~2.3% and 0.38~0.42, respectively. Then, the Penglai 19-3 oil spill accident in the Bohai Sea is used to verify the oil spill model by comparing the time and location of spilled oil arrival at shore, and the simulation results are consistent with the actual situation. The numerical results show that the wind, because of its large velocity and small change in direction, plays a major role on the drift of spilled oil floating on the water surface, and the drift direction of spilled oil is roughly the same with the main direction of wind. The effect of the tidal current, which is dominated by reversing current and rotary current with a very weak residual current, to the drift of spilled oil in Bohai Sea is less than the effect of the wind. The evaporation rate of spilled oil is quite large in the initial stage after oil spill, but significantly reduced in the later stage as a result of the reduction of the volatile fraction and the increase of the water content caused by emulsification. The reduction of floating oil is mainly caused by evaporation and sticking that are also the main fates of spilled oil.

    Bohai Sea; oil spill; transport and dispersion; numerical simulation

    2015-02-06

    國家海洋公益項目(201305003)

    匡翠萍(1966—),女,教授,博士生導(dǎo)師,工學(xué)博士,主要研究方向為河口海岸工程和環(huán)境.E-mail:cpkuang@#edu.cn

    X55; P76

    A

    猜你喜歡
    蓬萊溢油擴散系數(shù)
    蓬萊迎曦
    寶藏(2020年10期)2020-11-19 01:47:32
    近岸溢油漂移擴散預(yù)測方法研究——以膠州灣溢油事件為例
    海洋通報(2020年2期)2020-09-04 09:22:48
    基于GF-1衛(wèi)星的海上溢油定量監(jiān)測——以青島溢油事故為例
    海洋通報(2020年2期)2020-09-04 09:22:46
    煙臺 身在蓬萊就是仙
    蓬萊凝翠
    寶藏(2018年1期)2018-04-18 07:39:26
    岱山五云縹緲隔蓬萊
    中國三峽(2017年1期)2017-06-09 11:09:41
    基于Sauer-Freise 方法的Co- Mn 體系fcc 相互擴散系數(shù)的研究
    上海金屬(2015年5期)2015-11-29 01:13:59
    FCC Ni-Cu 及Ni-Mn 合金互擴散系數(shù)測定
    上海金屬(2015年6期)2015-11-29 01:09:09
    非時齊擴散模型中擴散系數(shù)的局部估計
    對白茆沙水域溢油事故后修復(fù)治理的思考
    中國水利(2015年4期)2015-02-28 15:12:23
    亚洲乱码一区二区免费版| 日韩伦理黄色片| 五月伊人婷婷丁香| 嫩草影院精品99| 丝袜美腿在线中文| 国产高清有码在线观看视频| 全区人妻精品视频| 日日摸夜夜添夜夜添av毛片| 国产成人精品一,二区| 免费看不卡的av| 久久久久久久久中文| 中文字幕人妻熟人妻熟丝袜美| 亚洲av成人精品一区久久| 日韩一区二区视频免费看| 久久鲁丝午夜福利片| 99久久九九国产精品国产免费| 国产精品一区二区三区四区免费观看| 人妻一区二区av| 久久久久免费精品人妻一区二区| 精品国内亚洲2022精品成人| 精品亚洲乱码少妇综合久久| 亚洲图色成人| 别揉我奶头 嗯啊视频| 久久久久久久亚洲中文字幕| 精品久久久噜噜| av专区在线播放| 精品亚洲乱码少妇综合久久| 街头女战士在线观看网站| 女人久久www免费人成看片| 国产探花极品一区二区| 不卡视频在线观看欧美| 看十八女毛片水多多多| av在线蜜桃| 欧美日韩精品成人综合77777| 国产精品三级大全| 五月伊人婷婷丁香| 久久久色成人| 久久久精品94久久精品| 久久精品综合一区二区三区| 毛片一级片免费看久久久久| 99久国产av精品国产电影| 日本欧美国产在线视频| 熟女人妻精品中文字幕| 国产精品美女特级片免费视频播放器| 成人无遮挡网站| 大片免费播放器 马上看| 在线观看美女被高潮喷水网站| 免费观看性生交大片5| 男女视频在线观看网站免费| 午夜精品在线福利| 日韩欧美国产在线观看| 麻豆久久精品国产亚洲av| av又黄又爽大尺度在线免费看| 18禁裸乳无遮挡免费网站照片| 国产乱人偷精品视频| 国产男女超爽视频在线观看| 国产精品国产三级国产av玫瑰| 国模一区二区三区四区视频| 别揉我奶头 嗯啊视频| .国产精品久久| 亚洲av成人精品一二三区| av卡一久久| 色视频www国产| 国产成人午夜福利电影在线观看| 街头女战士在线观看网站| 视频中文字幕在线观看| 国产久久久一区二区三区| 看十八女毛片水多多多| 亚洲精品色激情综合| 亚洲欧美清纯卡通| 日韩电影二区| 国产亚洲5aaaaa淫片| 亚洲av福利一区| 国产精品熟女久久久久浪| 亚洲国产精品国产精品| 国产淫语在线视频| 国产伦在线观看视频一区| 亚洲在久久综合| 国内精品美女久久久久久| 男女啪啪激烈高潮av片| 日韩制服骚丝袜av| 国产 亚洲一区二区三区 | 丝瓜视频免费看黄片| 青春草国产在线视频| 亚洲婷婷狠狠爱综合网| 老司机影院成人| 午夜激情福利司机影院| 18禁动态无遮挡网站| 免费黄网站久久成人精品| 免费看不卡的av| av.在线天堂| av免费观看日本| 久久久久久久久大av| 亚洲国产精品专区欧美| 亚洲18禁久久av| 精品国产一区二区三区久久久樱花 | 一个人观看的视频www高清免费观看| 国产精品一区www在线观看| 干丝袜人妻中文字幕| 亚洲在线观看片| 日韩,欧美,国产一区二区三区| 国产午夜福利久久久久久| 国产 一区精品| 丰满人妻一区二区三区视频av| 亚洲乱码一区二区免费版| 大陆偷拍与自拍| 中文字幕av成人在线电影| 少妇高潮的动态图| 熟女电影av网| 日本午夜av视频| 午夜激情久久久久久久| 国产亚洲最大av| 欧美性猛交╳xxx乱大交人| 午夜福利在线观看吧| 伦精品一区二区三区| 99久国产av精品| 青春草视频在线免费观看| 80岁老熟妇乱子伦牲交| 人妻制服诱惑在线中文字幕| 麻豆久久精品国产亚洲av| 亚洲熟女精品中文字幕| www.色视频.com| 日日干狠狠操夜夜爽| 国产精品一二三区在线看| 内地一区二区视频在线| av线在线观看网站| 国产大屁股一区二区在线视频| 亚洲图色成人| 亚洲av在线观看美女高潮| 亚洲久久久久久中文字幕| 在线观看人妻少妇| 特级一级黄色大片| 看十八女毛片水多多多| 日韩精品青青久久久久久| 美女黄网站色视频| 亚洲av成人av| 久久久国产一区二区| 69av精品久久久久久| 91aial.com中文字幕在线观看| 两个人的视频大全免费| 亚洲av中文av极速乱| 亚洲美女搞黄在线观看| 国产av不卡久久| 国产成人精品一,二区| 亚洲国产精品国产精品| 日韩 亚洲 欧美在线| 久久精品国产亚洲av天美| av专区在线播放| 欧美日韩一区二区视频在线观看视频在线 | 国产精品女同一区二区软件| 免费观看的影片在线观看| 高清欧美精品videossex| 精品久久久久久久末码| 一区二区三区免费毛片| 久热久热在线精品观看| 一个人看视频在线观看www免费| 国产精品无大码| 少妇裸体淫交视频免费看高清| 又爽又黄无遮挡网站| 日韩,欧美,国产一区二区三区| av国产免费在线观看| 亚洲精品久久午夜乱码| 久久精品久久久久久久性| 亚洲无线观看免费| 亚洲国产精品sss在线观看| 男女国产视频网站| 色视频www国产| 欧美+日韩+精品| 国产一区二区在线观看日韩| 欧美日韩视频高清一区二区三区二| 好男人在线观看高清免费视频| 国产亚洲av片在线观看秒播厂 | 大陆偷拍与自拍| 99久久中文字幕三级久久日本| 久久久久九九精品影院| 老师上课跳d突然被开到最大视频| 大香蕉久久网| 国产精品无大码| 天天躁日日操中文字幕| 超碰97精品在线观看| 国产精品.久久久| av播播在线观看一区| 日韩成人av中文字幕在线观看| 成年av动漫网址| 大陆偷拍与自拍| 又粗又硬又长又爽又黄的视频| 日韩亚洲欧美综合| 成年女人在线观看亚洲视频 | 亚洲av国产av综合av卡| 最近视频中文字幕2019在线8| 非洲黑人性xxxx精品又粗又长| 久久精品熟女亚洲av麻豆精品 | 国产精品久久久久久精品电影| 欧美激情在线99| 国产精品伦人一区二区| 少妇的逼水好多| 寂寞人妻少妇视频99o| 搡女人真爽免费视频火全软件| 一区二区三区四区激情视频| 伊人久久国产一区二区| 亚洲欧美一区二区三区黑人 | 国产成人福利小说| 在现免费观看毛片| 在线观看人妻少妇| 国产精品不卡视频一区二区| 神马国产精品三级电影在线观看| 国产精品爽爽va在线观看网站| 免费看av在线观看网站| 日本熟妇午夜| 亚洲精品日韩在线中文字幕| 亚洲18禁久久av| 六月丁香七月| 国产精品国产三级国产专区5o| 亚洲在线自拍视频| 九草在线视频观看| 男女啪啪激烈高潮av片| 美女黄网站色视频| 卡戴珊不雅视频在线播放| 禁无遮挡网站| 亚洲精品乱码久久久v下载方式| 色5月婷婷丁香| 在线观看一区二区三区| 国产av国产精品国产| 18禁动态无遮挡网站| 国产伦精品一区二区三区四那| 精品熟女少妇av免费看| 亚洲国产精品专区欧美| 国产人妻一区二区三区在| 久久久精品免费免费高清| 亚洲精品视频女| 观看美女的网站| 亚洲国产成人一精品久久久| av天堂中文字幕网| 精品久久久久久电影网| 91久久精品国产一区二区成人| 欧美日韩亚洲高清精品| 国产毛片a区久久久久| 哪个播放器可以免费观看大片| 校园人妻丝袜中文字幕| 午夜福利成人在线免费观看| 国产一区二区三区综合在线观看 | 久久久色成人| 国产日韩欧美在线精品| 亚洲精品aⅴ在线观看| 国产精品蜜桃在线观看| 久久久久免费精品人妻一区二区| 久久精品人妻少妇| 男女边吃奶边做爰视频| 国产美女午夜福利| 亚洲精华国产精华液的使用体验| 免费高清在线观看视频在线观看| 一个人免费在线观看电影| 丝瓜视频免费看黄片| 亚洲成色77777| 久久久久久伊人网av| 欧美成人一区二区免费高清观看| 日韩欧美国产在线观看| 在线观看av片永久免费下载| 亚洲精品国产成人久久av| 非洲黑人性xxxx精品又粗又长| 网址你懂的国产日韩在线| xxx大片免费视频| 水蜜桃什么品种好| 熟妇人妻久久中文字幕3abv| 久久这里只有精品中国| 日本av手机在线免费观看| 婷婷色av中文字幕| 亚洲经典国产精华液单| 97人妻精品一区二区三区麻豆| 亚洲成人一二三区av| 伦精品一区二区三区| 日韩中字成人| 男女视频在线观看网站免费| 别揉我奶头 嗯啊视频| 高清午夜精品一区二区三区| 国内揄拍国产精品人妻在线| av国产久精品久网站免费入址| 别揉我奶头 嗯啊视频| 日本-黄色视频高清免费观看| 亚洲精品色激情综合| 男插女下体视频免费在线播放| 啦啦啦韩国在线观看视频| 1000部很黄的大片| 乱人视频在线观看| 简卡轻食公司| 国产亚洲精品av在线| 免费播放大片免费观看视频在线观看| 美女高潮的动态| 我的老师免费观看完整版| 性插视频无遮挡在线免费观看| 天堂网av新在线| 26uuu在线亚洲综合色| 久久精品国产亚洲网站| 免费观看在线日韩| 欧美日韩视频高清一区二区三区二| 亚洲天堂国产精品一区在线| 国产老妇女一区| 亚洲无线观看免费| 亚洲精品成人久久久久久| 亚洲aⅴ乱码一区二区在线播放| 午夜福利成人在线免费观看| 波野结衣二区三区在线| 美女国产视频在线观看| 精品久久久久久久末码| 蜜桃久久精品国产亚洲av| 国产精品综合久久久久久久免费| 亚洲av成人精品一二三区| 亚洲综合精品二区| 大又大粗又爽又黄少妇毛片口| av在线蜜桃| 综合色av麻豆| 免费看av在线观看网站| av.在线天堂| 青春草亚洲视频在线观看| 一级黄片播放器| 亚洲国产日韩欧美精品在线观看| 亚洲欧洲国产日韩| 99久国产av精品国产电影| 人人妻人人看人人澡| 26uuu在线亚洲综合色| 97热精品久久久久久| 精品久久久噜噜| 亚洲真实伦在线观看| 好男人视频免费观看在线| av在线老鸭窝| 蜜臀久久99精品久久宅男| 久久6这里有精品| 一级爰片在线观看| 国产精品无大码| 亚洲av免费高清在线观看| 97精品久久久久久久久久精品| 超碰97精品在线观看| 欧美最新免费一区二区三区| 少妇的逼好多水| 国产精品伦人一区二区| 大话2 男鬼变身卡| 久久久精品94久久精品| 免费黄频网站在线观看国产| 国产女主播在线喷水免费视频网站 | 亚洲精品影视一区二区三区av| 日韩欧美一区视频在线观看 | 国产老妇伦熟女老妇高清| 国产一区二区三区综合在线观看 | 中文乱码字字幕精品一区二区三区 | 日韩中字成人| 久久久久精品性色| 久久久久久久久久人人人人人人| 久久久久久国产a免费观看| 成人鲁丝片一二三区免费| 国产精品一区二区性色av| 狂野欧美激情性xxxx在线观看| 一本一本综合久久| 国产精品爽爽va在线观看网站| 91精品国产九色| 三级男女做爰猛烈吃奶摸视频| 色综合站精品国产| 免费黄频网站在线观看国产| 久久精品夜夜夜夜夜久久蜜豆| 久久人人爽人人片av| 天天一区二区日本电影三级| 2018国产大陆天天弄谢| 国产精品久久久久久久久免| 欧美xxxx性猛交bbbb| 久久久久久久大尺度免费视频| .国产精品久久| 深爱激情五月婷婷| 亚洲精品乱码久久久v下载方式| 精品国产露脸久久av麻豆 | 久久久国产一区二区| 三级毛片av免费| 少妇高潮的动态图| 亚洲国产欧美人成| 精品久久久久久成人av| 国产午夜福利久久久久久| 一级毛片电影观看| 特大巨黑吊av在线直播| 三级国产精品欧美在线观看| 免费观看无遮挡的男女| 国产激情偷乱视频一区二区| av国产免费在线观看| 中文乱码字字幕精品一区二区三区 | 自拍偷自拍亚洲精品老妇| 午夜福利在线观看吧| 日韩欧美一区视频在线观看 | 青春草亚洲视频在线观看| 国产色爽女视频免费观看| 五月伊人婷婷丁香| 中文乱码字字幕精品一区二区三区 | 能在线免费看毛片的网站| 国产91av在线免费观看| 久久久亚洲精品成人影院| 成人无遮挡网站| 日日啪夜夜爽| 2021天堂中文幕一二区在线观| 欧美成人a在线观看| 亚洲无线观看免费| 亚洲国产av新网站| 男女视频在线观看网站免费| 美女被艹到高潮喷水动态| 在线天堂最新版资源| 国产欧美另类精品又又久久亚洲欧美| 久久久久久久国产电影| 汤姆久久久久久久影院中文字幕 | 成人二区视频| 精品人妻偷拍中文字幕| 大片免费播放器 马上看| 日韩在线高清观看一区二区三区| 又大又黄又爽视频免费| 久久人人爽人人片av| 成年女人看的毛片在线观看| 免费大片18禁| 秋霞伦理黄片| 国产成人免费观看mmmm| 欧美xxxx性猛交bbbb| 熟女电影av网| 青青草视频在线视频观看| 成年av动漫网址| 免费看a级黄色片| 国产精品熟女久久久久浪| 久久久久久久久久久免费av| 高清av免费在线| 别揉我奶头 嗯啊视频| 亚洲成人中文字幕在线播放| 青春草国产在线视频| 黄色日韩在线| 国精品久久久久久国模美| 国产伦理片在线播放av一区| 床上黄色一级片| 国产91av在线免费观看| 国产午夜精品久久久久久一区二区三区| 春色校园在线视频观看| 欧美不卡视频在线免费观看| 欧美性猛交╳xxx乱大交人| 亚洲在线观看片| 亚洲伊人久久精品综合| 又爽又黄无遮挡网站| av在线蜜桃| 成人无遮挡网站| 日本熟妇午夜| 欧美精品国产亚洲| 99久久人妻综合| 欧美bdsm另类| 秋霞在线观看毛片| 两个人的视频大全免费| 亚洲综合精品二区| 成人亚洲精品一区在线观看 | 亚洲丝袜综合中文字幕| 在线免费观看的www视频| 极品少妇高潮喷水抽搐| 免费播放大片免费观看视频在线观看| 免费看av在线观看网站| 亚洲av国产av综合av卡| 国产淫片久久久久久久久| 久久亚洲国产成人精品v| 久久久久久久亚洲中文字幕| 日本wwww免费看| 国产高清三级在线| 久久精品综合一区二区三区| 中文资源天堂在线| 在线观看免费高清a一片| 99视频精品全部免费 在线| 国产精品久久久久久精品电影小说 | 亚洲美女搞黄在线观看| 天堂√8在线中文| 午夜免费男女啪啪视频观看| 日日摸夜夜添夜夜爱| 亚洲成色77777| 亚洲精品影视一区二区三区av| 十八禁网站网址无遮挡 | 三级经典国产精品| 99久国产av精品国产电影| 男人舔奶头视频| 蜜桃亚洲精品一区二区三区| 免费电影在线观看免费观看| 波多野结衣巨乳人妻| 日本免费a在线| 亚洲色图av天堂| 午夜福利视频1000在线观看| 国产av在哪里看| 欧美高清性xxxxhd video| 日日撸夜夜添| kizo精华| 好男人视频免费观看在线| 边亲边吃奶的免费视频| 欧美成人精品欧美一级黄| 七月丁香在线播放| 18+在线观看网站| 亚洲av一区综合| 国产精品.久久久| 久久久色成人| 亚洲熟女精品中文字幕| 又黄又爽又刺激的免费视频.| 97在线视频观看| 久久精品国产亚洲av涩爱| 午夜福利在线观看吧| 纵有疾风起免费观看全集完整版 | 亚洲欧美一区二区三区黑人 | av免费在线看不卡| 免费看不卡的av| 久久6这里有精品| 亚洲欧美精品自产自拍| 99热6这里只有精品| 久久国产乱子免费精品| 一级爰片在线观看| 视频中文字幕在线观看| 全区人妻精品视频| 99热全是精品| 99热网站在线观看| 联通29元200g的流量卡| 舔av片在线| 日韩国内少妇激情av| 一二三四中文在线观看免费高清| 国产女主播在线喷水免费视频网站 | 国产精品1区2区在线观看.| 久久久久网色| 日韩av在线免费看完整版不卡| 99热6这里只有精品| 免费看a级黄色片| 99久国产av精品| 国产精品无大码| 观看免费一级毛片| 久久人人爽人人爽人人片va| 亚洲性久久影院| 亚洲欧美清纯卡通| 街头女战士在线观看网站| 99久久人妻综合| a级毛片免费高清观看在线播放| 白带黄色成豆腐渣| 狂野欧美白嫩少妇大欣赏| 国内揄拍国产精品人妻在线| 青春草视频在线免费观看| 日日干狠狠操夜夜爽| 国产免费一级a男人的天堂| 99久国产av精品| 午夜激情久久久久久久| www.av在线官网国产| 亚洲国产精品成人久久小说| 男女下面进入的视频免费午夜| 亚洲丝袜综合中文字幕| 久久韩国三级中文字幕| 亚洲精品乱码久久久久久按摩| 狠狠精品人妻久久久久久综合| 中文精品一卡2卡3卡4更新| 国产一区二区三区综合在线观看 | 18禁裸乳无遮挡免费网站照片| 日日摸夜夜添夜夜爱| 1000部很黄的大片| 狠狠精品人妻久久久久久综合| 国产精品三级大全| 国产精品一区二区在线观看99 | 寂寞人妻少妇视频99o| 国产一级毛片在线| 一级毛片 在线播放| 国产伦在线观看视频一区| 国产亚洲av嫩草精品影院| 亚洲精品,欧美精品| 国产精品麻豆人妻色哟哟久久 | 最近中文字幕2019免费版| 免费在线观看成人毛片| 在线免费观看的www视频| 亚洲在久久综合| 国产淫片久久久久久久久| 少妇丰满av| 国产激情偷乱视频一区二区| 好男人视频免费观看在线| av在线观看视频网站免费| 亚洲天堂国产精品一区在线| 综合色丁香网| 欧美xxⅹ黑人| 成人亚洲精品av一区二区| 麻豆国产97在线/欧美| 波野结衣二区三区在线| 久久99热6这里只有精品| 亚洲精品色激情综合| 人妻一区二区av| 丝袜美腿在线中文| 亚洲电影在线观看av| 国产麻豆成人av免费视频| 国产伦精品一区二区三区视频9| 国产淫语在线视频| 亚洲精品影视一区二区三区av| 国产精品国产三级专区第一集| 日韩欧美国产在线观看| 校园人妻丝袜中文字幕| 午夜福利高清视频| 又爽又黄a免费视频| 亚洲av一区综合| 国产女主播在线喷水免费视频网站 | 婷婷色麻豆天堂久久| av在线老鸭窝| 69av精品久久久久久| freevideosex欧美| 最近最新中文字幕大全电影3| 久久精品久久精品一区二区三区| 日韩人妻高清精品专区| 中文在线观看免费www的网站| 亚洲av一区综合| 成人毛片a级毛片在线播放| 亚洲美女搞黄在线观看| 欧美性感艳星| 日韩欧美 国产精品| 国产乱人偷精品视频| 亚洲一级一片aⅴ在线观看| 十八禁网站网址无遮挡 | 日韩精品有码人妻一区| 80岁老熟妇乱子伦牲交| 男插女下体视频免费在线播放| 国产一区二区在线观看日韩| 国产精品人妻久久久久久| 男人和女人高潮做爰伦理| 少妇人妻一区二区三区视频| 亚洲av中文字字幕乱码综合| 午夜福利视频1000在线观看| 又大又黄又爽视频免费| 国产伦理片在线播放av一区|