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

    一種基于雙界面函數(shù)的界面捕捉方法1)

    2017-12-18 13:23:48何志偉駱龍山田保林2
    力學學報 2017年6期
    關(guān)鍵詞:正弦重構(gòu)數(shù)值

    李 康 劉 娜 何志偉 駱龍山 田保林2)

    *(北京應用物理與計算數(shù)學研究所,北京100088)

    ?(中物院高性能數(shù)值模擬軟件中心,北京100088)

    一種基于雙界面函數(shù)的界面捕捉方法1)

    李 康*,?劉 娜*,?何志偉*駱龍山*,?田保林*2)

    *(北京應用物理與計算數(shù)學研究所,北京100088)

    ?(中物院高性能數(shù)值模擬軟件中心,北京100088)

    基于代數(shù)重構(gòu)思想,發(fā)展了一種新的雙界面函數(shù)重構(gòu)方法,并采用雙正弦函數(shù)構(gòu)造了雙正弦界面重構(gòu)方法(double sine interface capturing,DSINC).為驗證不同界面函數(shù)對界面捕捉效果的影響,用數(shù)值方法求解了可壓縮五方程模型,其中對流項的離散采用五階WENO(weighted essentially non-oscillatory method)格式,時間積分采用三階Runge--Kutta方法,通量計算分別考慮了HLL和HLLC方法,而狀態(tài)方程采用Mie-Grneisen狀態(tài)方程.在數(shù)值計算中,在界面附近,采用DSINC來獲得體積分數(shù)的重構(gòu),而在遠離界面的區(qū)域采用WENO格式來獲得高階插值狀態(tài).相比采用單界面函數(shù)的方法,如雙曲正切界面重構(gòu)方法(tangent of hyperbola for interface capturing,THINC),DSINC方法同樣具有界面重構(gòu)算法簡單,在程序中添加方便等特點,兩者區(qū)別在于,DSINC方法在重構(gòu)過程中未知函數(shù)更易于求解,而無需求解復雜的非線性超越方程,這就使其具有易于向多維擴展的能力.一些典型的兩相流動問題,如圓形水柱對流問題,兩相三波點問題和激波--界面不穩(wěn)定性問題等被用作不同界面函數(shù)對界面捕捉效果的影響對比.對比分析發(fā)現(xiàn),DSINC與THINC在界面捕捉效果上大致保持一致,并在計算中表現(xiàn)出了較好的穩(wěn)定性.雙界面函數(shù)重構(gòu)思想可以為多相流動界面的代數(shù)重構(gòu)提供了一種新的思路.

    多相流動,界面捕捉,雙曲正切界面重構(gòu)方法,雙界面函數(shù)重構(gòu),雙正弦界面重構(gòu)方法

    引言

    界面分割不同流體的研究在許多科學和工程領(lǐng)域都有重要的應用價值,如超燃沖壓發(fā)動機、天體物理Rayleigh-Taylor(RT)不穩(wěn)定性和慣性約束核聚變等[1].在界面問題的數(shù)值模擬過程中,界面的演變問題受到格外的關(guān)注.如何更精確捕捉界面特征,降低數(shù)值方法帶來的界面耗散以及更真實反映流動特征,一直是當前界面研究中重點關(guān)注的問題.

    從方法類型上來區(qū)分,界面捕捉方法有兩種,一是界面追蹤方法(interface tracking method,IT),另一個是界面捕捉方法(interface capturing method,IC).在界面追蹤方法中,界面位置可以顯式給出,如陣面追蹤方法 (front-tracking method)[2-5]和標記法 (marker method)[6-7].然而,界面追蹤方法雖可以有效地標識界面的位置,但在處理大變形界面和拓撲變化界面時存在困難.在界面捕捉方法中,界面位置通過隱式方法給出,如水平集方法(level-set method,LS)[8-9]和流體體積分數(shù)方法(volume-of- fl uid method,VOF)[10].LS方法通過引入一個帶符號的距離函數(shù)來確定界面位置,界面位于函數(shù)等于零處,所捕捉的界面十分陡峭,但其不足之處是不滿足離散守恒(體積守恒或質(zhì)量守恒).VOF方法采用體積分數(shù)來分辨界面的位置,該方法最大的優(yōu)點是可以保證體積守恒,但由于耗散性,界面的分辨率不高,即使采用高精度方法也無法完全去除界面耗散[11].Sussman等[12]結(jié)合LS和VOF方法的優(yōu)點,采用耦合的LS-VOF方法研究多相流動界面問題時,界面模擬的可信度明顯提高.So等[13]將反擴散 (anti-di ff usion)方法用于多相流的求解,通過引入反擴散項來降低界面的耗散,隨后將反擴散方法拓展至可壓縮多相流動的數(shù)值模擬中[14].這些數(shù)值模擬技術(shù)的提出在一定程度上能夠有效解決界面問題研究中的一些問題.

    除上述數(shù)值模擬技術(shù)外,界面重構(gòu)技術(shù)在界面問題的研究中也占有重要地位.從構(gòu)造過程來看,界面重構(gòu)分為代數(shù)重構(gòu)和幾何重構(gòu).Youngs等[15]提出了分段線性界面重構(gòu)(piecewise linear interface calculation,PLIC)的方法,采用幾何重構(gòu)方法通過分段的線性函數(shù)來重構(gòu)界面.而近年來,代數(shù)重構(gòu)方法在獲得清晰界面方面同樣受到了廣泛關(guān)注.Xiao等[16-17]采用雙曲正切函數(shù)重構(gòu)界面,從而使得數(shù)值模擬中的界面更為陡峭,由于采用的雙曲正切重構(gòu)(tangent of hyperbola interface capturing,THINC)方法不同于 PLIC幾何重構(gòu)方法[15]的復雜過程,THINC過程簡單,在程序中添加方便.從構(gòu)造思路來看,THINC方法的主要針對一維問題,在處理實際的多維問題時需要采用分裂方法或者多次積分的方法.Li等[18]考慮到采用分裂的方法將一維THINC方法拓展至多維情況可能會降低界面模擬的可信度,將THINC方法通過多次積分的方法拓展至多維情況.針對其中的雙曲正切函數(shù)多次積分的困難問題,其采用了高斯積分克服了困難.考慮不同的界面函數(shù)會對界面捕捉效果產(chǎn)生不同的影響,Cassidy等[19]分別采用分段的直線、正弦函數(shù)作為界面函數(shù)來構(gòu)造界面得到分段線性界面方法(piecewise linear interface capturing,PLINC)、分段正弦界面方法(piecewise sine interface capturing,PSINC),結(jié)果表明:不同的界面函數(shù)(線性或正弦)均可以獲得較為陡峭的界面,但針對不同的問題時略有差異.對以上代數(shù)界面重構(gòu)方法的對比發(fā)現(xiàn),這些方法均采用單界面函數(shù)來重構(gòu)界面,因此本文中稱之為單界面函數(shù)代數(shù)重構(gòu)方法.

    對于單界面函數(shù)重構(gòu)方法來說,為了保證單元體內(nèi)的質(zhì)量守恒,即使界面重構(gòu)函數(shù)較為簡單,在實際構(gòu)造過程中由于需要對界面函數(shù)進行守恒積分,單界面函數(shù)在守恒積分時的求解過程也甚為復雜.為了進一步研究界面重構(gòu)函數(shù)對數(shù)值模擬結(jié)果的差異以及進一步評估多維界面重構(gòu)時的簡易方法,本文提出一種基于雙界面函數(shù)的重構(gòu)方法,并采用雙界面重構(gòu)思想構(gòu)造了雙正弦界面重構(gòu)方法,對比了不同類型界面重構(gòu)函數(shù)對數(shù)值結(jié)果的影響差別.本文的組織結(jié)構(gòu)如下.第1部分介紹所采用的數(shù)值方法.第2部分詳細給出雙界面重構(gòu)的思想和雙正弦界面函數(shù)構(gòu)造方法.數(shù)值計算結(jié)果在第3部分給出.第4部分介紹研究結(jié)論.

    1 數(shù)值方法

    1.1 控制方程:可壓縮多相流五方程模型

    對于不可溶、無黏的可壓縮兩相流動來說,流動過程可用五方程來描述

    式中,ρ為混合密度,ρi為第i相的密度,u=(u,v)為流動速度,p為壓力,E為總能,αi為第i相的體積分數(shù).五方程模型除了包含質(zhì)量、動量和能量三大守恒方程外,還包括組分質(zhì)量守恒和體積分數(shù)演化方程.界面重構(gòu)后,通過對流場中的密度和總能做相應的修正,保證在界面處的質(zhì)量、動量和能量守恒[11].對流項的數(shù)值離散采用常用的五階WENO格式和HLL或HLLC方法求解[20-21].

    在等壓假設(shè)下,方程(2)可以導出

    在式 (2)和式 (3)中,ek為分相k的內(nèi)能,Γk(ρk)=(1/ρk)(?pk/?ek)|ρk,p∞,k和e∞,k定義壓力和內(nèi)能的參考狀態(tài).

    1.2 原始的THINC界面重構(gòu)方法

    以一維情況為例,先介紹THINC的構(gòu)造思想[16],然后以此為基礎(chǔ),進一步介紹采用雙界面函數(shù)的界面重構(gòu)方法的構(gòu)造思想.設(shè)為t時刻單元體i的平均體積分數(shù),αi(x,t)為t時刻體積分數(shù)在單元體i中的分布.因此,對于單元體i來說,即x∈ [xi?1/2,xi+1/2],有

    式中,?xi=xi+1/2?xi?1/2為單元體的網(wǎng)格尺寸.根據(jù)體積分數(shù)的定義,有

    在實際計算中,需要定義界面單元,這時體積分數(shù)的取值范圍滿足

    式中,ε為一個小量,可取為10?3.條件(a)定義了界面處體積分數(shù)的取值區(qū)間,條件(b)定義了界面的單調(diào)性.

    在界面處體積分數(shù)從0變化到1,THINC方法選用雙曲正切函數(shù)來模擬界面的變化趨勢,如圖1所示.對于固定的時間t和空間單元x∈[xi?1/2,xi+1/2],有從0到1的曲線

    圖1 雙曲函數(shù)的變化規(guī)律Fig.1 Hyperbolic tangent function with x

    求解方程(8)可得

    于是可知,包括體積分數(shù)重構(gòu)在內(nèi),任何物理量的重構(gòu)過程如下

    式中

    2 基于雙界面函數(shù)的界面重構(gòu)方法

    由上述可知,THINC方法的構(gòu)造過程保證了重構(gòu)后體積分數(shù)的守恒性,但一次積分后雙曲正切函數(shù)變得較為復雜.本文的研究在考察不同界面函數(shù)對界面重構(gòu)效果的同時,也提供一種向多維擴展的可能性.與單界面函數(shù)不同,雙界面函數(shù)在一個子網(wǎng)格內(nèi)采用兩個函數(shù)來重構(gòu)界面.對于雙界面重構(gòu)方法,多種界面函數(shù)可供選擇,本文以雙正弦界面函數(shù)的界面重構(gòu)方法DSINC為例來詳細介紹雙界面函數(shù)的構(gòu)造過程.如圖2所示,DSINC方法采用兩條連續(xù)的正弦函數(shù)來重構(gòu)界面,兩個函數(shù)可以存在一階間斷,間斷點位于單元體的中心.圖2中,HL(x)定義為左正弦函數(shù),HR(x)定義為右正弦函數(shù),Cj為一階間斷點.下面給出界面重構(gòu)后左右狀態(tài)的推導過程.

    圖2 DSINC界面重構(gòu)示意圖Fig.2 The illustration of interface functions used by DSINC

    由于HL(x)和HR(x)連結(jié)了0和1,易知其函數(shù)形式為

    若設(shè)定i為單元體的中心索引,在保證界面函數(shù)連續(xù)的情況下,對于上升界面,即當時,則

    重構(gòu)后的界面需保證體積分數(shù)守恒,即有

    將方程(12)和方程(13)代入方程(14)可以求得連結(jié)點Ci的位置.經(jīng)求解可得,上升界面(σi>0)與下降界面(σi<0)具有相同的積分結(jié)果

    于是可知,采用DSINC方法重構(gòu)界面后,左右狀態(tài)分別為:

    當σi>0時,有

    需要指出的是,由Cj的意義可知,0

    此時,界面重構(gòu)后的左右狀態(tài)分別為

    此時,界面重構(gòu)后的左右狀態(tài)分別為

    此時,界面重構(gòu)后的左右狀態(tài)分別為

    此時,界面重構(gòu)后的左右狀態(tài)分別為

    基于以上重構(gòu)思想,同樣可以得到雙線性界面重構(gòu)方法 (double linear interface capturing,DLINC).對于單界面函數(shù)重構(gòu)方法和雙界面函數(shù)重構(gòu)方法,圖 3給出了單界面函數(shù)以及雙界面函數(shù)在不同單元體中心體積分數(shù)值時的分布規(guī)律,其中ξ=(x?xi?1/2)/(xi+1/2?xi?1/2),而αcell為計算所得單元體中心值.對于THINC方法來說,反映界面斜率的參數(shù)β取2.3.從界面函數(shù)的分布規(guī)律來看,PSINC方法的界面最為陡峭,PLINC次之,THINC最為平緩,而DSINC的界面函數(shù)平緩度位于PLINC和THINC之間,界面函數(shù)的特性在一定程度上決定了界面捕捉的效果.

    圖3 界面重構(gòu)后界面函數(shù)α隨空間位置的變化Fig.3 The distribution of interface function in a cell under di ff erent cell value

    方法測試的數(shù)值計算基于作者所在課題組開發(fā)的可壓縮流體有限差分程序(code of fi nite di ff erence forcompressible fl uiddynamics,CFD2)開展.CFD2程序包含多種高階方法,如WENO,MP和MUSCL等,通量的求解可選用Steger-Warming(SW),Lax-Friedrichs(LF)或者HLL黎曼求解器.考慮到本文所研究的主要內(nèi)容為數(shù)值方法的對比,選用五階WENO方法來求解左、右特征狀態(tài),用HLL方法來求解通量.在界面方法的驗證中,在界面位置,WENO方法用THINC或DSINC方法來代替,流場中其他位置的求解保持不變.

    3 數(shù)值結(jié)果與分析

    3.1 圓形水柱對流問題

    為了研究所構(gòu)造的界面方法對界面的捕捉效果,首先模擬了典型的圓形水柱的對流問題.計算區(qū)域(x,y)∈[0,2]m×[0,2]m,圓形水柱的初始位置為(0.5,0.5)m,水柱的半徑R為0.25m,對流速度設(shè)為常數(shù)(ux,uy)=(100,100)m/s.此問題中,水的密度為ρ1=1000kg/m3,周圍空氣的密度為ρ2=1kg/m3.對于狀態(tài)方程,參數(shù)選取為:對于水相,γ1=4.4,c01=1624.8m/s,ρ01=1000kg/m3;對于氣相,γ2=1.4,c02=0,ρ02=1kg/m3.邊界條件設(shè)為周期邊界.在計算中,x,y方向的網(wǎng)格數(shù)均為100,計算CFL數(shù)取為0.2.

    圖4~圖6分別給出了計算300時間步后無界面重構(gòu)方法、THINC方法和DSINC方法所得水柱界面體積分數(shù)分布.在不采用界面重構(gòu)方法時,界面隨著計算時間的延長而逐漸變寬,即使采用高階格式,界面依然在逐漸耗散,如圖4所示.對于THINC格式來說,無論是低階還是高階格式,界面均可以保持一種銳利(sharp)的效果.圖6所給出的DSINC的計算結(jié)果,在低階時,界面可以保持銳利,而高階時界面捕捉效果更好.

    圖4 無界面重構(gòu)時圓柱對流數(shù)值模擬結(jié)果Fig.4 Numerical results of water column problems without interface reconstruction

    圖5 THINC界面重構(gòu)下圓柱對流數(shù)值模擬結(jié)果Fig.5 Numerical results of water column problems with THINC

    圖6 DSINC界面重構(gòu)下界面重構(gòu)數(shù)值結(jié)果Fig.6 Numerical results of water column problems with DSINC

    3.2 兩相三波點問題

    兩相三波點問題的計算區(qū)域為(x,y)∈[0,7]×[0,3],計算網(wǎng)格數(shù)為210×90,其他計算參數(shù)依照文獻[11]來設(shè)置.圖7給出了一階精度下流場中的密度和體積分數(shù)分布.結(jié)果顯示,在無界面重構(gòu)方法時,流場中密度和體積分數(shù)分布逐漸變寬,以至于界面變得十分模糊.在采用界面重構(gòu)方法時,對于體積分數(shù)來說,DSINC和THINC方法均能獲得較為清晰的界面.同時,對比密度和體積分數(shù)的分布表明,體積分數(shù)的耗散來源于非守恒方程,而密度的耗散來自于連續(xù)方程,即使計算中采用體積分數(shù)去修正了流場中的狀態(tài)分布,密度的耗散依然十分嚴重.在高階格式下,所得密度和體積分數(shù)分布如圖8(a)和圖8(b)所示.而對于體積分數(shù)分布來說,采用高階格式可以獲得較為清晰的界面,但界面依然在逐漸耗散.這一結(jié)論與圓形水柱對流的計算結(jié)果較為一致.HLL通量方法本身具有一定的耗散性,為了分析這種特性,圖8(c)和圖8(d)給出了采用HLLC方法模擬三波點問題的數(shù)值結(jié)果.結(jié)果顯示,對于所模擬的問題來說,結(jié)果僅在半圓形的界面處存在有限的改善.這表明,在兩相流問題的數(shù)值模擬中,獲得清晰的界面捕捉效果仍需對界面采用界面捕捉格式進行處理.

    采用五階精度數(shù)值格式和THINC/DSINC界面方法來模擬的兩相三波點問題如圖9所示.在高階格式下,對比圖7發(fā)現(xiàn),采用DSINC和THINC依然均可以使得界面保持銳利的捕捉效果.進一步分析發(fā)現(xiàn),相比THINC方法,DSINC方法略帶耗散.Xiao等[16]分析了界面函數(shù)變化趨勢以及在界面函數(shù)下所產(chǎn)生的耗散性,依據(jù)其分析結(jié)果可知,DSINC所產(chǎn)生的略大的界面耗散性是由界面函數(shù)的特點所引起.

    圖7 一階精度下兩相三波點問題密度和體積分數(shù)分布Fig.7 Density and volume fraction distributions for two-phase triple point problems using 1st order scheme

    圖8 兩相三波問題五階精度下采用HLL和HLLC所得密度和體積分數(shù)分布的對比((a)和(b)采用HLL格式計算,(c)和(d)采用HLLC格式計算)Fig.8 Density and volume fraction distributions using 5th WENO+HLL/HLLC methods(HLL is used in(a)and(b),HLLC is used in(c)and(d))

    圖8 兩相三波問題五階精度下采用HLL和HLLC所得密度和體積分數(shù)分布的對比((a)和(b)采用HLL格式計算,(c)和(d)采用HLLC格式計算)(續(xù))Fig.8 Density and volume fraction distributions using 5th WENO+HLL/HLLC methods(HLL is used in(a)and(b),HLLC is used in(c)and(d))(continued)

    圖9 五階精度下兩相三波問題DSINC和THINC所得密度和體積分數(shù)分布的對比((a)和(b)采用HLL格式計算,(c)和(d)采用HLLC格式計算)Fig.9 Density and volume fraction distributions using 5th WENO and DSINC/THINC methods(HLL is used in(a)and(b),HLLC is used in(c)and(d))

    3.3 激波與氣泡相互作用問題

    激波界面相互作用問題是兩相流數(shù)值模擬中另一個典型的問題[22],其中激波氣泡相互作用問題經(jīng)常用作數(shù)值驗證和實驗對比分析[23-24].在本文的數(shù)值模擬中,氣泡內(nèi)的氣體為R22,周圍為空氣.在初始時刻,馬赫數(shù)為1.22的左行激波位于x=275mm處,而氣泡的半徑R=25mm,位于(x,y)=(225,44.5)mm之處.R22氣體采用理想氣體模型:γ1=1.249,ρ01=3.863kg/m3;空氣的相應參數(shù)為:γ2=1.4,ρ02=1.225kg/m3.在氣泡內(nèi),流場初始變量設(shè)為

    激波后區(qū)域

    其他區(qū)域

    在此問題的模擬中,計算區(qū)域的大小為(x,y)∈[0,445]×[0,89]mm,計算網(wǎng)格數(shù)為1780×178.數(shù)值模擬結(jié)果如圖10所示.計算中除所采用的界面方法不同外,其他條件均相同.從圖中可以看出,對于相同的時間(t=0.3ms)的計算結(jié)果,DSINC和THINC方法所得結(jié)果略有差異.從計算結(jié)果來看,較長時間的數(shù)值模擬之后,密度所表示的界面依然可以清晰分辨出來,兩種方法所計算的體積分數(shù)等值線較為相似,界面較為清晰.由此可知,對于本文所構(gòu)造的雙正弦界面重構(gòu)方法來說,可以獲得與THINC相當?shù)挠嬎阈Ч?

    圖10 采用DSINC和THINC界面方法所得密度、體積分數(shù)和壓力等值線圖對比Fig.10 Density,volume fraction and pressure distributions using DSINC and THINC interface methods

    圖10 采用DSINC和THINC界面方法所得密度、體積分數(shù)和壓力等值線圖對比(續(xù))Fig.10 Density,volume fraction and pressure distributions using DSINC and THINC interface methods(continued)

    4 結(jié)論

    本文的研究發(fā)展了一種新的雙界面函數(shù)重構(gòu)方法,并基于此思想采用雙正弦函數(shù)構(gòu)造了雙正弦界面重構(gòu)方法DSINC.隨后在數(shù)值方法的驗證中,采用作者所在課題組發(fā)展的CFD2程序,分別對比了一階、五階精度數(shù)值模擬中無界面重構(gòu)、DSINC界面重構(gòu)和THINC界面重構(gòu)計算結(jié)果的不同.從圓柱對流問題、兩相三波點問題和激波與界面相互作用問題的對比發(fā)現(xiàn),DSINC與THINC出現(xiàn)了基本一致的界面重構(gòu)效果.另外,雙界面函數(shù)重構(gòu)時,由于函數(shù)的易于積分的特性,左、右特征變量并沒有出現(xiàn)復雜的表達式,這一特點使其更易于向多維擴展.

    1 趙寧,王東紅.多介質(zhì)流體界面問題的數(shù)值模擬.北京:科學出版社,2016(Zhao Ning,Wang Donghong.Numerical Simulation of Multimaterial Fluid Interface Problems.Beijing:Science Press,2016(in Chinese))

    2 Jesus WC,Roma AM,Pivello MR,et al.A 3D front-tracking approach for simulation of a two-phase fl uid with insoluble surfactant.Journal of Computational Physics,2015,281(c):403-420

    3 Terashima H,Tryggvason G.A front-tracking/ghost- fl uid method for fl uid interfaces in compressible fl ows.Journal of Computational Physics,2009,228(11):4012-4037

    4 Tryggvason G,Bunner B,Esmaeeli A,et al.A front-tracking method for the computations of multiphase fl ow.Journal of Computational Physics,2001,169(2):708-759

    5 Lu H,Zhao N,Wang D.A front tracking method for the simulation of compressible multimedium fl ows.Communications in Computational Physics,2016,19(1):124-142

    6 ScardovelliR,ZaleskiS.Directnumericalsimulationoffree-surface and interfacial fl ow.Annu Rev Fluid Mech,1999,31(1):567-603

    7 Torres DJ,Brackbill JU.The point-set mothod:front-tracking without connectivity.Journal of Computational Physics,2000,165(2):620-644

    8 Fedkiw RP,Sapiro G,Shu C-W.Shock capturing,Level sets and PDE based methods in computer vision and image processing:A review of Osher’s constributions.Journal of Computational Physics,2003,185(2):309-341

    9 Sussman M,Smereka P,Osher S.A level set approach for computing solutions to incompressible two-phase fl ow.Journal of Computational Physics,1994,114(1):146-159

    10 Fleckenstein S,Bothe D.A volume-of- fl uid-based numerical method for multi-component mass transfer with local volume changes.Journal of Computational Physics,2015,301(c):35-58

    11 Shyue K-M,Xiao F.An Eulerian interface sharpening algorithm for compressible two-phase fl ow:The algebraic THINC approach.Journal of Computational Physics,2014,268(2):326-354

    12 Susman M,Puckett EG.A coupled level set and volume-of- fl uid method for computing 3d and axisymmetric incompressible twophase fl ows.Journal of Computational Physics,2000,162(2):301-337

    13 So KK,Hu XY,Adams NA.Anti-di ff usion method for interface steepening in two-phase incompressible fl ow.Journal of Computational Physics,2011,230(13):5155-5177

    14 So KK,Hu XY,Adams NA.Anti-di ff usion interface sharpening technique for two-phase compressible fl ow simulations.Journal of Computational Physics,2012,231(11):4304-4323

    15 Youngs DL.An interface tracking method for a 3D Eulerian hydrodynamics code.Technical Report 44/92/35,AWRE,1984

    16 Xiao F,Honma Y,Kono T.A simple algebraic interface capturing scheme using hyperbolic tangent function.International Journal for Numerical Methods in Fluids,2005,48(9):1023-1040

    17 Xiao F,Li S,Chen C.Revisit to the THINC scheme:A simple algebraic VOF algorithm.Journal of Computational Physics,2011,230(19):7086-7092

    18 Li S,Sugiyama K,Takeuchi S,et al.An interface capturing method with a continuous function:The THINC method with multidimensional reconstruction.Journal of Computational Physics,2012,231(5):2328-2358

    19 Cassidy DA,Edwards JR,Tian M.An investigation of interfacesharpening schemes for multi-phase mixture fl ows.Journal of Computational Physics,2009,228(16):5628-5649

    20 Toro EF.Riemann Solvers and Numerical Methods for Fluid Dynamics(Third Edition).Berlin,Heidelberg:Springer-Verlag,2009

    21 Johnsen E,Colonius T.Implementation of WENO schemes in compressible multicomponent fl ow problems.Journal of Computational Physics,2006,219(2):715-732

    22 蔣華,董剛,陳霄.小擾動界面形態(tài)對RM不穩(wěn)定性影響的數(shù)值分析.力學學報,2014,46(4):544-552(Jiang Hua,Dong Gang,Chen Xiao.Numerical study on the e ff ects of small-amplitude initial perturbations on RM instability.Chinese Journal of Theoretical and Applied Mechanics,2014,46(4):544-552(in Chinese))

    23 王革,關(guān)奔.激波作用下R22氣泡射流現(xiàn)象研究.力學學報,2013,45(5):707-715(Wang Ge,Guan Ben.A study on jet phenomenon of R22 gas cylinder under the impact of shock wave.Chinese Journal of Theoretical and Applied Mechanics,2013,45(5):707-715(in Chinese))

    24 王顯圣,司廷,羅喜勝等.反射激波沖擊重氣柱的RM不穩(wěn)定性數(shù)值研究.力學學報,2012,44(4):666-674(Wang Xiansheng,Si Ting,Luo Xisheng,et al.Numerical study on the RM instability of a heavy-gas cylinder interacted with reshock.Chinese Journal of Theoretical and Applied Mechanics,2012,44(4):664-672(in Chinese))

    A NEW INTERFACE CAPTURING METHOD BASED ON DOUBLE INTERFACE FUNCTIONS1)

    Li Kang*,?Liu Na*,?He Zhiwei*Luo Longshan*,?Tian Baolin*,2)
    *(Institute of Applied Physics and Computational Mathematics,Beijing100088,China)
    ?(CAEP Software Center for High Performance Numerical Simulation,Beijing100088,China)

    Wedescribeanoveldouble-interface-function(DIF)reconstructionmethodforefficientnumericalresolutionof a compressible two-phase fl ow.Based on the new method,double sine interface capturing scheme(DSINC)is obtained.Five-equation model is solved to analyze the e ff ect of di ff erent interface functions such as DIF and Single Interface function(SIF)on the interfaces captured numerically.Near the interfaces,the algorithm uses the DIF or SIF as a basis for the reconstruction of a sub-grid discontinuity of volume fractions.In regions away from the interfaces,WENO is used to reconstruct the convective term,and time integration of the algorithm is done by employing the TVD Runge-Kutta method.Comparing with tangent of hyperbola for interface capturing(THINC)using SIF method,the left and right states reconstructed by DSINC is simpler and we need not solve a transcendental equation.Numerical results are shown with the Mie-Grneisen equation of state(EOS)for sample problems such as discontinuous advection,two-phase triple problem and shock-bubble interaction problem with THINC and DSINC.It can be found that DSINC is able to get as efficient resolution interface as THINC and shows to be more stable in the simulation.

    multi-phase fl ow,interface capturing,THINC,double interface functions reconstruction,DSINC

    O359+.1

    A doi:10.6052/0459-1879-17-210

    2017–09–04 收稿,2017–09–12 錄用,2017–09–13 網(wǎng)絡(luò)版發(fā)表.

    1)國家自然科學基金資助項目(11472059,U1630247,1730111).

    2)田保林,研究員,主要研究方向:計算流體力學數(shù)值方法.E-mail:tian_baolin@iapcm.ac.cn

    李康,劉娜,何志偉,駱龍山,田保林.一種基于雙界面函數(shù)的界面捕捉方法.力學學報,2017,49(6):1290-1300

    Li Kang,Liu Na,He Zhiwei,Luo Longshan,Tian Baolin.A new interface capturing method based on double interface functions.Chinese Journal of Theoretical and Applied Mechanics,2017,49(6):1290-1300

    猜你喜歡
    正弦重構(gòu)數(shù)值
    用固定數(shù)值計算
    例說正弦定理的七大應用
    正弦、余弦定理的應用
    長城敘事的重構(gòu)
    攝影世界(2022年1期)2022-01-21 10:50:14
    數(shù)值大小比較“招招鮮”
    北方大陸 重構(gòu)未來
    “美”在二倍角正弦公式中的應用
    北京的重構(gòu)與再造
    商周刊(2017年6期)2017-08-22 03:42:36
    論中止行為及其對中止犯的重構(gòu)
    基于Fluent的GTAW數(shù)值模擬
    焊接(2016年2期)2016-02-27 13:01:02
    亚洲人成网站高清观看| 国产精品女同一区二区软件| 免费大片18禁| 国产精品三级大全| 视频中文字幕在线观看| 国产淫片久久久久久久久| 午夜精品国产一区二区电影 | 好男人视频免费观看在线| 午夜爱爱视频在线播放| 嘟嘟电影网在线观看| 国产色爽女视频免费观看| 免费av观看视频| 日本-黄色视频高清免费观看| 久久国产乱子免费精品| 丰满乱子伦码专区| 亚洲精品影视一区二区三区av| 久久久久性生活片| 久久久久久久亚洲中文字幕| 熟女电影av网| 在线观看一区二区三区| 国产精品三级大全| 亚洲av免费高清在线观看| 中文字幕制服av| 中文字幕av在线有码专区| 高清午夜精品一区二区三区| 国产精品久久久久久精品电影| 久久久精品欧美日韩精品| 午夜福利在线观看吧| 国产亚洲一区二区精品| 精品一区二区三卡| 一个人看的www免费观看视频| av在线观看视频网站免费| 日韩大片免费观看网站| 国产黄片视频在线免费观看| 三级国产精品片| 精品一区二区三卡| 免费观看在线日韩| 91av网一区二区| 丰满乱子伦码专区| 国产精品久久久久久精品电影小说 | 色综合亚洲欧美另类图片| 欧美一区二区亚洲| 好男人在线观看高清免费视频| 欧美成人午夜免费资源| 伊人久久精品亚洲午夜| 国产黄色免费在线视频| 777米奇影视久久| 国产成人a∨麻豆精品| 色视频www国产| 可以在线观看毛片的网站| 插逼视频在线观看| 草草在线视频免费看| 美女xxoo啪啪120秒动态图| 好男人在线观看高清免费视频| 国产在视频线在精品| 久久鲁丝午夜福利片| 久久久午夜欧美精品| 国产男人的电影天堂91| 亚洲精品国产av成人精品| 亚洲精品乱码久久久v下载方式| 伦精品一区二区三区| 男的添女的下面高潮视频| 国产 亚洲一区二区三区 | 免费观看无遮挡的男女| 亚洲精品中文字幕在线视频 | 亚洲av福利一区| 大片免费播放器 马上看| 男人狂女人下面高潮的视频| 免费看光身美女| 日产精品乱码卡一卡2卡三| 高清av免费在线| 春色校园在线视频观看| 尤物成人国产欧美一区二区三区| 大陆偷拍与自拍| 日韩伦理黄色片| 国产精品女同一区二区软件| 国内精品宾馆在线| 高清欧美精品videossex| 欧美最新免费一区二区三区| 欧美97在线视频| 欧美高清成人免费视频www| 啦啦啦韩国在线观看视频| 久久99热6这里只有精品| 国产中年淑女户外野战色| 国产麻豆成人av免费视频| 亚洲国产最新在线播放| 2021少妇久久久久久久久久久| 免费观看的影片在线观看| 午夜福利视频1000在线观看| 午夜免费激情av| 国产精品熟女久久久久浪| 乱系列少妇在线播放| 精品人妻熟女av久视频| 国产精品久久久久久精品电影小说 | 一级爰片在线观看| av天堂中文字幕网| 精品国产三级普通话版| 男女视频在线观看网站免费| videos熟女内射| 黄色日韩在线| 亚洲一级一片aⅴ在线观看| 国产精品久久久久久av不卡| 亚洲美女搞黄在线观看| 久久久久久久久久人人人人人人| 国产精品国产三级国产av玫瑰| 91久久精品国产一区二区成人| 免费看日本二区| 国产伦在线观看视频一区| freevideosex欧美| 丰满乱子伦码专区| 建设人人有责人人尽责人人享有的 | 综合色丁香网| 婷婷色麻豆天堂久久| 精品人妻视频免费看| 亚洲最大成人av| 成年人午夜在线观看视频 | 久久久久精品久久久久真实原创| 国产精品久久久久久av不卡| 午夜精品一区二区三区免费看| 精品熟女少妇av免费看| 亚洲精品色激情综合| 美女cb高潮喷水在线观看| 国产亚洲av片在线观看秒播厂 | 日韩成人av中文字幕在线观看| 又爽又黄无遮挡网站| 免费看光身美女| 少妇被粗大猛烈的视频| av卡一久久| 精品99又大又爽又粗少妇毛片| 黄色配什么色好看| 成人国产麻豆网| 亚洲国产精品成人综合色| 成人亚洲精品av一区二区| 黄色配什么色好看| 亚洲伊人久久精品综合| 一本久久精品| 国产欧美日韩精品一区二区| 国产黄频视频在线观看| 欧美激情在线99| 午夜福利成人在线免费观看| 九九在线视频观看精品| 久久久国产一区二区| 国产美女午夜福利| 国产成人一区二区在线| 中文欧美无线码| 观看美女的网站| 蜜桃亚洲精品一区二区三区| 人妻制服诱惑在线中文字幕| 亚洲av.av天堂| 毛片女人毛片| 国产午夜精品一二区理论片| 亚洲自拍偷在线| 国产黄色小视频在线观看| 国产精品久久久久久久久免| 精品不卡国产一区二区三区| 国产精品国产三级国产av玫瑰| 成人无遮挡网站| 在线观看人妻少妇| 两个人的视频大全免费| 美女国产视频在线观看| 亚洲国产欧美人成| 亚洲18禁久久av| 一级毛片久久久久久久久女| 一区二区三区免费毛片| 国产在视频线精品| 麻豆精品久久久久久蜜桃| 久久精品国产亚洲网站| 久久久a久久爽久久v久久| 国产激情偷乱视频一区二区| 一个人看视频在线观看www免费| 日本黄大片高清| 亚洲精品久久久久久婷婷小说| 男女下面进入的视频免费午夜| 亚洲一级一片aⅴ在线观看| 国产精品精品国产色婷婷| 插阴视频在线观看视频| 高清毛片免费看| 大陆偷拍与自拍| 亚洲18禁久久av| 国产麻豆成人av免费视频| 国产伦一二天堂av在线观看| 国产精品国产三级国产av玫瑰| 久久久久精品性色| 中文资源天堂在线| 国产不卡一卡二| 欧美xxxx黑人xx丫x性爽| 国产大屁股一区二区在线视频| 成人毛片60女人毛片免费| 日韩强制内射视频| 特大巨黑吊av在线直播| 夫妻性生交免费视频一级片| 国国产精品蜜臀av免费| 嘟嘟电影网在线观看| av天堂中文字幕网| 国国产精品蜜臀av免费| 女的被弄到高潮叫床怎么办| 国产午夜福利久久久久久| 国产亚洲91精品色在线| 99热全是精品| 看十八女毛片水多多多| 亚洲18禁久久av| 一级黄片播放器| 伊人久久精品亚洲午夜| 亚洲精品色激情综合| 午夜激情欧美在线| 高清视频免费观看一区二区 | 国产一区有黄有色的免费视频 | 亚洲欧美清纯卡通| 日韩电影二区| eeuss影院久久| 午夜老司机福利剧场| 又黄又爽又刺激的免费视频.| 91午夜精品亚洲一区二区三区| 99九九线精品视频在线观看视频| 99热这里只有精品一区| 久久久欧美国产精品| 亚洲av在线观看美女高潮| 爱豆传媒免费全集在线观看| 高清日韩中文字幕在线| 日韩强制内射视频| 亚洲av一区综合| 简卡轻食公司| 极品少妇高潮喷水抽搐| 国产午夜精品论理片| 天堂中文最新版在线下载 | 日韩国内少妇激情av| 国产午夜精品论理片| av网站免费在线观看视频 | 天天躁日日操中文字幕| 青春草视频在线免费观看| 日韩av在线免费看完整版不卡| 国精品久久久久久国模美| 国产日韩欧美在线精品| 神马国产精品三级电影在线观看| videossex国产| 成年女人在线观看亚洲视频 | 尾随美女入室| 亚洲欧美成人综合另类久久久| 亚洲va在线va天堂va国产| 成人一区二区视频在线观看| 久久精品国产自在天天线| 99久久精品国产国产毛片| 人人妻人人看人人澡| 日韩中字成人| 国产精品一区www在线观看| 久久久久久久久久成人| www.色视频.com| 能在线免费观看的黄片| 国产精品一区二区三区四区免费观看| 男女边吃奶边做爰视频| 亚洲久久久久久中文字幕| 日本黄大片高清| 日本av手机在线免费观看| 亚洲在线自拍视频| 日韩欧美精品v在线| 最后的刺客免费高清国语| 免费看不卡的av| 日韩精品有码人妻一区| 在线免费十八禁| 免费看光身美女| 亚洲av福利一区| 极品教师在线视频| 永久网站在线| 91久久精品国产一区二区成人| 久久精品熟女亚洲av麻豆精品 | 久久久久性生活片| 亚洲第一区二区三区不卡| 国产精品精品国产色婷婷| 午夜久久久久精精品| 男人舔奶头视频| 男女视频在线观看网站免费| 国产伦在线观看视频一区| 亚洲一级一片aⅴ在线观看| 99热这里只有是精品在线观看| 精品国产一区二区三区久久久樱花 | 看黄色毛片网站| 久久久久精品久久久久真实原创| 建设人人有责人人尽责人人享有的 | 熟女人妻精品中文字幕| 三级经典国产精品| 国产亚洲最大av| 色尼玛亚洲综合影院| 久久久午夜欧美精品| 午夜老司机福利剧场| av在线亚洲专区| av免费在线看不卡| 婷婷色av中文字幕| 纵有疾风起免费观看全集完整版 | 亚洲人成网站在线观看播放| 亚洲18禁久久av| 欧美日韩亚洲高清精品| 日韩一本色道免费dvd| 男女视频在线观看网站免费| 熟妇人妻久久中文字幕3abv| 麻豆国产97在线/欧美| 久久久久久久国产电影| 亚洲精品亚洲一区二区| 亚洲av免费高清在线观看| 日韩,欧美,国产一区二区三区| xxx大片免费视频| 午夜精品国产一区二区电影 | 搡女人真爽免费视频火全软件| 水蜜桃什么品种好| 最近中文字幕2019免费版| 综合色av麻豆| 一级毛片aaaaaa免费看小| 男人舔女人下体高潮全视频| 国产69精品久久久久777片| 亚洲电影在线观看av| 国产综合精华液| 日本一本二区三区精品| 婷婷六月久久综合丁香| 久久人人爽人人片av| 天堂中文最新版在线下载 | av女优亚洲男人天堂| 国产毛片a区久久久久| 国产精品三级大全| 日本一本二区三区精品| 日韩制服骚丝袜av| 国产v大片淫在线免费观看| 天天躁夜夜躁狠狠久久av| 能在线免费看毛片的网站| 国产探花极品一区二区| a级毛片免费高清观看在线播放| 日日啪夜夜爽| 3wmmmm亚洲av在线观看| 久久久久性生活片| ponron亚洲| 欧美xxxx黑人xx丫x性爽| 熟妇人妻久久中文字幕3abv| 我要看日韩黄色一级片| 久久这里有精品视频免费| 九九久久精品国产亚洲av麻豆| 国产精品国产三级国产专区5o| 成年女人看的毛片在线观看| 大片免费播放器 马上看| 国产伦一二天堂av在线观看| 激情 狠狠 欧美| 国产色婷婷99| 我的女老师完整版在线观看| 精品久久久久久成人av| 91aial.com中文字幕在线观看| 国产精品一区www在线观看| 国产精品久久久久久久电影| 国产成人a∨麻豆精品| 麻豆成人av视频| 久久久国产一区二区| 天堂影院成人在线观看| 最新中文字幕久久久久| 极品教师在线视频| 国产精品国产三级专区第一集| 国产老妇伦熟女老妇高清| 看免费成人av毛片| 我要看日韩黄色一级片| 中文字幕av在线有码专区| 亚洲av福利一区| av黄色大香蕉| 国产三级在线视频| 欧美成人a在线观看| 夜夜爽夜夜爽视频| 99视频精品全部免费 在线| 草草在线视频免费看| 毛片一级片免费看久久久久| 天堂俺去俺来也www色官网 | 亚洲熟女精品中文字幕| 丰满乱子伦码专区| 麻豆成人午夜福利视频| 久久久久久久久中文| 国产免费福利视频在线观看| 国产亚洲午夜精品一区二区久久 | 国产v大片淫在线免费观看| 色综合色国产| 日韩伦理黄色片| 最后的刺客免费高清国语| 少妇高潮的动态图| 亚洲aⅴ乱码一区二区在线播放| 国产伦一二天堂av在线观看| 深爱激情五月婷婷| 91狼人影院| 久久精品综合一区二区三区| 成年女人在线观看亚洲视频 | 男插女下体视频免费在线播放| 久久久久免费精品人妻一区二区| 一级毛片电影观看| 亚洲丝袜综合中文字幕| 中文精品一卡2卡3卡4更新| 精品亚洲乱码少妇综合久久| 岛国毛片在线播放| 久久久久久久午夜电影| 成人欧美大片| 一边亲一边摸免费视频| 精品一区二区三卡| 国产av不卡久久| 久久久a久久爽久久v久久| 亚洲伊人久久精品综合| 午夜久久久久精精品| 亚洲国产av新网站| 又爽又黄a免费视频| 日本wwww免费看| 如何舔出高潮| 哪个播放器可以免费观看大片| 久久久久性生活片| 人人妻人人澡欧美一区二区| 国产69精品久久久久777片| 亚洲精品视频女| 亚洲国产精品成人久久小说| 国产精品综合久久久久久久免费| 成年人午夜在线观看视频 | 亚洲高清免费不卡视频| 国产白丝娇喘喷水9色精品| 国产精品一及| 亚洲av不卡在线观看| 黄片无遮挡物在线观看| 久久99热6这里只有精品| 亚洲精品久久久久久婷婷小说| 99热6这里只有精品| av网站免费在线观看视频 | 国产有黄有色有爽视频| 一本久久精品| 日本黄大片高清| 亚洲四区av| 一级毛片黄色毛片免费观看视频| 国产亚洲精品久久久com| 美女脱内裤让男人舔精品视频| 成人无遮挡网站| 人体艺术视频欧美日本| 丰满少妇做爰视频| 国产av码专区亚洲av| 亚洲欧美清纯卡通| 免费观看av网站的网址| 国产亚洲av嫩草精品影院| 亚洲第一区二区三区不卡| 2021少妇久久久久久久久久久| 国产成人精品久久久久久| 亚洲精品国产av成人精品| 日本免费在线观看一区| 国产亚洲av嫩草精品影院| 国产成人精品久久久久久| 老司机影院成人| 日韩欧美 国产精品| 人妻少妇偷人精品九色| 久久精品人妻少妇| 在线免费观看不下载黄p国产| 亚洲精品aⅴ在线观看| 日韩亚洲欧美综合| 中文天堂在线官网| 国产黄色视频一区二区在线观看| 男女啪啪激烈高潮av片| 亚洲熟妇中文字幕五十中出| 午夜爱爱视频在线播放| 免费av不卡在线播放| 嫩草影院入口| 人体艺术视频欧美日本| 一级毛片我不卡| 精品国产一区二区三区久久久樱花 | 精品人妻偷拍中文字幕| 97精品久久久久久久久久精品| 男女边摸边吃奶| 99热这里只有是精品在线观看| 亚洲自拍偷在线| 特大巨黑吊av在线直播| 九九久久精品国产亚洲av麻豆| 99久久人妻综合| 丰满人妻一区二区三区视频av| av在线老鸭窝| 又粗又硬又长又爽又黄的视频| 大又大粗又爽又黄少妇毛片口| 久久久久久久午夜电影| 国产成人一区二区在线| 亚洲国产精品成人综合色| 在线a可以看的网站| 免费观看的影片在线观看| 午夜福利网站1000一区二区三区| 少妇熟女欧美另类| 日本与韩国留学比较| 国产精品伦人一区二区| 国产一区二区三区av在线| 99久国产av精品| 少妇丰满av| 中文字幕人妻熟人妻熟丝袜美| 美女黄网站色视频| 免费观看精品视频网站| 丰满乱子伦码专区| 精品少妇黑人巨大在线播放| 大片免费播放器 马上看| 肉色欧美久久久久久久蜜桃 | 成年女人在线观看亚洲视频 | 欧美3d第一页| av国产免费在线观看| 亚洲av在线观看美女高潮| 97超视频在线观看视频| 波多野结衣巨乳人妻| 免费大片18禁| .国产精品久久| 亚洲av免费高清在线观看| 日韩一本色道免费dvd| 欧美3d第一页| 亚洲国产av新网站| 天堂网av新在线| av一本久久久久| 十八禁网站网址无遮挡 | 91久久精品电影网| 韩国高清视频一区二区三区| 日韩中字成人| av在线老鸭窝| 日韩强制内射视频| 黄色一级大片看看| 国产亚洲精品久久久com| 中文字幕av成人在线电影| 婷婷色av中文字幕| 美女国产视频在线观看| 一级毛片aaaaaa免费看小| 成人漫画全彩无遮挡| 欧美不卡视频在线免费观看| 一个人观看的视频www高清免费观看| 欧美日韩一区二区视频在线观看视频在线 | 99九九线精品视频在线观看视频| 少妇熟女aⅴ在线视频| 两个人视频免费观看高清| 免费无遮挡裸体视频| 日韩欧美精品免费久久| 天天躁夜夜躁狠狠久久av| 超碰97精品在线观看| 免费播放大片免费观看视频在线观看| 久久久久久久亚洲中文字幕| 久久久精品欧美日韩精品| 国产精品不卡视频一区二区| 久久久久久久久久黄片| 99热这里只有精品一区| 91aial.com中文字幕在线观看| 午夜爱爱视频在线播放| 一个人看的www免费观看视频| 欧美日韩视频高清一区二区三区二| 日韩一区二区视频免费看| 在线免费观看的www视频| 在现免费观看毛片| 亚洲国产高清在线一区二区三| 欧美激情在线99| 午夜免费男女啪啪视频观看| 亚洲精品色激情综合| 中文字幕制服av| 日韩欧美精品v在线| 国产黄色视频一区二区在线观看| 天堂中文最新版在线下载 | 97在线视频观看| 成人美女网站在线观看视频| eeuss影院久久| 国产一级毛片在线| 午夜日本视频在线| 简卡轻食公司| 婷婷色综合大香蕉| 丝袜喷水一区| 综合色av麻豆| 一个人看的www免费观看视频| 成年人午夜在线观看视频 | 最近的中文字幕免费完整| 日韩欧美精品v在线| 亚洲天堂国产精品一区在线| 男女国产视频网站| 成人亚洲精品一区在线观看 | 久久精品综合一区二区三区| 九草在线视频观看| 精品欧美国产一区二区三| 别揉我奶头 嗯啊视频| 国产精品三级大全| 国产毛片a区久久久久| 好男人在线观看高清免费视频| 亚洲成人精品中文字幕电影| 国产综合懂色| 99九九线精品视频在线观看视频| 亚洲成人av在线免费| 免费黄色在线免费观看| 看非洲黑人一级黄片| 成人综合一区亚洲| 久久综合国产亚洲精品| 你懂的网址亚洲精品在线观看| 日韩一区二区三区影片| 久久久久久伊人网av| 人人妻人人看人人澡| 嫩草影院入口| 日韩av在线免费看完整版不卡| 欧美高清性xxxxhd video| 麻豆国产97在线/欧美| 大香蕉久久网| 亚洲欧美中文字幕日韩二区| 亚洲av国产av综合av卡| 亚洲精品一区蜜桃| 国产亚洲最大av| 老司机影院毛片| 国产中年淑女户外野战色| 九九久久精品国产亚洲av麻豆| 午夜视频国产福利| 性插视频无遮挡在线免费观看| 天天躁夜夜躁狠狠久久av| 欧美高清性xxxxhd video| 国国产精品蜜臀av免费| 一个人看的www免费观看视频| 亚洲欧美一区二区三区国产| 久久这里只有精品中国| 少妇的逼水好多| videossex国产| 一区二区三区高清视频在线| 国产69精品久久久久777片| 亚洲av免费在线观看| 亚洲在久久综合| 在线 av 中文字幕| 亚洲熟女精品中文字幕| 黄色欧美视频在线观看| 日韩成人av中文字幕在线观看| 日本黄大片高清| 国产精品蜜桃在线观看| 精品久久久久久久久av| 国产黄片美女视频| 非洲黑人性xxxx精品又粗又长| 国产精品久久久久久精品电影|