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

    水/氣多介質(zhì)問題的界面處理方法*

    2015-06-07 11:38:12王春武王東紅
    爆炸與沖擊 2015年3期
    關(guān)鍵詞:激波邊界條件流場

    徐 爽,趙 寧,王春武,王東紅

    (1.南京 航空航天大 學(xué)航空宇航 學(xué)院,江蘇 南 京 210016; 2.南京航空航天大學(xué)理學(xué)院,江蘇 南京210016)

    水/氣多介質(zhì)問題的界面處理方法*

    徐 爽1,趙 寧1,王春武2,王東紅2

    (1.南京 航空航天大 學(xué)航空宇航 學(xué)院,江蘇 南 京 210016; 2.南京航空航天大學(xué)理學(xué)院,江蘇 南京210016)

    針對不可壓縮可壓縮水/氣多介質(zhì)問題,提出一種新的界面處理方法。在可壓縮水/氣界面處構(gòu)造Riemann問題,在水中設(shè)音速趨于無窮大,求解 Riemann問題得到不可壓縮可壓縮水/氣界面處流體的準(zhǔn)確流動(dòng)狀態(tài);然后以此狀態(tài)結(jié)合 GFM(ghost fluid method)方法分別為2種流體定義界面邊界條件,將兩相流 問題轉(zhuǎn)化為單相流問題計(jì)算,通過求解level set方程來跟蹤界面的位置。對各種不同的界面邊界條件定義方法進(jìn)行了比較,數(shù)值模擬結(jié)果表明算法能準(zhǔn)確地捕捉各類間斷的位置,證明了算法的有效性和穩(wěn)健性。

    爆炸力學(xué);Riemann問題;GFM 方法;level set;水/氣界面 ;多介質(zhì)流 動(dòng)

    在多介質(zhì)問題的數(shù)值模擬中,由于水/氣物質(zhì)屬性的巨大差異,使得水/氣多介質(zhì)問題成為多介質(zhì)問題的難點(diǎn)之一。針對具有大密度比的水/氣多 介質(zhì) 問 題,傳 統(tǒng) 方 法 主 要 有 3 類[1]:(1)水/氣 皆 被 作 為 可壓縮流體處理,其中水的狀態(tài)方程為 Tait方程;(2)水/氣皆作為不可壓縮流體,這類方法往往用來處理低速水/氣多介質(zhì)運(yùn)動(dòng);(3)水作為不可壓流體,氣作為可壓縮流體,相應(yīng)的分別采用不可壓方程和可壓縮歐拉方程作為水/氣的流場控制方程。對于第1類方法,由于界面兩邊的物質(zhì)屬性差異較大,過大的密度比和狀態(tài)方程差異,容易在界面處產(chǎn)生非物理震蕩,并且由于水的剛性很強(qiáng),聲速較大會(huì)導(dǎo)致計(jì)算時(shí)間步長過小,降低計(jì)算效率。第2類方法在處理氣體中有激波存在以及高速問題時(shí)會(huì)遇到困難,通常用來模擬低速多介質(zhì)運(yùn)動(dòng),如水中氣泡上升,靜止水滴下落等。采用第3類方法時(shí),充分考慮到水/氣物理屬性差異,可以處理氣體中有激波等高速運(yùn)動(dòng)問題,問題的關(guān)鍵在于如何給定合適的水/氣界面邊界條件,使得界面邊界條件能更合理的反映出水/氣界面處真實(shí)的流動(dòng)狀態(tài)。

    R.P.Fedkiw 等[2]提出了用 GFM 方法來處理多介質(zhì)問題界面邊界條件,最初只是用來處理氣氣多介質(zhì)問題,有效的 抑制了 在界面 處產(chǎn)生 的非物 理震蕩 。T.G.Liu等[3-4]在 其 基 礎(chǔ) 上 發(fā) 展 了 MGFM 方 法(modified ghost fluid method),通過在界面處定義 Riemann問題,并利用 Riemann問題的解定義了界面處虛擬流體點(diǎn)的速度和壓力值,通過等熵修正定義了真實(shí)流體點(diǎn)和虛擬流體點(diǎn)的密度值,MGFM 方法可以有效的處理各種氣氣多介質(zhì)問題,但是在處理類似激波阻尼的臨界問題時(shí)會(huì)遇到困難;C.W. Wang等[5-6]提出了 RGFM 方法(real ghost fluid method),利用 Riemann問題的解不 僅 定 義 了 虛 擬 流體點(diǎn)的值,而且對真實(shí)流體中點(diǎn)的密度,速度和壓力都進(jìn)行了更新,使得在界面處的流體狀態(tài)滿足接觸間斷的性質(zhì),得到了更加準(zhǔn)確的界面邊界條件。利用上述幾種GFM 方法在處理水/氣多介質(zhì)問題時(shí),把水作為可壓縮流體處理,狀態(tài)方程采用 Tait方程。而當(dāng)將水作為不可壓流體,氣體作為可壓縮流體處理時(shí),在界面處需要采用更加合理的界面邊界條件,才能滿足由于界面兩邊不同的控制方程和物質(zhì)屬性差異對界面處真實(shí)流動(dòng) 狀 態(tài) 的 影響。針 對 水/氣多 介 質(zhì) 問題,R.Caiden 等[1]提 出 了 new GFM 方法,界面處的速度采用界面附近水的速度值,而將界面附近氣體的壓力直接作為界面處的壓力。數(shù)值結(jié)果表明,該方法在一定程度上可以有效地捕捉物理現(xiàn)象。

    new GFM 方法給出的界面邊界條件過于簡單,并不能充分反應(yīng)出界面處真實(shí)的流動(dòng)狀態(tài),特別是當(dāng)界面兩邊水/氣的物質(zhì)屬性差異很大時(shí)。由于流體的流動(dòng)性質(zhì)很大程度上是物質(zhì)屬性相關(guān)的,因此界面處的流動(dòng)狀態(tài)通常非線性的依賴于界面兩邊流體的物質(zhì)屬性。在文獻(xiàn)[7-9]中有以下結(jié)論:在有限區(qū)域內(nèi),在馬赫數(shù)趨于零的情況下,無黏不可壓縮方程是歐拉方程的收斂極限。在這個(gè)結(jié)論的基礎(chǔ)上,考慮在界面處定義新的水/氣 Riemann問題,其中可以把水視為聲速趨于無窮大的可壓縮流體,并求解Riemann問題的解,從而得到界面處的真實(shí)流體狀態(tài)。在本文中,利用得到的 Riemann問題的解分別采用 MGFM 與RGFM 這2種方法定義水/氣界面邊界條件。通過一維數(shù)值模擬實(shí)驗(yàn),以及與利用new GFM 方法的計(jì)算結(jié)果比較分析,2種定義界面邊界條件的方法均能有效的捕捉到激波與界面的位置,同時(shí)分析3種方法的計(jì)算結(jié)果的差異及其原因。

    1 控制方程

    1.1 可壓縮流體控制方程

    一維無黏可壓縮流體控制方程為歐拉方程:

    式 中 :ρ為 密 度 ,u為 速 度 ,p為 壓 力 ,為 體 積 總 能 ,e為 質(zhì) 量 內(nèi) 能 。 理 想 氣 體 狀 態(tài) 方 程 為,其 中

    1.2 不可壓縮流體控制方程

    采 用 投 影 法 求 解 不 可 壓 控 制 方 程,對 式(2)關(guān) 于x方 向 求 偏 導(dǎo) 得 到,由 于 密 度 為 常 數(shù),則有。 表 明p為 常 數(shù),那 么 在 不 可 壓 流 體 區(qū) 域 內(nèi) 的 壓 力 滿 足 線 性 關(guān) 系 式,即 有,其 中xpr和pl分 別 為 不 可 壓 流 體 區(qū) 域 左 右 兩 個(gè) 端 點(diǎn) 上 的 的 壓 力 值,l為 不 可 壓 區(qū) 域 的 長 度 。

    1.3 Level set方 程

    通過求解level set方程隱式追蹤界面:

    Level set函 數(shù)初始 化為距 離符號 函數(shù),空間方向離散采用五階 WENO 格式[10],時(shí)間方向離散則采用三階 TVD Runge-Kutta方法。

    2 界面邊界條件

    在 new GFM 方法[1]中,只是簡單地直接取界面附近真實(shí)氣體的壓力作為界面處的壓力 值,界面的速度取臨近界面處水的速度值。事實(shí)上界面處的流動(dòng)狀態(tài)都是與物質(zhì)相關(guān)的,界面處的流動(dòng)狀態(tài)非線性的依賴于界面兩邊介質(zhì)的物質(zhì)屬性,由于水/氣的物質(zhì)屬性差異很大,應(yīng)考慮到其對界面處真實(shí)流動(dòng)狀態(tài)的影響。在 MGFM 和 RGFM 方法中,水作為可壓縮流體處理,采用剛性較強(qiáng)的 Tait方程作為其狀態(tài)方程,通過在界面處定義氣氣 Riemann問題,求解給出了界面處的流體狀態(tài),體現(xiàn)了界面兩邊水/氣的物質(zhì)屬性對于界面處流動(dòng)狀態(tài)的影響。但是考慮到水作為不可壓流體時(shí),具有密度不變等性質(zhì),不能求解氣氣 Riemann問題,上述方法則無法成立。當(dāng)將水作為不可壓流體處理時(shí),在文獻(xiàn)[9]中,有以下結(jié)論:隨著流體馬赫數(shù)趨向于零,可壓縮歐拉方程組的解會(huì)收斂到無黏不可壓流體方程組的解。另外,R.Agemi[7]證明了在有限區(qū)域內(nèi),當(dāng)馬赫數(shù)趨向于零時(shí),無 黏可壓縮流體會(huì) 收 斂 為相應(yīng)的無黏 不 可壓流體。基于以上結(jié)論,考慮將不可壓流體作為聲速無窮大的可壓縮流體處理,在這一基礎(chǔ)上,定義新的針對水/氣的 Riemann問題,并求解出 Riemann問題的解,利用其定義更合理的水/氣界面邊界條件。

    2.1 水/氣 Riemann問題

    在 水/氣 界 面 處 ,定 義 Riemann問 題 ,初 始 條 件ρl,ul,pl和ρr,ur,pr分 別 為 界 面 左 右 兩 側(cè) 的 密 度 ,速度和壓力。為不失一般性,假設(shè)左側(cè)為氣體,右側(cè)為水。根據(jù)以上結(jié)論:無黏不可壓流體是可壓縮流體在聲速趨向于無窮大時(shí)的收斂極限。在求解水/氣 Riemann問題時(shí)可以將水作為聲速趨向于無窮大的可壓流體處理。可壓縮流體中的 Riemann問題可以用雙激波近似方法求解,雙激波近似求解 Riemann問題有以下關(guān)系式:

    2.2 定義界面邊界條件

    在處理多介質(zhì)問題時(shí),GFM 通過定義虛擬流體點(diǎn)來給定合適的界面邊界條件。首先,在水/氣界面處建立求解水/氣 Riemann問題,利用得到的 Riemann問題的解,定義界面處真實(shí)流體點(diǎn)和虛擬流體點(diǎn)的流體狀態(tài)值,從而給定合適的界面邊界條件。在得到的水/氣 Riemann問題解的基礎(chǔ)上,本文中采用兩種方式定義界面邊界條件,分別為 MGFM 方法和RGFM 方法。

    如圖1所示,界面在網(wǎng)格點(diǎn)i與i+1之間,氣體在界面左側(cè),水在界面的右側(cè),取網(wǎng)格點(diǎn)和2的流場狀態(tài)分別作為 Riemann 問 題的 初 始 左 右 狀 態(tài)值,求 得 相 應(yīng) 的 解。接 觸 間 斷 左 側(cè) 狀 態(tài) 為,右 側(cè) 狀 態(tài) 為。 在 MGFM 方 法 中[3],定 義 氣 體 界 面 邊 界 條 件 ,利 用 得 到 的 Rie-mann問題的解左側(cè)狀態(tài) 定義網(wǎng)格點(diǎn)i+1上的虛擬流體點(diǎn)的速度和壓力,網(wǎng)格點(diǎn)i+2和i+3處的速度和壓力取當(dāng)?shù)卣鎸?shí)流體點(diǎn)的值,在求得界面左側(cè)熵值的基礎(chǔ)上利用等熵修正得到網(wǎng)格點(diǎn)i上真實(shí)流體 以 及 網(wǎng) 格 點(diǎn)i+1、i+2、i+3 上 出 虛 擬 流 體 的 密 度 值 ,上 述 則 給 定 了 氣 體 的 界 面 邊 界 條 件 。 可 以 采 用類似方法定義水的界面邊界條件,值得注意的是,Riemann問題的解中界面右側(cè)狀態(tài)中的密度也是保持不變的,滿足水的密度不變的性質(zhì),虛擬點(diǎn)上的密度都為水的密度。

    圖1 利用 MGFM方法對氣體定義界面邊界條件Fig.1 The definition of interface boundary condition for gas by MGFM method

    如 圖 2 所 示 ,在 RGFM 方 法 中[6],對 于 氣 體 ,得 到 的 Riemann 問 題 解 的 左 側(cè) 狀 態(tài)直接定義了網(wǎng)格點(diǎn)i+1、i+2和i+3上的虛擬流體狀態(tài)及網(wǎng) 格 點(diǎn)i上 的真實(shí)流體 狀 態(tài),上述即給定 了 氣體的界面邊界條件。對于水可進(jìn)行類似操作,同樣的,虛擬點(diǎn)上的密度都為水的密度。

    上述2種方法給出了水和氣體的界面邊界條件,在這個(gè)基礎(chǔ)上,分別采用不同的離散方法求解水和氣的控制方程,并通過求解level set方程更新level set函數(shù)隱式追蹤界面。

    圖2 利用RGFM方法對氣體定義界面邊界條件Fig.2 The definition of interface boundary condition for gas by RGFM method

    3 數(shù)值模擬

    對一維水/氣多介質(zhì)問題進(jìn)行數(shù)值模擬,時(shí)間方向上采用三階 TVD Runge-Kutta方法,空間上用三階ENO格式離散歐拉方程組,對于不可壓方程組采用標(biāo)準(zhǔn)的二階投影算法離散求解。

    3.1 水在空氣中高速運(yùn)動(dòng)

    在長度為1 m,網(wǎng)格點(diǎn)為200的區(qū)域內(nèi),中心處有長度為0.2 m 的水被周圍的靜止的氣體包圍,其中氣 體 的 狀 態(tài) 為γ=1.4,ρ=1.226 kg/m3,u=0,p=100 kPa,水 的 狀 態(tài) 為ρ=1 000 kg/m3,u= 100 m/s,p=100 k Pa。水在靜止的空氣中突然以一定速度向右運(yùn)動(dòng),導(dǎo)致水右側(cè)的氣體中產(chǎn)生激波,在水的左側(cè)氣體中產(chǎn)生稀疏波。圖3~6所示為在時(shí)間t=0.75 ms時(shí),利用 MGFM 和 RGFM 和 newGFM 這3種方法[1]計(jì)算得到的流場狀態(tài)結(jié)果對比圖。

    由圖3~6中可以看出,3種方法都能較準(zhǔn)確地捕捉到界面和激波的位置,在激波處,由于密度變化過小而很難在圖中體現(xiàn)。3種方法計(jì)算結(jié)果對于界面位置的捕捉結(jié)果基本相同,這是由于水/氣 Riemann問題的解給出的界面處的速度為界面附近水的速度,這一結(jié)論與new GFM 方法中界面的速度直接取界面附近水的速度相同。而 MGFM 方法和 RGFM 方法計(jì)算結(jié)果中,激波位置比new GFM 方法領(lǐng)先1個(gè)網(wǎng)格步長,這是由于界面處的壓力是由水/氣 Riemann問題得到,沒有簡單的取界面附近的氣體的壓力,充分考慮了界面兩側(cè)水/氣物質(zhì)屬性差異對界面處流動(dòng)狀態(tài)的非線性影響。此外,通過網(wǎng)格加密方法對 RGFM 方法進(jìn)行了收斂速度測試,在200、400 和 800 個(gè)網(wǎng)格點(diǎn)上計(jì)算得到在時(shí)間t= 0.75 ms時(shí)水的速度依次為99.693、99.689 和 99.687 m/s,計(jì)算得到收斂速度為0.992 7。

    圖3 密度1 000 kg/m3的水在空氣中向右運(yùn)動(dòng)時(shí)流場密度Fig.3 Density profile of the water movement in air while water density is 1 000 kg/m3

    圖4 密度1 000 kg/m3的水在空氣中向右運(yùn)動(dòng)時(shí)流場速度Fig.4 Velocity profile of the water movement in air while water density is 1 000 kg/m3

    圖5 密度1 000 kg/m3的水在空氣中向右運(yùn)動(dòng)時(shí)流場壓力Fig.5 Pressure profile of the water movement in air while water density is 1 000 kg/m3

    圖6 流場壓力細(xì)節(jié)對比圖Fig.6 Detail comparison of pressure profile

    為了進(jìn)一步考察水的密度對于計(jì)算結(jié)果的影響,將上述算例中的水密度變?yōu)?0 kg/m3,圖7~10給出了在時(shí)間t=0.75 ms時(shí),3種方法得到的計(jì)算結(jié)果對比。

    對于界面位置的捕捉,3種方法計(jì)算的結(jié)果基本相同,在激波處,密度變化明顯可以分辨。由于水的密度變小,更容易被空氣減慢速度,使得在水與右側(cè)的激波間,和水與左側(cè)的稀疏波之間都產(chǎn)生1個(gè)二次稀疏波。在這2個(gè)二次稀疏波中,new GFM 方法計(jì)算結(jié)果的速度大于 RGFM 和 MGFM 方法的計(jì)算結(jié)果,激波位置落后1個(gè)網(wǎng)格步長,MGFM 和 RGFM 方法只在稀疏波區(qū)域有細(xì)微差異,在區(qū)域其他點(diǎn)上基本相同。利用網(wǎng)格加密方法對 RGFM 方法進(jìn)行收斂速度測試,在200、400和800個(gè)網(wǎng)格點(diǎn)上計(jì)算得在t=0.75 ms時(shí),水的速度分別為73.589、73.315和73.171 m/s,計(jì)算出收斂速度為0.953 1。

    3.2 激波與水/氣界面相互作用

    圖7 密度10 kg/m3的水在空氣中向右運(yùn)動(dòng)時(shí)流場密度Fig.7 Density profile of the water movement in air while water density is 10 kg/m3

    圖8 密度10 kg/m3的水在空氣中向右運(yùn)動(dòng)時(shí)流場速度Fig.8 Velocity profile of the water movement in air while water density is 10 kg/m3

    圖9 密度10 kg/m3的水在空氣中向右運(yùn)動(dòng)時(shí)流場壓力Fig.9 Pressure profile of the water movement in airwhile water density is 10 kg/m3

    圖10 流場壓力細(xì)節(jié)對比圖Fig.10 Detail comparison of pressure profile

    在長度為1 m,網(wǎng)格數(shù) 為200 的區(qū)域 內(nèi),有 狀 態(tài) 為ρ=1 000 kg/m3,u=0,p=98.067 k Pa,長 度 為0.2 m的 水 在 區(qū) 域 中 心 ,周 圍 的 空 氣 狀 態(tài) 為ρ=1.58 kg/m3,u=0,p=98.067 kPa,在0.1 m 處 的 空 氣 中有一向右運(yùn)動(dòng)的激波,其 中波 后 狀 態(tài) 為ρ=2.124 kg/m3,u=89.98 m/s,p=148.407 k Pa。 圖 11~14所示為在時(shí)間t=1.75 ms時(shí),分別采用 MGFM、RGFM 和 new GFM 這3種方法的激波與水/氣界面相互作用的計(jì)算結(jié)果對比圖。激波與水/氣界面相互作用時(shí),由于水的剛性很強(qiáng)會(huì)產(chǎn)生1個(gè)反射激波,和1個(gè)很弱的入射波,入射波以極快的速度穿過水的區(qū)域,進(jìn)入右側(cè)的空氣,入射波由于強(qiáng)度太小而很難在圖中分辨。MGFM 和 RGFM 方法的計(jì)算結(jié)果中,激波的位置領(lǐng)先new GFM 方法中激波位置1個(gè)網(wǎng)格步長。由于入射波的作用,水滴產(chǎn)生1個(gè)向右的速度,同時(shí)由于水密度較大,導(dǎo)致產(chǎn)生的速度過小。通過網(wǎng)格加密方法對 RGFM 方法進(jìn)行收斂速度測試,在200、400和800個(gè)網(wǎng)格點(diǎn)上計(jì)算得到在時(shí)間t=1.75 ms時(shí)水的速度分別為0.531、0.538和0.542 m/s,計(jì)算得到收斂速度為0.9835。

    將上述算例中的水密度變?yōu)棣?10 kg/m3,圖15~18所示為采用3種方法計(jì)算得到結(jié)果對比圖。由結(jié)果可以看出,激波與水/氣界面相互作用以后,產(chǎn)生1個(gè)向左的反射激波,在這個(gè)反射激波與水之間形成了1個(gè)二次稀疏波,于此同時(shí),向右運(yùn)動(dòng)的入射波以極快的速度穿過水進(jìn)入右側(cè)的氣體中形成一個(gè)稀疏波,這一結(jié)果也符合波在水中的傳播速度無窮大的性質(zhì)。

    圖12 水密度1 000 kg/m3時(shí)激波與水/氣界面相互作用后的流場速度Fig.12 Velocity profile of shock impact with water-gas interface while water density is 1 000 kg/m3

    圖13 水密度1 000 kg/m3時(shí)激波與水/氣界面相互作用后的流場壓力Fig.13 Pressure profile of shock impact with water-gas interface while water density is 1 000 kg/m3

    圖14 流場壓力細(xì)節(jié)對比圖Fig.14 Detail comparison of pressure profile

    圖15 水密度10 kg/m3時(shí)激波與水/氣界面相互作用后的流場密度Fig.15 Density profile of shock impact with water-gas interface while water density is 10 kg/m3

    圖16 水密度10 kg/m3時(shí)激波與水/氣界面相互作用后的流場速度Fig.16 Velocity profile of shock impact with water-gas interface while water density is 10 kg/m3

    由圖15~18中的對比結(jié)果可知,在激波位置方面,MGFM 和 RGFM 方法的計(jì)算結(jié)果仍然比new GFM 方法的計(jì)算結(jié)果要領(lǐng)先,在二次稀疏波范圍內(nèi),new GFM 方法的計(jì)算結(jié)果中,速度比另兩種方法的結(jié)果稍小,而壓力比另外2種方法的結(jié)果大。

    由于 RGFM 方法中,不僅用 Riemann問題的解定義了界面處虛擬點(diǎn)的狀態(tài)值,同時(shí)也修改了界面處真實(shí)流體的速度和壓力,而 MGFM 方法中,只是修改了虛擬流體點(diǎn)的速度和壓力,對于界面處真實(shí)流體點(diǎn)只是通過等熵修正了密度,并沒有修改壓力和速度,因此在得到的結(jié)果中,RGFM 方法中水的速度相對另外兩種方法要略大。通過網(wǎng)格加密方法對 RGFM 方法進(jìn)行了收斂速度測試,在200、400和800個(gè)網(wǎng)格點(diǎn)上計(jì)算得到在時(shí)間t=1.75 ms時(shí)水的速度分別為40.117、40.431和40.589 m/s,計(jì)算得到收斂速度為0.9922。

    圖17 水密度10 kg/m3時(shí)激波與水/氣界面相互作用后的流場壓力Fig.17 Pressure profile of shock impact with water-gas interface while water density is 10 kg/m3

    圖18 流場壓力細(xì)節(jié)對比圖Fig.18 Detail comparison of pressure profile

    4 結(jié) 論

    針對水/氣多介質(zhì)問題,考慮將水作為不可壓流體,氣體作為可壓縮流體,在基于不可壓無黏流體是可壓縮流體在馬赫數(shù)趨向于零的收斂極限的結(jié)論上,定義并求解了水/氣 Riemann問題,Riemann問題的解反映出界面處真實(shí)的流動(dòng)狀態(tài),并采用 MGFM 以及RGFM 這2種方法給出了更合理的反應(yīng)水/氣物質(zhì)屬性差異的界面邊界條件。通過一維數(shù)值算例結(jié)果以及與相應(yīng)的new GFM 方法計(jì)算結(jié)果的對比分析,證明了算法的有效性及合理性,同時(shí)對 RGFM 方法進(jìn)行了收斂速度測試,充分考慮了水的不可壓縮性質(zhì)以及氣體的可壓縮性質(zhì),適合針對水/氣具有不同的運(yùn)動(dòng)特性的水/氣多介質(zhì)問題,例如空氣處于高速運(yùn)動(dòng)或者有強(qiáng)激波存在而水處于低速運(yùn)動(dòng)狀態(tài),此時(shí)不能被忽略水的不可壓性質(zhì),同樣的適合于需要考慮水中波的傳播速度等問題,如水下爆炸。

    [1]Caiden R,Fedkiw R P,Anderson C.A numerical method for two-phase flow consisting of separate compressible and incompressible regionsg[J].Journal of Computational Physics,2001,166(1):1-27.

    [2]Fedkiw R P,Aslam T,Merriman B,et al.A non-oscillatory Eulerian approach to interfaces in multimaterial flows:The ghost fluid method[J].Journal of Computational Physics,1999,152(2):457-492.

    [3]Liu T G,Khoo B C,Yeo K S.Ghost fluid method for strong shock impacting on material interface[J].Journal of Computational Physics,2003,190(2):651-681.

    [4]Liu T G,Khoo B C,Wang C W.The ghost fluid method for gas-water simulation[J].Journal of Computational Physics,2005,204(1):193-221.

    [5]Wang C W,Liu T G,Khoo B C.A real ghost fluid method for the simulation of multimedium compressible flow [J].SIAM Journal on Scientific Computing,2006,28(1):278-302.

    [6]王春 武,趙 寧.基于求 解 Riemann問題的界面 處理方法[J].計(jì)算物理 ,2006,22(4):306-310. Wang Chun-wu,Zhao Ning.An interface treating method based on Riemann problems[J].Chinese Journal of Computational Physics,2006,22(4):306-310.

    [7]Agemi R.The incompressible limit of compressible fluid motion in a bounded domain[C]∥Proceedings of the Japan Academy Series A:Mathematical Sciences.1981.

    [8]Steve S.The compressible Euler equations in a bounded domain:Existence of solutions and the incompressible limit [J].Communications in Mathematical Physics,1986,104(1):49-75.

    [9]Asano K.On the incompressible limit of the compressible Euler equation[J].Japan Journal of Applied Mathematics,1987,4(3):455-488.

    [10]Jiang G S,Shu C W.Efficient implementation of weighted ENO schemes[J].Journal of Computational Physics, 1996,126(1):202-228.

    Interface treating methods for the gas-water multi-phase flows

    Xu Shuang1,Zhao Ning1,Wang Chun-wu2,Wang Dong-hong2
    (1.College of Aerospace Engineering,Nanjing University of Aeronautics and Astronautics, Nanjing 210016,Jiangsu,China;
    2 College of Science,Nanjing University of Aeronautics and Astronautics, Nanjing 210016,Jiangsu,China)

    A new interface treating method is presented for the compressible-incompressible gas-water multi-phase flow.The Riemann problem is constructed at the compressible gas-water interface,and then solved according to the hypothesis that the sound speed tends to infinity in the water.The solution of Riemann problem provides the fluid states for compressible gas and incompressible water at the interface.Those states can then be used to define the interface boundary condition by coupling the ghost fluid method.The level set method is employed to track the interface.The numerical examples of one-dimension case are given in this paper,furthermore,several comparisons are made with other results to verify the algorithm.Numerical results show that the provided algorithm can capture the discontinuities accurately,which demonstrates the robustness and efficiency.

    mechanics of explosion;Riemann problem;ghost fluid method;level set;gas-water interface;multi-phase flows

    O382.1國標(biāo)學(xué)科代碼:13035

    :A

    10.11883/1001-1455-(2015)03-0326-09

    (責(zé)任編輯 王易難)

    2013-11-14;

    2014-02-28

    國家 自然科學(xué)基 金項(xiàng)目(11271188,91130030);北京理工大學(xué)爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目(KFJJ11-4 M)

    徐 爽(1984— ),男,博士研 究生,shuangxu@nuaa.edu.cn。

    猜你喜歡
    激波邊界條件流場
    大型空冷汽輪發(fā)電機(jī)轉(zhuǎn)子三維流場計(jì)算
    一類帶有Stieltjes積分邊界條件的分?jǐn)?shù)階微分方程邊值問題正解
    帶有積分邊界條件的奇異攝動(dòng)邊值問題的漸近解
    一種基于聚類分析的二維激波模式識別算法
    基于HIFiRE-2超燃發(fā)動(dòng)機(jī)內(nèi)流道的激波邊界層干擾分析
    斜激波入射V形鈍前緣溢流口激波干擾研究
    轉(zhuǎn)杯紡排雜區(qū)流場與排雜性能
    適于可壓縮多尺度流動(dòng)的緊致型激波捕捉格式
    基于HYCOM的斯里蘭卡南部海域溫、鹽、流場統(tǒng)計(jì)分析
    基于瞬態(tài)流場計(jì)算的滑動(dòng)軸承靜平衡位置求解
    无人区码免费观看不卡| 露出奶头的视频| 在线永久观看黄色视频| 国产在线观看jvid| 亚洲自拍偷在线| 国产xxxxx性猛交| 男女午夜视频在线观看| 一级a爱片免费观看的视频| 一区福利在线观看| 午夜免费观看网址| www.999成人在线观看| www.999成人在线观看| 久久久久久大精品| 日本免费a在线| 国产精品久久视频播放| 亚洲色图 男人天堂 中文字幕| 日本一区二区免费在线视频| 丝袜美足系列| 亚洲美女黄片视频| 亚洲九九香蕉| 成年女人毛片免费观看观看9| 精品国产国语对白av| 麻豆久久精品国产亚洲av | 制服诱惑二区| 人人妻人人添人人爽欧美一区卜| 黄色 视频免费看| 午夜免费观看网址| 香蕉国产在线看| 伦理电影免费视频| 国产99白浆流出| 国产成人系列免费观看| 久久99一区二区三区| 性欧美人与动物交配| 国产成人免费无遮挡视频| 色哟哟哟哟哟哟| 美女高潮喷水抽搐中文字幕| 国产精华一区二区三区| 免费在线观看日本一区| 韩国精品一区二区三区| 久久中文字幕一级| 国产精品亚洲av一区麻豆| 精品久久久精品久久久| 午夜福利欧美成人| 丝袜人妻中文字幕| 久久青草综合色| 丝袜在线中文字幕| 国产免费现黄频在线看| 亚洲人成伊人成综合网2020| 超碰成人久久| 国产不卡一卡二| 大型av网站在线播放| 大码成人一级视频| 亚洲精品一卡2卡三卡4卡5卡| 男人舔女人下体高潮全视频| 免费看十八禁软件| 多毛熟女@视频| 亚洲av五月六月丁香网| 亚洲,欧美精品.| 国产一区二区三区综合在线观看| tocl精华| 在线国产一区二区在线| 波多野结衣av一区二区av| 久久中文字幕一级| 婷婷丁香在线五月| 国产视频一区二区在线看| 亚洲激情在线av| 18禁裸乳无遮挡免费网站照片 | 亚洲一码二码三码区别大吗| 高清在线国产一区| 国产成人欧美| 视频区欧美日本亚洲| 热99re8久久精品国产| 看片在线看免费视频| 老司机靠b影院| 天堂√8在线中文| 午夜精品久久久久久毛片777| 夜夜躁狠狠躁天天躁| 欧美日韩黄片免| 老汉色av国产亚洲站长工具| 亚洲成a人片在线一区二区| 国产亚洲精品一区二区www| 色在线成人网| 色精品久久人妻99蜜桃| 欧美最黄视频在线播放免费 | 午夜免费激情av| 国产精品一区二区在线不卡| 三上悠亚av全集在线观看| 国产亚洲欧美在线一区二区| 国产成人欧美在线观看| 在线观看免费视频日本深夜| 岛国视频午夜一区免费看| 亚洲色图 男人天堂 中文字幕| 国产激情欧美一区二区| 无遮挡黄片免费观看| а√天堂www在线а√下载| 成人三级做爰电影| 丰满饥渴人妻一区二区三| 亚洲熟妇熟女久久| 首页视频小说图片口味搜索| 欧美成人午夜精品| 日本黄色日本黄色录像| 国产成人欧美在线观看| 一区二区三区激情视频| 黑丝袜美女国产一区| 亚洲免费av在线视频| 两性午夜刺激爽爽歪歪视频在线观看 | 天堂中文最新版在线下载| 婷婷丁香在线五月| 身体一侧抽搐| 午夜亚洲福利在线播放| 熟女少妇亚洲综合色aaa.| 国产精品1区2区在线观看.| 午夜免费观看网址| 日本vs欧美在线观看视频| 中文字幕精品免费在线观看视频| av在线播放免费不卡| 麻豆av在线久日| 黄色怎么调成土黄色| 纯流量卡能插随身wifi吗| 嫩草影院精品99| 亚洲成国产人片在线观看| 国产乱人伦免费视频| 亚洲激情在线av| 久9热在线精品视频| 热re99久久国产66热| 亚洲成人免费av在线播放| 国产av精品麻豆| 欧美人与性动交α欧美精品济南到| 国产亚洲欧美98| 国产av一区在线观看免费| 一区二区三区国产精品乱码| 黄网站色视频无遮挡免费观看| 在线视频色国产色| 国产一区二区激情短视频| 国产精品久久视频播放| 亚洲一区高清亚洲精品| 怎么达到女性高潮| 黄色成人免费大全| 久久人人精品亚洲av| 亚洲av成人av| 亚洲五月天丁香| 成人国语在线视频| 真人做人爱边吃奶动态| 麻豆成人av在线观看| 天天躁狠狠躁夜夜躁狠狠躁| 精品久久久久久久毛片微露脸| 欧美人与性动交α欧美精品济南到| 久久 成人 亚洲| 超碰97精品在线观看| 国产主播在线观看一区二区| 亚洲九九香蕉| 18禁黄网站禁片午夜丰满| 夜夜躁狠狠躁天天躁| 国产乱人伦免费视频| 新久久久久国产一级毛片| 88av欧美| 国产黄色免费在线视频| 婷婷丁香在线五月| 国产熟女xx| 12—13女人毛片做爰片一| 欧美av亚洲av综合av国产av| 99re在线观看精品视频| 亚洲精品中文字幕一二三四区| 久99久视频精品免费| 黑人操中国人逼视频| avwww免费| 午夜精品在线福利| 欧美精品啪啪一区二区三区| 美女高潮到喷水免费观看| 大码成人一级视频| 91字幕亚洲| 亚洲aⅴ乱码一区二区在线播放 | 亚洲av电影在线进入| 亚洲精品国产一区二区精华液| 丰满饥渴人妻一区二区三| 国产精品爽爽va在线观看网站 | a级毛片在线看网站| 人人妻人人添人人爽欧美一区卜| 色在线成人网| 99香蕉大伊视频| 不卡av一区二区三区| 黄片小视频在线播放| 性少妇av在线| 精品国产超薄肉色丝袜足j| 极品教师在线免费播放| xxxhd国产人妻xxx| 两个人免费观看高清视频| 视频区欧美日本亚洲| 人妻久久中文字幕网| 不卡一级毛片| 母亲3免费完整高清在线观看| 亚洲 欧美一区二区三区| 美女大奶头视频| 黄片大片在线免费观看| 色婷婷久久久亚洲欧美| 日日爽夜夜爽网站| 亚洲成人精品中文字幕电影 | 欧美性长视频在线观看| 亚洲中文av在线| 精品日产1卡2卡| 久久久久久免费高清国产稀缺| 99国产综合亚洲精品| 国产片内射在线| 韩国精品一区二区三区| 国产精品成人在线| 99久久99久久久精品蜜桃| 久久亚洲精品不卡| 欧美日韩一级在线毛片| 美女国产高潮福利片在线看| 国产精品自产拍在线观看55亚洲| 欧美久久黑人一区二区| 久久精品aⅴ一区二区三区四区| 久久香蕉激情| 激情在线观看视频在线高清| 日韩视频一区二区在线观看| 亚洲 欧美一区二区三区| 在线观看一区二区三区| 两人在一起打扑克的视频| 亚洲熟女毛片儿| 嫩草影院精品99| 久久影院123| 最好的美女福利视频网| 亚洲一卡2卡3卡4卡5卡精品中文| 国产亚洲av高清不卡| av在线天堂中文字幕 | avwww免费| 国产极品粉嫩免费观看在线| 夜夜看夜夜爽夜夜摸 | 真人做人爱边吃奶动态| 国产免费av片在线观看野外av| 一级作爱视频免费观看| 97人妻天天添夜夜摸| 亚洲精品美女久久av网站| 国产成人欧美在线观看| 女人被躁到高潮嗷嗷叫费观| 久久人妻福利社区极品人妻图片| 久久人人爽av亚洲精品天堂| 国产成人系列免费观看| 欧美精品一区二区免费开放| 国产熟女xx| 欧美+亚洲+日韩+国产| 99国产精品免费福利视频| 国产黄色免费在线视频| 在线av久久热| 一级毛片女人18水好多| 欧美乱色亚洲激情| 国产男靠女视频免费网站| avwww免费| 大陆偷拍与自拍| 日韩有码中文字幕| 中文字幕最新亚洲高清| 在线观看免费视频日本深夜| 国产免费av片在线观看野外av| 国产精品久久视频播放| 最近最新中文字幕大全电影3 | 日本a在线网址| 亚洲精品国产一区二区精华液| 免费日韩欧美在线观看| 国产又色又爽无遮挡免费看| 在线观看午夜福利视频| 精品久久久久久久久久免费视频 | 国产一区二区三区在线臀色熟女 | 午夜两性在线视频| 女生性感内裤真人,穿戴方法视频| 国产黄色免费在线视频| 亚洲色图综合在线观看| 亚洲 国产 在线| 伦理电影免费视频| 咕卡用的链子| 亚洲专区中文字幕在线| 亚洲成av片中文字幕在线观看| 成年人黄色毛片网站| 久久 成人 亚洲| 午夜免费成人在线视频| 久久人妻福利社区极品人妻图片| 欧美丝袜亚洲另类 | 高清黄色对白视频在线免费看| 亚洲avbb在线观看| 欧美日韩av久久| 亚洲成人国产一区在线观看| 一区在线观看完整版| 国产99久久九九免费精品| 99国产极品粉嫩在线观看| 欧美在线一区亚洲| 中文字幕精品免费在线观看视频| 大陆偷拍与自拍| 国产精品电影一区二区三区| 老汉色av国产亚洲站长工具| 黄色a级毛片大全视频| 欧美黑人精品巨大| 一级片'在线观看视频| 免费看a级黄色片| 欧美日韩一级在线毛片| 天堂俺去俺来也www色官网| 男女下面插进去视频免费观看| 午夜91福利影院| 黄色视频不卡| 日韩大尺度精品在线看网址 | 黄片大片在线免费观看| 亚洲人成网站在线播放欧美日韩| 99国产精品一区二区蜜桃av| 热99国产精品久久久久久7| 免费高清在线观看日韩| 一区二区三区国产精品乱码| 国产一区二区三区综合在线观看| 精品第一国产精品| 日韩欧美三级三区| 久久人妻av系列| 婷婷六月久久综合丁香| 亚洲avbb在线观看| 亚洲精品久久成人aⅴ小说| 狂野欧美激情性xxxx| 人人妻,人人澡人人爽秒播| 国产熟女午夜一区二区三区| 中文亚洲av片在线观看爽| 亚洲aⅴ乱码一区二区在线播放 | 亚洲一卡2卡3卡4卡5卡精品中文| 久久婷婷成人综合色麻豆| 18禁裸乳无遮挡免费网站照片 | 日韩欧美免费精品| 国产伦一二天堂av在线观看| 午夜福利在线免费观看网站| 亚洲中文日韩欧美视频| 美女高潮到喷水免费观看| 黄色怎么调成土黄色| 黄片小视频在线播放| 色尼玛亚洲综合影院| 一a级毛片在线观看| 久久精品国产99精品国产亚洲性色 | 国产乱人伦免费视频| 亚洲国产精品一区二区三区在线| 日韩精品中文字幕看吧| 亚洲狠狠婷婷综合久久图片| 色播在线永久视频| 91成年电影在线观看| 黄色毛片三级朝国网站| 一二三四社区在线视频社区8| 久久国产精品人妻蜜桃| 男女床上黄色一级片免费看| 日韩大尺度精品在线看网址 | 99在线视频只有这里精品首页| 天堂俺去俺来也www色官网| 免费看十八禁软件| 成人特级黄色片久久久久久久| 这个男人来自地球电影免费观看| 日韩中文字幕欧美一区二区| 一个人观看的视频www高清免费观看 | 免费一级毛片在线播放高清视频 | 黄色片一级片一级黄色片| 国产一区二区三区视频了| 国产精品偷伦视频观看了| 亚洲成人免费电影在线观看| 在线观看免费日韩欧美大片| 一区二区三区激情视频| 国产欧美日韩一区二区精品| 成人av一区二区三区在线看| av福利片在线| www.自偷自拍.com| 桃色一区二区三区在线观看| 丁香欧美五月| 一级片'在线观看视频| 美女午夜性视频免费| 婷婷丁香在线五月| 一区二区日韩欧美中文字幕| 国产人伦9x9x在线观看| 88av欧美| 国产成人免费无遮挡视频| 女生性感内裤真人,穿戴方法视频| 啦啦啦在线免费观看视频4| 国产黄a三级三级三级人| av国产精品久久久久影院| 久久久水蜜桃国产精品网| 欧美最黄视频在线播放免费 | 成人18禁高潮啪啪吃奶动态图| 亚洲激情在线av| 亚洲精品国产色婷婷电影| 国产精品久久久人人做人人爽| 国产精品永久免费网站| 亚洲av美国av| 国产精品野战在线观看 | 曰老女人黄片| 这个男人来自地球电影免费观看| 国产无遮挡羞羞视频在线观看| 国产色视频综合| 亚洲 欧美 日韩 在线 免费| 精品熟女少妇八av免费久了| 亚洲伊人色综图| 极品人妻少妇av视频| 99久久人妻综合| 极品人妻少妇av视频| 国产三级在线视频| 欧美一级毛片孕妇| 欧美在线黄色| 黑人猛操日本美女一级片| 在线永久观看黄色视频| 50天的宝宝边吃奶边哭怎么回事| 国产精品日韩av在线免费观看 | 国产av在哪里看| 亚洲一区高清亚洲精品| 一a级毛片在线观看| 日韩欧美在线二视频| 丰满的人妻完整版| 日韩大尺度精品在线看网址 | 精品一区二区三区四区五区乱码| 亚洲欧美日韩无卡精品| 欧美另类亚洲清纯唯美| 亚洲第一青青草原| 一级毛片精品| 大陆偷拍与自拍| 老司机深夜福利视频在线观看| 狠狠狠狠99中文字幕| a级毛片黄视频| 中亚洲国语对白在线视频| 变态另类成人亚洲欧美熟女 | 午夜免费鲁丝| 国产成人精品在线电影| 国产精品久久久久成人av| 男女下面插进去视频免费观看| 国产精品1区2区在线观看.| 18禁裸乳无遮挡免费网站照片 | 亚洲 欧美一区二区三区| 琪琪午夜伦伦电影理论片6080| 精品卡一卡二卡四卡免费| 少妇裸体淫交视频免费看高清 | 叶爱在线成人免费视频播放| 手机成人av网站| 日韩高清综合在线| 国产欧美日韩一区二区三区在线| 免费不卡黄色视频| 欧美在线黄色| 国产精品一区二区在线不卡| 国产精品九九99| 精品国产国语对白av| 51午夜福利影视在线观看| 777久久人妻少妇嫩草av网站| 午夜免费激情av| 欧美日韩亚洲高清精品| 99久久精品国产亚洲精品| 国产真人三级小视频在线观看| 久久伊人香网站| 男女床上黄色一级片免费看| 国产亚洲精品综合一区在线观看 | 动漫黄色视频在线观看| 午夜老司机福利片| 精品国产一区二区三区四区第35| 丰满迷人的少妇在线观看| 欧美激情极品国产一区二区三区| 一本综合久久免费| 视频在线观看一区二区三区| 国产高清视频在线播放一区| 两人在一起打扑克的视频| 国产又色又爽无遮挡免费看| 国产在线精品亚洲第一网站| 国产一区二区三区在线臀色熟女 | 欧美精品一区二区免费开放| 美女国产高潮福利片在线看| av网站免费在线观看视频| 天堂俺去俺来也www色官网| 在线观看日韩欧美| 丁香欧美五月| xxx96com| 亚洲狠狠婷婷综合久久图片| 精品福利观看| 久久久国产成人免费| 免费看a级黄色片| 亚洲男人的天堂狠狠| 1024香蕉在线观看| 男男h啪啪无遮挡| 欧美日韩乱码在线| 黑人巨大精品欧美一区二区mp4| 中文字幕人妻熟女乱码| 国产在线观看jvid| 久久人妻av系列| 国产免费现黄频在线看| 手机成人av网站| 制服诱惑二区| 无遮挡黄片免费观看| 男女高潮啪啪啪动态图| 欧美激情极品国产一区二区三区| 欧美性长视频在线观看| 在线观看免费日韩欧美大片| 国产av在哪里看| 夫妻午夜视频| 女生性感内裤真人,穿戴方法视频| 欧美乱色亚洲激情| 桃红色精品国产亚洲av| 午夜免费成人在线视频| 免费av中文字幕在线| 久久精品影院6| 高清欧美精品videossex| 免费看十八禁软件| 波多野结衣高清无吗| 亚洲avbb在线观看| 在线天堂中文资源库| 中文字幕精品免费在线观看视频| 超色免费av| 乱人伦中国视频| 日韩欧美在线二视频| 99在线人妻在线中文字幕| 美国免费a级毛片| 免费看十八禁软件| 久久久久亚洲av毛片大全| 欧美日韩亚洲国产一区二区在线观看| 一二三四社区在线视频社区8| 国产伦一二天堂av在线观看| 嫁个100分男人电影在线观看| 男人舔女人的私密视频| 国产激情久久老熟女| 天堂√8在线中文| 超碰成人久久| tocl精华| 免费少妇av软件| 久久久久久久久免费视频了| 久久中文字幕一级| 日韩中文字幕欧美一区二区| 亚洲免费av在线视频| 国产精品久久久人人做人人爽| 久久国产亚洲av麻豆专区| 国产成人精品无人区| 麻豆久久精品国产亚洲av | 精品乱码久久久久久99久播| 女人爽到高潮嗷嗷叫在线视频| 1024视频免费在线观看| 久久久久亚洲av毛片大全| 免费av中文字幕在线| 久久久久久久午夜电影 | 久久香蕉激情| 久久狼人影院| 国产精品野战在线观看 | 悠悠久久av| 婷婷丁香在线五月| 国产一区二区三区综合在线观看| 又黄又爽又免费观看的视频| 欧美日韩福利视频一区二区| 男女下面进入的视频免费午夜 | 国产欧美日韩综合在线一区二区| 精品一品国产午夜福利视频| 99国产精品99久久久久| 最近最新免费中文字幕在线| 乱人伦中国视频| 国产成人啪精品午夜网站| 亚洲av美国av| 国产精品乱码一区二三区的特点 | 亚洲 欧美一区二区三区| 满18在线观看网站| 悠悠久久av| 亚洲男人的天堂狠狠| 女性生殖器流出的白浆| 国产精品免费视频内射| 色婷婷av一区二区三区视频| 精品乱码久久久久久99久播| 欧美 亚洲 国产 日韩一| 成人手机av| 在线视频色国产色| 脱女人内裤的视频| 母亲3免费完整高清在线观看| 自线自在国产av| 亚洲精品国产色婷婷电影| 99国产综合亚洲精品| 精品国产亚洲在线| 女人被狂操c到高潮| 91精品三级在线观看| 欧美激情久久久久久爽电影 | 丝袜美足系列| 亚洲人成电影观看| 波多野结衣av一区二区av| 一级毛片女人18水好多| 夜夜爽天天搞| 乱人伦中国视频| 丰满迷人的少妇在线观看| 18禁美女被吸乳视频| 亚洲成人久久性| 校园春色视频在线观看| 五月开心婷婷网| 九色亚洲精品在线播放| av片东京热男人的天堂| 男女床上黄色一级片免费看| 色老头精品视频在线观看| 国产成人av教育| 18禁国产床啪视频网站| 日韩国内少妇激情av| 99在线视频只有这里精品首页| 国产主播在线观看一区二区| netflix在线观看网站| 国产色视频综合| 国产精品香港三级国产av潘金莲| 中文字幕高清在线视频| 在线观看午夜福利视频| 亚洲国产精品sss在线观看 | 午夜免费观看网址| 超碰成人久久| 99久久国产精品久久久| 午夜免费观看网址| 欧美 亚洲 国产 日韩一| 在线观看免费午夜福利视频| 久久天堂一区二区三区四区| 99国产精品一区二区三区| 日本一区二区免费在线视频| 亚洲精品一二三| 国产精品久久久久成人av| 99热国产这里只有精品6| 精品国产美女av久久久久小说| 亚洲自偷自拍图片 自拍| 操出白浆在线播放| 宅男免费午夜| 久久久久国产精品人妻aⅴ院| 操出白浆在线播放| a级片在线免费高清观看视频| 老司机深夜福利视频在线观看| 国产精品美女特级片免费视频播放器 | 中文字幕人妻丝袜一区二区| 女性被躁到高潮视频| 亚洲人成伊人成综合网2020| 十八禁网站免费在线| 少妇 在线观看|