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

    一種滑移區(qū)氣體流動的格子Boltzmann曲邊界處理新格式?

    2017-08-09 07:34:22顧娟黃榮宗劉振宇吳慧英
    物理學(xué)報 2017年11期
    關(guān)鍵詞:偏移量無量邊界條件

    顧娟 黃榮宗 劉振宇 吳慧英

    (上海交通大學(xué)機(jī)械與動力工程學(xué)院,上海 200240)

    一種滑移區(qū)氣體流動的格子Boltzmann曲邊界處理新格式?

    顧娟 黃榮宗 劉振宇 吳慧英?

    (上海交通大學(xué)機(jī)械與動力工程學(xué)院,上海 200240)

    (2016年12月18日收到;2017年4月5日收到修改稿)

    針對滑移區(qū)復(fù)雜氣-固邊界存在速度滑移現(xiàn)象,提出了一種基于格子Boltzmann方法的非平衡態(tài)外推與有限差分相結(jié)合的曲邊界處理新格式.該格式具有可考慮實際物理邊界與網(wǎng)格線偏移量的優(yōu)勢,較傳統(tǒng)half-way DBB(di ff usive bounce-back)格式更能準(zhǔn)確反映實際邊界情況,同時還可獲取壁面處氣體宏觀量及其法向梯度等信息.采用本文所提曲邊界處理格式模擬分析了滑移區(qū)氣體平直/傾斜微通道Poiseuille流、微圓柱繞流和同心微圓柱面旋轉(zhuǎn)Couette流問題.研究結(jié)果表明,采用曲邊界處理新格式所得結(jié)果與理論值以及文獻(xiàn)結(jié)果符合良好,適用于滑移區(qū)氣體流動的復(fù)雜邊界處理,且比half-way DBB格式具有更高的精度,較修正DBB格式具有更好的適應(yīng)性.

    曲邊界,格子Boltzmann方法,滑移區(qū),微氣體流動

    1 引言

    隨著微機(jī)電系統(tǒng)(micro-electro-mechanical system,MEMS)技術(shù)在生物醫(yī)療、信息通訊及國防軍工等諸多領(lǐng)域的逐步應(yīng)用,微尺度條件下新的物理現(xiàn)象引起了越來越多學(xué)者的重視[1,2].其中對于微尺度氣體流動[3,4]問題,因氣體分子平均自由程λ與特征尺寸H兩者數(shù)量級相當(dāng),氣體稀薄效應(yīng)已不可忽略.氣體稀薄程度可用Knudsen數(shù)表示,Knudsen數(shù)(Kn)的定義為Kn=λ/H.滑移區(qū)(0.001 6 Kn 6 0.1)氣體在氣-固邊界處將出現(xiàn)明顯的速度滑移現(xiàn)象,基于連續(xù)介質(zhì)假定的Navier-Stokes方程需結(jié)合滑移邊界條件方可描述該問題.有別于傳統(tǒng)計算流體動力學(xué)方法,格子Boltzmann(LB)方法不受連續(xù)介質(zhì)假定限制,具有微觀粒子特性和介觀物理背景,成為研究微尺度氣體流動問題的理想數(shù)值方法[5?7].

    目前,針對常規(guī)尺度曲邊界處理的研究已有較多文獻(xiàn)報道[8?22],其中非平衡外推格式[8]已被證明可有效處理復(fù)雜曲邊界問題[13,20],前人也對Ladd[10]提出的考慮邊界速度的half-way反彈格式進(jìn)行了修正以推廣到曲面邊界條件.Yin和Zhang[15]基于Ladd[10]的half-way反彈格式將虛擬節(jié)點(實際物理邊界附近流體格點與固體格點連線的中點)處速度取代原格式中的壁面速度,提高了模擬復(fù)雜邊界問題的精度.在此基礎(chǔ)上,史冬巖等[16]對離流體格點較近的虛擬節(jié)點速度進(jìn)行修正,進(jìn)一步提高了計算精度.盡管以上方法在常規(guī)尺度曲邊界處理中表現(xiàn)出較好的優(yōu)越性,但由于其不能反映因微尺度氣體-壁面相互作用引起的邊界速度滑移,無法適用于微尺度滑移區(qū)流動問題中復(fù)雜曲邊界的處理.Szalmas[23]采用貼體曲坐標(biāo)下SRB(specular-re fl ection-bounceback)格式有效地模擬了微尺度復(fù)雜氣-固邊界速度滑移問題,但貼體網(wǎng)格需要求解Jacobi變換矩陣,增加了編程復(fù)雜度.另一類處理曲邊界的方法是采用笛卡爾網(wǎng)格系統(tǒng),不需要進(jìn)行網(wǎng)格變換,方便處理各種復(fù)雜物理邊界,但問題在于網(wǎng)格邊界與物理邊界不重合,因此邊界條件的處理尤為重要[24].針對這個問題,Guo等[25]采用half-way DBB格式處理同心微圓柱面旋轉(zhuǎn)Couette流,但該方法僅當(dāng)虛擬節(jié)點與固體壁面重合時,才能體現(xiàn)其計算精度優(yōu)勢.Suga[26]將插值格式分別引入到BB(bounce-back)格式和DS(di ff usescattering)格式中組合為IDBB(interpolation di ff usive bounce-back)格式以處理微尺度流動問題,但其組合系數(shù)的選取仍參考傳統(tǒng)half-way DBB格式.Tao和Guo[24]在DBB格式中引入實際物理邊界與網(wǎng)格線偏移量這一參數(shù),并應(yīng)用其對模型松弛時間τq和參數(shù)r進(jìn)行修正,將適用于平直邊界的DBB格式推廣至滑移曲邊界情況.當(dāng)偏移量為半個格子時,該格式則還原為half-way DBB格式,但為保證參數(shù)τq和r取值的合理性,增加了滑移邊界條件中氣體-壁面相互作用參數(shù)取值的難度.Yang等[27]提出一種考慮速度滑移的浸入邊界格子Boltzmann方法,并用其模擬微氣體擠壓膜阻尼效應(yīng)引起的空腔內(nèi)剛性物體振動,但僅驗證了該方法在不考慮稀薄效應(yīng)情況下的準(zhǔn)確度.目前針對滑移區(qū)曲邊界的速度滑移問題,精度高、適用性廣的方法鮮有研究報道.

    本文提出了一種綜合Guo等[8]非平衡態(tài)外推及Le等[17]有限差分格式的曲邊界處理格式,該格式考慮了實際物理邊界與網(wǎng)格線的偏移量,較halfway DBB格式更能準(zhǔn)確地描述實際物理邊界的情況.由于不涉及參數(shù)取值合理性問題,本文所提格式較修正DBB格式更能適應(yīng)氣體-壁面相互作用參數(shù)的選取,同時還可獲得壁面處氣體宏觀量及其法向梯度等信息.基于本文所提格式對滑移區(qū)氣體平直/傾斜微通道Poiseuille流、微圓柱繞流和同心微圓柱面旋轉(zhuǎn)Couette流問題進(jìn)行LB數(shù)值模擬研究,將模擬結(jié)果與解析解及已有文獻(xiàn)結(jié)果進(jìn)行對比討論,以分析本文所提格式的準(zhǔn)確度與適用性.

    2 LB模型與曲邊界格式

    2.1 微流動的LB基本模型

    LB模型通過格點上微觀粒子分布函數(shù)的演化過程來描述流體的流動.本文采用單松弛(LBGK)模型[28],其演化方程為

    其中,ci為粒子在i方向的速度,δt為時間步長,fi(x,t)為t時刻x位置處以速度ci運動的粒子的分布函數(shù),τ為無量綱松弛時間.(x,t)為i方向粒子的平衡態(tài)分布函數(shù),

    Fi(x,t)為離散外力項[29]

    式中,ωi為權(quán)系數(shù),cs為格子聲速,G為外力.宏觀物理量密度ρ和速度u的計算式為

    對于二維問題,本文選用二維九速(D2Q9)模型,相應(yīng)的離散速度為

    其中,粒子遷移速率c=δx/δt(δx為格子步長),格子聲速權(quán)系數(shù)ω0=4/9,ω1=ω2=ω3=ω4=1/9,ω5=ω6=ω7=ω8=1/36.通過Chapman-Enskog分析可知,模型對應(yīng)的流體壓力P=,流體的運動黏度ν=(τ?0.5)δt.

    對于微尺度氣體流動問題,需在模型中引入Knudsen數(shù)的影響.由動理學(xué)理論可知,氣體分子平均自由程λ與動力黏度μ、壓力P、溫度T的關(guān)系式為結(jié)合Knudsen數(shù)的定義式Kn=λ/H及LB模型中運動黏度ν與無量綱松弛時間τ的關(guān)系式ν=c2s(τ?0.5)δt,可得模型中無量綱松弛時間與Knudsen數(shù)的關(guān)系式為

    2.2微尺度流動的曲邊界處理

    對于微尺度氣體流動問題,壁面處滑移邊界條件的處理至關(guān)重要.本文在Guo等[8]的非平衡態(tài)外推格式的基礎(chǔ)上,提出了一種復(fù)雜曲邊界上壁面滑移速度的邊界處理格式.如圖1所示的曲邊界,xf為流體格點,xs為固體內(nèi)部格點,xint為固體邊界格點.流體格點xf處i方向的未知分布函數(shù)fi(xf,t+δt)由固體格點xs處的分布函數(shù)fi(xs,t)演化得到,而fi(xs,t)采用Guo等[8]的非平衡態(tài)外推格式重構(gòu)得到.為此,首先將固體格點處的分布函數(shù)分解為平衡態(tài)和非平衡態(tài)兩部分,即

    式中us可由u(xint),u(xf)及u(xff)外推確定

    需要指出的是,對于微尺度氣體流動問題,壁面處氣體速度ug(即u(xint))不等于固體壁面速度uw,需進(jìn)一步考慮滑移速度uslip的影響.本文采用微尺度一階滑移曲邊界條件[30],

    其中,uτ為氣體在壁面處的切向速度矢量,r為半徑,σv為動量協(xié)調(diào)系數(shù).對于平直邊界(r→∞)也適用,此時(11)式括號內(nèi)第二項為0,與微尺度經(jīng)典的Maxwell-Smoluchowski一階滑移直邊界條件[1]表達(dá)式一致.對于完全漫反射固體壁面,有σv=1[3?7,23?27].針對復(fù)雜曲邊界,本文采用Le等[17]提出求解壁面法向溫度梯度的有限差分格式求解壁面處流體速度的法向梯度.如圖1所示,沿壁面法向n在流體內(nèi)部選取兩個虛擬點x′和x′′,滿足|x′′?x′|=|x′?xint|=δ.為保證緊鄰x′(x′′)的四個網(wǎng)格點均為流體格點,本文取δ=1.5δx.對于虛擬點x′處的流體宏觀速度,通過其緊鄰的流體網(wǎng)格點采用雙線性插值法確定,即有

    其中A1,A2,A3和A4分別為x′所在網(wǎng)格內(nèi)四塊子區(qū)域的面積(如圖1所示).采用相同方法可求得虛擬點x′′處的流體宏觀速度u′′.由虛擬點x′和x′′處的流體宏觀速度,結(jié)合(11)式,可求得壁面處氣體速度為

    其中第二項即為微尺度壁面處氣體滑移速度,R為固壁的曲率半徑.

    圖1 曲邊界示意圖Fig.1.Schematic illustration of the curved boundary.

    3 結(jié)果與討論

    為驗證本文提出的滑移區(qū)曲邊界滑移速度邊界處理格式的適用性,本節(jié)依次模擬研究平直/傾斜微通道Poiseuille流、微圓柱繞流和同心微圓柱面旋轉(zhuǎn)Couette流問題.

    3.1 平直微通道Poiseuille流

    微尺度Poiseuille流在MEMS中廣泛存在且有解析解,可用于驗證滑移邊界處理方法的準(zhǔn)確度.本文模擬了網(wǎng)格線與平直微通道壁面偏移量qδx(q=0.25,0.5,0.75,1.0)不同情況下的滑移區(qū)外力驅(qū)動Poiseuille流.微通道內(nèi)氣體受到沿流動方向恒定外力Gx=10?4的驅(qū)動,進(jìn)出口采用周期性邊界條件,固體邊界上的滑移速度邊界條件((11)式)采用本文提出的邊界處理格式處理.計算過程中,計算網(wǎng)格取為Nx×Ny=4×32.對于該微尺度Poiseuille流,沿截面方向速度分布的解析解為

    其中H為微通道的寬度.

    圖2(a)和圖2(b)分別給出了當(dāng)Kn=0.0194時,采用本文格式和half-way DBB格式計算得到的截面無量綱速度分布與解析解的對比情況.由圖2可知,本文格式得到的截面無量綱速度分布和解析解在不同q情況下均符合良好;而half-way DBB格式得到的結(jié)果僅在q=0.5時才與解析解一致,這是因為half-way DBB格式基于half-way反彈格式,當(dāng)偏移量為0.5δx時,half-way DBB格式中的虛擬邊界與實際物理邊界重合,固體邊界可由虛擬邊界準(zhǔn)確刻畫,故此時計算精度最高.但當(dāng)偏移量發(fā)生改變時,half-way DBB格式因默認(rèn)偏移量為0.5δx,故其模擬結(jié)果并未有相應(yīng)改變,仍對應(yīng)q=0.5條件下的計算結(jié)果.故對于偏移量不為0.5δx,half-way DBB格式已不能再準(zhǔn)確描述固體邊界的實際情況.

    為定量比較本文格式與half-way DBB格式的計算優(yōu)勢,定義了如下誤差量:1)最大相對誤差,

    2)相對誤差,

    其中,(i,j)為計算域上的網(wǎng)格點,u(i,j)為流體格點上的速度值,uexact(i,j)為流體格點上的速度的解析值.表1給出了不同q條件下(Kn=0.0194),本文格式與half-way DBB格式計算結(jié)果誤差對比情況.從表中可知,如前所述half-way DBB格式僅在q=0.5時誤差才達(dá)最小值,本文格式在其余不同偏移量條件下較half-way DBB格式計算精度都有明顯提高.

    圖2 (網(wǎng)刊彩色)不同q時平直微通道截面無量綱速度分布(Kn=0.0194)(a)本文格式模擬結(jié)果;(b)half-way DBB格式模擬結(jié)果Fig.2.(color online)Normalized velocity pro fi les of the aligned microchannel fl ow with di ff erent values of q(Kn=0.0194):(a)The present scheme;(b)the half-way DBB scheme.

    3.2 傾斜微通道Poiseuille流

    為進(jìn)一步驗證本文提出的格式能有效處理滑移區(qū)復(fù)雜的幾何邊界,對圖3所示的傾斜微通道內(nèi)外力驅(qū)動Poiseuille流進(jìn)行了模擬.微通道傾斜角度θ依次取為arctan(0.46),arctan(1.20)和arctan(3.72),計算網(wǎng)格取為Nx×NAB=100×100,對應(yīng)的y方向格點數(shù)為Ny=NAB+Nxtanθ.對于該傾斜微通道Poiseuille流,截面速度分布解析解仍如(14)式所示,其中y代表到微通道下壁面的垂直距離,H為微通道的寬度.

    圖3傾斜微通道示意圖Fig.3.Schematic illustration of the inclined microchannel.

    圖4 (a)和圖4(b)分別給出了當(dāng)Kn=0.0194時,采用本文格式和half-way DBB格式計算得到的截面無量綱速度分布與解析解對比情況.由圖4可知,在不同傾斜角度情況下,本文格式的模擬結(jié)果和解析解始終符合良好,而half-way DBB格式的模擬結(jié)果在微通道壁面處將顯著偏離理論值.這是因為halfway DBB格式實際物理邊界與網(wǎng)格線偏移量恒為0.5δx,使得傾斜的微通道壁面成為階梯狀.

    表2給出了不同傾斜角條件下(Kn=0.0194),本文格式與half-way DBB格式計算結(jié)果誤差對比情況.從表2可知,本文格式較half-way DBB格式在模擬不同傾斜角時最大相對誤差及相對誤差均明顯降低.

    圖4 (網(wǎng)刊彩色)不同傾斜角θ時,傾斜微通道截面無量綱速度分布(Kn=0.0194)(a)本文格式模擬結(jié)果;(b)half-way DBB格式模擬結(jié)果Fig.4.(color online)Normalized velocity pro fi les of the inclined microchannel fl ow with di ff erent inclined angles θ(Kn=0.0194):(a)The present scheme;(b)the half-way DBB scheme.

    表2 不同傾斜角條件下本文格式與half-way DBB格式計算誤差對比(Kn=0.0194)Table 2.Computational errors of the present scheme and the half-way DBB scheme with di ff erent inclined angles θ(Kn=0.0194).

    表1 不同q條件下本文格式與half-way DBB格式計算誤差對比(Kn=0.0194)Table 1.Computational errors of the present scheme and the half-way DBB scheme with di ff erent values of q(Kn=0.0194).

    3.3微圓柱氣體繞流

    圓柱繞流是流體力學(xué)中經(jīng)典的復(fù)雜外流問題,常被用于檢測曲邊界的處理能力.本節(jié)對微圓柱氣體繞流問題進(jìn)行模擬,并將模擬結(jié)果與文獻(xiàn)結(jié)果[31]進(jìn)行對比驗證.對于微圓柱繞流問題,Reynolds數(shù)的定義為

    其中,U∞為均勻來流速度,D為微圓柱直徑.在模擬過程中,計算區(qū)域選取為40D×20D(D劃分為30個格子),微圓柱放置在離進(jìn)口及上下邊界均為10D的位置.計算區(qū)域邊界條件設(shè)置如下:

    圖5 (網(wǎng)刊彩色)微圓柱繞流(Re=20,Kn=0.015)(a)微圓柱附近流線;(b)微圓柱表面無量綱滑移速度分布Fig.5.(color online)Flow past a microcylinder(Re=20,Kn=0.015):(a)Streamlines around the microcylinder;(b)normalized slip velocity pro fi les on the surface of the microcylinder.

    微圓柱表面則采用(11)式所示的一階滑移曲邊界條件.

    圖5(a)給出了Re=20時,考慮稀薄效應(yīng)(Kn=0.015)與不考慮稀薄效應(yīng)(Kn=0)情況下微圓柱附近流線分布的對比.從圖中可以看出,考慮稀薄效應(yīng)時回流區(qū)長度明顯縮短,這是由于氣體-壁面滑移速度使得壁面附近流體微團(tuán)動能增加,逆壓梯度減小,從而分離渦變小.圖5(b)給出了Kn=0.015,Re=20時,采用本文格式及half-way DBB格式模擬得到的微圓柱表面無量綱滑移速度分布同Beskok和Karniadakis[31]采用譜元方法計算得到結(jié)果的對比.由圖可知,本文格式得到的模擬結(jié)果與譜元方法計算結(jié)果更為接近,在分離點之前幾乎完全一致,在分離點之后存在一定偏差;而half-way DBB格式得到的模擬結(jié)果整體偏小,最大滑移速度較文獻(xiàn)[31]的結(jié)果存在較大偏差.

    3.4 同心微圓柱面Couette流

    上述算例均為靜止壁面上的滑移速度邊界條件,為驗證本文提出的格式亦能處理運動壁面上的滑移速度邊界條件,本文模擬了同心微圓柱面旋轉(zhuǎn)Couette流問題.如圖6所示,內(nèi)微圓柱的半徑為R1(R1=30δx),以恒定的角速度ω逆時針旋轉(zhuǎn);外微圓柱的半徑為R2,保持為靜止?fàn)顟B(tài);內(nèi)外微圓柱的半徑比R1:R2=3:5,其表面均為完全漫反射表面(σv1=σv2=1).

    對于同心微圓柱面旋轉(zhuǎn)Couette流,圓柱坐標(biāo)系下的Navier-Stokes方程可簡化為[23?26,30]

    其中r為半徑,uθ為切向速度.內(nèi)外圓柱壁面處的氣體速度邊界條件為

    由此可得速度分布解析解為

    其中系數(shù)A1,A2為

    圖6同心微圓柱面Couette流示意圖Fig.6.Schematic illustration of the concentric microcylinder Couette fl ow.

    圖7 給出了Kn=0.05和Kn=0.1時,采用本文提出的格式、half-way DBB格式及Tao和Guo[24]提出的修正DBB格式計算得到的徑向方向無量綱速度分布與解析解對比.由圖7可知,本文格式的模擬結(jié)果始終與解析解符合良好;而halfway DBB格式和修正DBB格式模擬結(jié)果分別在內(nèi)微圓柱附近和外微圓柱附近存在較大偏差;且Kn數(shù)越大,half-way DBB格式和修正DBB格式在邊界處的偏差越明顯.

    表3給出了不同Kn數(shù)條件下,本文格式同half-way DBB格式和修正DBB格式計算結(jié)果誤差對比情況.從表中可知,本文格式在模擬不同Kn數(shù)時最大相對誤差及相對誤差最小.且由表3可知,修正DBB格式在邊界處最大相對誤差較half-way DBB格式優(yōu)勢并不明顯,但由于在其余區(qū)域符合較好,故在相對誤差方面修正DBB格式比half-way DBB格式計算精度有所提高.

    此外,為檢驗新格式的收斂性能,定義了流場的計算收斂判據(jù):

    其中,u(t+1)(i,j)為當(dāng)前時刻流體格點的速度值,u(t)(i,j)為前一時刻流體格點的速度值.圖8給出了Kn=0.05時,本文格式同half-way DBB格式和修正DBB格式收斂過程對比情況.由圖8可知,三種邊界處理格式的計算收斂時間相近.基于上述計算精度及收斂性能分析比較,本文提出的格式能更好地處理滑移區(qū)運動曲邊界問題.

    圖7 (網(wǎng)刊彩色)微圓柱間Couette流無量綱速度分布(a)Kn=0.05;(b)Kn=0.1Fig.7.(color online)Normalized velocity pro fi les of the microcylinder Couette fl ow:(a)Kn=0.05;(b)Kn=0.1.

    表3 不同Kn數(shù)條件下本文格式與half-way DBB格式及修正DBB格式計算誤差對比Table 3.Computational errors of the present scheme,the half-way DBB scheme and the modi fi ed DBB scheme with di ff erent Kn numbers.

    圖8 (網(wǎng)刊彩色)不同格式的收斂曲線(Kn=0.05)Fig.8.(color online)Convergence pro fi les of di ff erent schemes(Kn=0.05).

    4 結(jié)論

    本文基于非平衡態(tài)外推格式處理曲邊界與有限差分格式獲取曲邊界處宏觀物理量及其法向梯度相結(jié)合的思想,提出了一種適用于滑移區(qū)復(fù)雜氣-固邊界的格子Boltzmann曲邊界處理格式.該格式在處理復(fù)雜邊界時考慮了實際物理邊界與格線偏移量的影響,較half-way DBB格式具有更高的準(zhǔn)確度;且該格式對氣體-壁面相互作用參數(shù)選取無任何限定,相較修正DBB格式有更好的適應(yīng)性.基于該格式對滑移區(qū)氣體平直/傾斜微通道內(nèi)Poiseuille流、微圓柱繞流和同心微圓柱面旋轉(zhuǎn)Couette流問題進(jìn)行計算分析,得到以下結(jié)論:

    1)本文所提邊界處理格式能夠準(zhǔn)確模擬網(wǎng)格線與實際物理邊界不重合的平直以及傾斜微通道外力驅(qū)動Poiseuille流;而half-way DBB格式僅在特定偏移量條件下才能準(zhǔn)確模擬平直微通道內(nèi)Poiseuille流,且在模擬傾斜微通道內(nèi)Poiseuille流時在壁面附近較理論值存在明顯偏離;

    2)通過對微圓柱氣體繞流計算分析,發(fā)現(xiàn)考慮稀薄效應(yīng)時回流區(qū)長度明顯縮短;本文所提格式得到的微圓柱表面無量綱滑移速度與文獻(xiàn)報道結(jié)果更為接近,而half-way DBB格式得到的結(jié)果整體偏小,且最大滑移速度存在較大偏差;

    3)在模擬滑移區(qū)同心微圓柱面Couette流時,本文格式預(yù)測的徑向方向無量綱速度分布與解析解符合良好,模擬精度較half-way DBB格式及修正DBB格式有明顯提高.

    綜上可知,本文提出的曲邊界處理格式具有處理滑移區(qū)復(fù)雜流固邊界速度滑移問題的能力,且較half-way DBB格式和修正DBB格式分別有更高的精度和更好的適應(yīng)性.此外,本文提出的曲邊界處理格式,亦可推廣處理微氣體流動傳熱問題復(fù)雜流-固邊界溫度跳躍邊界條件.

    [1]Gad-el-Hak M 1999 ASME J.Fluids Eng.121 5

    [2]Ho C M,Tai Y C 1998 Annu.Rev.Fluid Mech.30 579

    [3]Zhang T T,Jia L,Wang Z C 2008 Chin.Phys.Lett.25 180

    [4]Sun X M,He F,Ding Y T 2003 Chin.Phys.Lett.20 2199

    [5]Liu C F,Ni Y S 2008 Chin.Phys.B 17 4554

    [6]Tao S,Wang L,Guo Z L 2014 Acta Phys.Sin.63 214703(in Chinese)[陶實,王亮,郭照立2014物理學(xué)報63 214703]

    [7]Wang Z,Liu Y,Zhang J Z 2016 Acta Phys.Sin.65 014703(in Chinese)[王佐,劉雁,張家忠2016物理學(xué)報65 014703]

    [8]Guo Z L,Zheng C G,Shi B C 2002 Phys.Fluids 14 2007

    [9]Mei R W,Luo L S,Shyy W 1999 J.Comput.Phys.155 307

    [10]Ladd A J C 1994 J.Fluid Mech.271 285

    [11]Noble D R,Torczynski J R 1998 Int.J.Mod.Phys.C 9 1189

    [12]Yu Z,Yang H,Fan L S 2011 Chem.Eng.Sci.66 3441[13]Lee C H,Huang Z H,Chiew Y M 2015 Eng.Appl.Comp.Fluid 9 370

    [14]Guo K K,Li L K,Xiao G,AuYeung N A,Mei R W 2015 Int.J.Heat Mass Tran.88 306

    [15]Yin X W,Zhang J F 2012 J.Comput.Phys.231 4295

    [16]Shi D Y,Wang Z K,Zhang A M 2014 Acta Phys.Sin.63 074703(in Chinese)[史冬巖,王志凱,張阿漫2014物理學(xué)報63 074703]

    [17]Le G G,Oulaid O,Zhang J F 2015 Phys.Rev.E 91 033306

    [18]Fu S C,Leung W W F,So R M C 2013 Commun.Comput.Phys.14 126

    [19]Gao D Y,Chen Z Q,Chen L H,Zhang D L 2017 Int.J.Heat Mass Tran.105 673

    [20]Izham M,Fukui T,Morinishi K 2011 J.Fluid Sci.Tech.6 812

    [21]Nie D M,Lin J Z 2010 Chin.Phys.Lett.27 104701

    [22]Fang H P,Wan R Z,Fan L W 2000 Chin.Phys.9 515

    [23]Szalmas L 2007 Int.J.Mod.Phys.C 18 15

    [24]Tao S,Guo Z L 2015 Phys.Rev.E 91 043305

    [25]Guo Z L,Shi B C,Zheng C G 2011 Comput.Math.Appl.61 3519

    [26]Suga K 2013 Fluid Dyn.Res.45 034501

    [27]Yang W L,Li H X,Zhang T J,Elfadel I M 2015 ASME 2015 13th International Conference on Nanochannels,Microchannels,and Minichannels San Francisco,CA,USA,July 6–9,2015 pV001T04A001

    [28]Buick J M,Greated C A 2000 Phys.Rev.E 61 5307

    [29]Guo Z L,Zheng C G,Shi B C 2002 Phys.Rev.E 65 046308

    [30]Sun Y H,Barber R W,Emerson D R 2005 Phys.Fluids 17 047102

    [31]Beskok A,Karniadakis G E 1994 J.Thermophys.Heat Tr.8 647

    PACS:47.11.Qr,47.15.Rq,02.60.CbDOI:10.7498/aps.66.114701

    A new curved boundary treatment in lattice Boltzmann method for micro gas fl ow in the slip regime?

    Gu JuanHuang Rong-ZongLiu Zhen-YuWu Hui-Ying?
    (School of Mechanical Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

    18 December 2016;revised manuscript

    5 April 2017)

    A new curved boundary treatment in lattice Boltzmann method is developed for micro gas fl ow in the slip regime.The proposed treatment is a combination of the nonequilibrium extrapolation scheme for curved boundary with noslip velocity condition and the counter-extrapolation method for the velocity and its normal gradient on the curved boundary.Taking into consideration the e ff ect of the o ff set between the physical boundary and the closest grid line,the new treatment is proved to be more accurate than the traditional half-way di ff usive bounce-back(DBB)scheme.The present treatment is also more applicable than the modi fi ed DBB scheme because the speci fi c gas-wall interaction parameters need to be determined to ensure the validation of the modi fi ed DBB scheme.

    The proposed boundary treatment is implemented to simulate the benchmark problems,which include a Poiseuille fl ow in the aligned/inclined micro-channel,a fl ow past a microcylinder and a microcylindrical Couette fl ow.The results and conclusions are summarized as follows.

    1)The force-driven Poiseuille fl ow in an aligned microchannel is simulated separately with di ff erent values of wallgrid o ff set qδx(q=0.25,0.5,0.75,1.0).With the consideration of the wall-grid o ff set,the numerical results with the new boundary treatment show good agreement with the analytical solutions.However,the results obtained by using the half-way DBB scheme only accord well with the analytical solutions under the condition of a fi xed wall-grid o ff set(q=0.5).

    2)To demonstrate the capability of the present treatment in dealing with gas fl ow in a more complex geometry,the force-driven Poiseuille fl ow in a micro-channel is investigated separately with di ff erent inclined angles.The present numerical results fi t well with the analytical solutions.However,the discrepancy between the results obtained with the half-way DBB scheme and the analytical solutions can be clearly observed near the inclined boundaries.

    3)The gas fl ow past a microcylinder is simulated to prove that the present treatment can deal with the curved boundary.The slip velocity pro fi le along the micro cylinder periphery obtained with the present treatment accords well with the available data in the published literature.However,the results obtained with the half-way DBB scheme show lower values than the results from the published work.

    4)In the simulations of the microcylindrical Couette fl ow between two coaxial rotating cylinders for di ff erent Knudsen numbers the results obtained by using the present treatment show excellent agreement with the analytical solutions.However,the results obtained with the half-way DBB scheme and the modi fi ed DBB scheme deviate obviously from the analytical solutions near the inner and outer cylindrical walls,respectively.

    In summary,the new boundary treatment proposed in this work is capable of dealing with the complex gas-solid boundary in the slip regime.The new treatment has a higher accuracy than the half-way DBB scheme and shows a better applicability than the modi fi ed DBB scheme.

    curved boundary,lattice Boltzmann method,slip regime,micro gas fl ow

    10.7498/aps.66.114701

    ?國家自然科學(xué)基金(批準(zhǔn)號:51536005,51521004)資助的課題.

    ?通信作者.E-mail:whysrj@sjtu.edu.cn

    ?2017中國物理學(xué)會Chinese Physical Society

    http://wulixb.iphy.ac.cn

    *Project supported by the National Natural Science Foundation of China(Grant Nos.51536005,51521004).

    ?Corresponding author.E-mail:whysrj@sjtu.edu.cn

    猜你喜歡
    偏移量無量邊界條件
    烏雷:無量之物
    基于格網(wǎng)坐標(biāo)轉(zhuǎn)換法的矢量數(shù)據(jù)脫密方法研究
    一類帶有Stieltjes積分邊界條件的分?jǐn)?shù)階微分方程邊值問題正解
    帶有積分邊界條件的奇異攝動邊值問題的漸近解
    劉少白
    藝術(shù)品(2020年8期)2020-10-29 02:50:02
    攪拌針不同偏移量對6082-T6鋁合金接頭勞性能的影響
    基于最小二乘平差的全極化SAR配準(zhǔn)偏移量估計方法
    測繪工程(2017年3期)2017-12-22 03:24:50
    論書絕句·評謝無量(1884—1964)
    炳靈寺第70 窟無量壽經(jīng)變辨識
    西藏研究(2017年3期)2017-09-05 09:45:07
    帶Robin邊界條件的2維隨機(jī)Ginzburg-Landau方程的吸引子
    √禁漫天堂资源中文www| av又黄又爽大尺度在线免费看| 99热国产这里只有精品6| 免费不卡的大黄色大毛片视频在线观看| 国产高清国产精品国产三级| 色播在线永久视频| 91精品国产国语对白视频| 亚洲精品国产色婷婷电影| 国语对白做爰xxxⅹ性视频网站| 久久精品久久久久久久性| 少妇精品久久久久久久| 国产一区二区三区综合在线观看| videosex国产| 午夜av观看不卡| 精品一区二区三区四区五区乱码 | 九九爱精品视频在线观看| 日韩 亚洲 欧美在线| 一区二区三区四区激情视频| 日本av手机在线免费观看| 婷婷色综合www| 黄色 视频免费看| 日韩在线高清观看一区二区三区| 亚洲av中文av极速乱| 熟女少妇亚洲综合色aaa.| 高清欧美精品videossex| 韩国精品一区二区三区| 日本免费在线观看一区| 各种免费的搞黄视频| 精品久久久精品久久久| 在线观看一区二区三区激情| 永久免费av网站大全| 最近最新中文字幕免费大全7| 黄色怎么调成土黄色| 曰老女人黄片| 日韩视频在线欧美| 国产免费福利视频在线观看| 一本色道久久久久久精品综合| 日韩成人av中文字幕在线观看| 曰老女人黄片| 精品人妻在线不人妻| 最近中文字幕2019免费版| 日本色播在线视频| av一本久久久久| 国产精品久久久久成人av| 久久久久久人人人人人| 久久久久精品人妻al黑| 男女边摸边吃奶| 亚洲精品乱久久久久久| 最近手机中文字幕大全| 精品亚洲成a人片在线观看| 国产探花极品一区二区| 如日韩欧美国产精品一区二区三区| 精品一区二区三区四区五区乱码 | 亚洲久久久国产精品| 精品人妻偷拍中文字幕| 亚洲精品国产一区二区精华液| 99久久综合免费| 亚洲精品在线美女| 国产毛片在线视频| 国产精品无大码| 欧美少妇被猛烈插入视频| 99久久精品国产国产毛片| 亚洲精品国产av蜜桃| 亚洲精品一二三| 大话2 男鬼变身卡| 日韩av免费高清视频| 咕卡用的链子| 老女人水多毛片| 看十八女毛片水多多多| 最新中文字幕久久久久| 亚洲成人手机| 久久久久精品性色| 久久久久网色| 午夜激情av网站| 亚洲欧美清纯卡通| 好男人视频免费观看在线| 国产精品一区二区在线不卡| www.熟女人妻精品国产| 777米奇影视久久| 国产亚洲最大av| 18禁裸乳无遮挡动漫免费视频| 欧美日韩亚洲国产一区二区在线观看 | 婷婷色综合大香蕉| 精品一区二区三区四区五区乱码 | 又黄又粗又硬又大视频| 极品少妇高潮喷水抽搐| 在线天堂最新版资源| 又大又黄又爽视频免费| 免费av中文字幕在线| 免费日韩欧美在线观看| 男人添女人高潮全过程视频| 国产精品 国内视频| av又黄又爽大尺度在线免费看| 一区二区三区乱码不卡18| 免费黄网站久久成人精品| 国产精品二区激情视频| 制服诱惑二区| 男男h啪啪无遮挡| 如何舔出高潮| 久久精品人人爽人人爽视色| 一级毛片 在线播放| 精品少妇内射三级| 精品99又大又爽又粗少妇毛片| 亚洲成av片中文字幕在线观看 | 久久ye,这里只有精品| 亚洲av免费高清在线观看| 国产亚洲精品第一综合不卡| 制服丝袜香蕉在线| 欧美另类一区| 超碰97精品在线观看| 黄片小视频在线播放| 婷婷色麻豆天堂久久| 多毛熟女@视频| 午夜福利视频精品| 国产亚洲精品第一综合不卡| 国产成人精品无人区| 国产精品人妻久久久影院| 一级毛片我不卡| 97在线人人人人妻| 国产福利在线免费观看视频| 欧美另类一区| 丝袜在线中文字幕| 90打野战视频偷拍视频| 女性生殖器流出的白浆| 一区二区av电影网| 春色校园在线视频观看| 国产亚洲精品第一综合不卡| 成人毛片a级毛片在线播放| 亚洲综合色网址| 超碰成人久久| 另类亚洲欧美激情| 91午夜精品亚洲一区二区三区| 国产精品欧美亚洲77777| 午夜91福利影院| 久久精品久久久久久久性| 国产欧美日韩一区二区三区在线| 九九爱精品视频在线观看| 精品人妻一区二区三区麻豆| 极品少妇高潮喷水抽搐| 亚洲精品一区蜜桃| 男女啪啪激烈高潮av片| 中文天堂在线官网| 日韩一卡2卡3卡4卡2021年| 久久久久久久久免费视频了| 黄色视频在线播放观看不卡| videos熟女内射| 日韩熟女老妇一区二区性免费视频| 国产 一区精品| 热re99久久国产66热| 国产日韩欧美在线精品| 欧美日韩视频精品一区| 不卡av一区二区三区| 国产精品国产av在线观看| 一区在线观看完整版| 国语对白做爰xxxⅹ性视频网站| 秋霞伦理黄片| av福利片在线| 欧美日韩精品网址| 九草在线视频观看| 国产激情久久老熟女| 天堂俺去俺来也www色官网| 久久精品人人爽人人爽视色| 国产高清国产精品国产三级| 久久久a久久爽久久v久久| 一级,二级,三级黄色视频| 看免费av毛片| 精品一品国产午夜福利视频| 天天躁夜夜躁狠狠久久av| 国产成人91sexporn| 好男人视频免费观看在线| 一个人免费看片子| 亚洲色图 男人天堂 中文字幕| 国产精品久久久久久av不卡| 久久久久久久亚洲中文字幕| 日本午夜av视频| 夜夜骑夜夜射夜夜干| 两个人免费观看高清视频| 亚洲色图 男人天堂 中文字幕| 天天躁日日躁夜夜躁夜夜| 欧美在线黄色| 精品人妻熟女毛片av久久网站| 性色av一级| 宅男免费午夜| 制服诱惑二区| 毛片一级片免费看久久久久| 国产乱人偷精品视频| 丝袜美腿诱惑在线| 国产一区亚洲一区在线观看| 亚洲国产欧美网| 男女午夜视频在线观看| 国产成人精品福利久久| av在线观看视频网站免费| 成人影院久久| 最近2019中文字幕mv第一页| 国产亚洲精品第一综合不卡| 色网站视频免费| 人人妻人人澡人人看| 观看美女的网站| 啦啦啦啦在线视频资源| 午夜免费男女啪啪视频观看| 欧美精品一区二区大全| 各种免费的搞黄视频| 一本大道久久a久久精品| 日日啪夜夜爽| 黑丝袜美女国产一区| 精品一区二区免费观看| 亚洲国产最新在线播放| 99久久人妻综合| 99国产综合亚洲精品| 一级毛片电影观看| 少妇精品久久久久久久| 亚洲欧洲国产日韩| 丰满饥渴人妻一区二区三| 国产有黄有色有爽视频| 黑人欧美特级aaaaaa片| 少妇精品久久久久久久| 精品国产一区二区久久| 日韩中文字幕欧美一区二区 | 久久久久久免费高清国产稀缺| 亚洲成人av在线免费| 男女国产视频网站| 午夜影院在线不卡| 少妇人妻久久综合中文| 老汉色∧v一级毛片| 亚洲一码二码三码区别大吗| 午夜日韩欧美国产| 欧美人与善性xxx| 又粗又硬又长又爽又黄的视频| 欧美成人精品欧美一级黄| 亚洲av日韩在线播放| 啦啦啦中文免费视频观看日本| 精品国产一区二区三区四区第35| 男人操女人黄网站| 97精品久久久久久久久久精品| 黑丝袜美女国产一区| 成人18禁高潮啪啪吃奶动态图| 国产日韩欧美亚洲二区| 日韩精品有码人妻一区| 丰满饥渴人妻一区二区三| 制服丝袜香蕉在线| 丰满少妇做爰视频| 精品酒店卫生间| 国产精品欧美亚洲77777| 97在线人人人人妻| 免费播放大片免费观看视频在线观看| 婷婷色av中文字幕| 最新的欧美精品一区二区| 最近手机中文字幕大全| 成年人免费黄色播放视频| 多毛熟女@视频| 极品少妇高潮喷水抽搐| 亚洲欧美成人综合另类久久久| 99热全是精品| 满18在线观看网站| 国产精品久久久久久av不卡| 亚洲精品美女久久av网站| 免费av中文字幕在线| 国产乱人偷精品视频| 超色免费av| 欧美国产精品va在线观看不卡| 久久鲁丝午夜福利片| 高清av免费在线| 美女视频免费永久观看网站| 日韩av在线免费看完整版不卡| 精品国产国语对白av| 亚洲欧美成人精品一区二区| 国产激情久久老熟女| 97在线人人人人妻| 精品国产一区二区久久| 999久久久国产精品视频| 欧美亚洲 丝袜 人妻 在线| 免费在线观看视频国产中文字幕亚洲 | 99热网站在线观看| 久热久热在线精品观看| 国产一区二区三区av在线| 亚洲欧美中文字幕日韩二区| 欧美xxⅹ黑人| 少妇猛男粗大的猛烈进出视频| 1024香蕉在线观看| 丰满乱子伦码专区| 国产日韩欧美亚洲二区| 日韩一卡2卡3卡4卡2021年| 丝袜在线中文字幕| 国产 一区精品| 精品少妇一区二区三区视频日本电影 | 国产一区二区三区综合在线观看| 久久久久久久久免费视频了| 国产伦理片在线播放av一区| 制服人妻中文乱码| 午夜福利一区二区在线看| 各种免费的搞黄视频| 韩国精品一区二区三区| 最新中文字幕久久久久| 亚洲精品国产av蜜桃| 国产精品成人在线| 丰满饥渴人妻一区二区三| 日韩一本色道免费dvd| 欧美国产精品va在线观看不卡| 国产男女超爽视频在线观看| 久久午夜福利片| 日韩,欧美,国产一区二区三区| 欧美国产精品一级二级三级| 日韩一区二区视频免费看| 国产成人免费观看mmmm| 亚洲精品一区蜜桃| 99热国产这里只有精品6| 热re99久久精品国产66热6| 精品人妻偷拍中文字幕| 五月天丁香电影| 香蕉丝袜av| 美女福利国产在线| 看十八女毛片水多多多| 99久久综合免费| 岛国毛片在线播放| 黄片播放在线免费| 十八禁网站网址无遮挡| 观看美女的网站| 亚洲精华国产精华液的使用体验| 亚洲精品国产色婷婷电影| 精品一品国产午夜福利视频| 国产精品欧美亚洲77777| 精品一区二区三卡| 国产亚洲欧美精品永久| 美女午夜性视频免费| 国产97色在线日韩免费| 亚洲国产毛片av蜜桃av| 亚洲精品乱久久久久久| 国产免费又黄又爽又色| 熟女少妇亚洲综合色aaa.| 天天影视国产精品| 亚洲综合色网址| 一级爰片在线观看| 国产一区亚洲一区在线观看| 精品一区二区三卡| 国产精品三级大全| 热99久久久久精品小说推荐| 成人漫画全彩无遮挡| 亚洲精华国产精华液的使用体验| 在线观看国产h片| 美女中出高潮动态图| xxx大片免费视频| 亚洲国产毛片av蜜桃av| 少妇的丰满在线观看| 乱人伦中国视频| 99热国产这里只有精品6| 晚上一个人看的免费电影| 国产精品国产三级专区第一集| 一二三四在线观看免费中文在| 欧美精品人与动牲交sv欧美| 亚洲欧美精品综合一区二区三区 | 亚洲欧美清纯卡通| 日本-黄色视频高清免费观看| 婷婷色综合www| 国产有黄有色有爽视频| 在现免费观看毛片| 国产精品 国内视频| 午夜激情av网站| 母亲3免费完整高清在线观看 | 国产高清不卡午夜福利| 亚洲欧美一区二区三区久久| av.在线天堂| 久久国内精品自在自线图片| 永久网站在线| 亚洲欧美中文字幕日韩二区| 国产黄色免费在线视频| 国产无遮挡羞羞视频在线观看| 黄片无遮挡物在线观看| 丝袜在线中文字幕| 国产极品粉嫩免费观看在线| 亚洲av在线观看美女高潮| 精品一区二区三区四区五区乱码 | 午夜福利影视在线免费观看| 久久精品国产亚洲av高清一级| av不卡在线播放| 满18在线观看网站| 久久久久国产精品人妻一区二区| 18在线观看网站| 麻豆乱淫一区二区| 欧美精品一区二区免费开放| 亚洲国产成人一精品久久久| 日韩一卡2卡3卡4卡2021年| 精品亚洲成a人片在线观看| 日韩大片免费观看网站| 精品视频人人做人人爽| 欧美精品av麻豆av| 精品人妻熟女毛片av久久网站| 亚洲色图综合在线观看| 亚洲人成77777在线视频| 午夜老司机福利剧场| 人体艺术视频欧美日本| 国产精品熟女久久久久浪| 成人毛片60女人毛片免费| 中文乱码字字幕精品一区二区三区| 男女下面插进去视频免费观看| 久久久久网色| 汤姆久久久久久久影院中文字幕| 婷婷色av中文字幕| 日本wwww免费看| 涩涩av久久男人的天堂| 高清av免费在线| 考比视频在线观看| 色视频在线一区二区三区| 十八禁网站网址无遮挡| 国产免费一区二区三区四区乱码| 亚洲欧美精品自产自拍| 亚洲欧美精品综合一区二区三区 | 亚洲国产色片| 欧美亚洲 丝袜 人妻 在线| 妹子高潮喷水视频| 一区二区av电影网| 日本欧美视频一区| 伦理电影免费视频| 久久久久精品人妻al黑| 久久精品国产综合久久久| 中文字幕人妻丝袜制服| 日韩 亚洲 欧美在线| 日产精品乱码卡一卡2卡三| 亚洲av综合色区一区| av福利片在线| 丝袜喷水一区| 国产精品 欧美亚洲| 国产又色又爽无遮挡免| 亚洲 欧美一区二区三区| 女性生殖器流出的白浆| 超碰成人久久| 国产黄色免费在线视频| 老汉色∧v一级毛片| 91午夜精品亚洲一区二区三区| 夜夜骑夜夜射夜夜干| 亚洲精品aⅴ在线观看| 99久国产av精品国产电影| 成人免费观看视频高清| 久久 成人 亚洲| 国产高清不卡午夜福利| 9热在线视频观看99| 亚洲激情五月婷婷啪啪| 女人被躁到高潮嗷嗷叫费观| 亚洲伊人久久精品综合| 精品视频人人做人人爽| 欧美日韩视频高清一区二区三区二| 五月开心婷婷网| 又黄又粗又硬又大视频| 亚洲av国产av综合av卡| 18禁动态无遮挡网站| 99热网站在线观看| 人妻 亚洲 视频| 丰满迷人的少妇在线观看| 久久鲁丝午夜福利片| 咕卡用的链子| 一级片免费观看大全| 人妻系列 视频| 成人黄色视频免费在线看| 亚洲av在线观看美女高潮| xxx大片免费视频| 日本vs欧美在线观看视频| 国产色婷婷99| 在线观看国产h片| 精品酒店卫生间| 中文字幕精品免费在线观看视频| 国产黄色免费在线视频| 久久人人97超碰香蕉20202| 中国三级夫妇交换| 亚洲精华国产精华液的使用体验| 91精品国产国语对白视频| 黄频高清免费视频| 黄网站色视频无遮挡免费观看| 日本午夜av视频| 丁香六月天网| 97在线视频观看| 亚洲久久久国产精品| videossex国产| 亚洲av欧美aⅴ国产| 久久久精品区二区三区| 亚洲综合色惰| 人妻 亚洲 视频| 亚洲国产最新在线播放| 制服丝袜香蕉在线| 啦啦啦啦在线视频资源| 桃花免费在线播放| 少妇被粗大的猛进出69影院| 亚洲,一卡二卡三卡| 国产精品av久久久久免费| a级毛片黄视频| 在线观看免费日韩欧美大片| 女人高潮潮喷娇喘18禁视频| 伦理电影免费视频| 亚洲经典国产精华液单| 欧美日韩av久久| 哪个播放器可以免费观看大片| 亚洲少妇的诱惑av| 国产精品 欧美亚洲| 高清视频免费观看一区二区| 午夜91福利影院| 男女边摸边吃奶| 亚洲精品国产av蜜桃| 国产男人的电影天堂91| 黄片小视频在线播放| 精品卡一卡二卡四卡免费| 伊人亚洲综合成人网| 老汉色∧v一级毛片| 一级毛片 在线播放| 日韩欧美一区视频在线观看| 校园人妻丝袜中文字幕| 国产成人精品久久久久久| 精品国产超薄肉色丝袜足j| 精品国产一区二区久久| av视频免费观看在线观看| 中文字幕av电影在线播放| 最近最新中文字幕大全免费视频 | 最近2019中文字幕mv第一页| 26uuu在线亚洲综合色| 巨乳人妻的诱惑在线观看| 亚洲国产精品一区三区| 色网站视频免费| 一二三四在线观看免费中文在| 免费在线观看黄色视频的| 日韩视频在线欧美| 亚洲 欧美一区二区三区| 男女免费视频国产| 国产精品蜜桃在线观看| 大码成人一级视频| 亚洲,欧美,日韩| av电影中文网址| 免费人妻精品一区二区三区视频| 久久精品国产鲁丝片午夜精品| 激情五月婷婷亚洲| 人人妻人人澡人人爽人人夜夜| 丝袜脚勾引网站| av免费观看日本| 亚洲国产看品久久| 国产av国产精品国产| 制服丝袜香蕉在线| 丝袜在线中文字幕| 亚洲精品国产一区二区精华液| 黄色配什么色好看| 男人爽女人下面视频在线观看| 美女主播在线视频| 午夜福利乱码中文字幕| 热re99久久国产66热| 日韩av在线免费看完整版不卡| 久久精品国产a三级三级三级| 国产片内射在线| 国产高清不卡午夜福利| 三级国产精品片| 成人国语在线视频| 亚洲在久久综合| 国产老妇伦熟女老妇高清| av福利片在线| 午夜福利一区二区在线看| 好男人视频免费观看在线| 汤姆久久久久久久影院中文字幕| 在线观看人妻少妇| 黑人欧美特级aaaaaa片| 亚洲国产精品一区二区三区在线| 中文字幕制服av| www.精华液| 午夜福利在线免费观看网站| 亚洲精品第二区| 国产97色在线日韩免费| 亚洲在久久综合| 美女国产高潮福利片在线看| av福利片在线| 国产精品国产av在线观看| 性色avwww在线观看| 国产一区二区三区av在线| 国产又爽黄色视频| xxxhd国产人妻xxx| 母亲3免费完整高清在线观看 | 欧美人与性动交α欧美精品济南到 | 寂寞人妻少妇视频99o| 国产成人欧美| 夫妻午夜视频| 波野结衣二区三区在线| 丰满饥渴人妻一区二区三| 亚洲av福利一区| 女性生殖器流出的白浆| 久久人人爽人人片av| 各种免费的搞黄视频| 看免费成人av毛片| 国产 精品1| 亚洲伊人久久精品综合| 国产精品二区激情视频| 久久久久精品性色| 性少妇av在线| 99re6热这里在线精品视频| 五月天丁香电影| 一区二区av电影网| 午夜福利,免费看| 国产片特级美女逼逼视频| 亚洲欧美成人综合另类久久久| 宅男免费午夜| www日本在线高清视频| 搡女人真爽免费视频火全软件| a级片在线免费高清观看视频| 日韩免费高清中文字幕av| 午夜av观看不卡| 久久狼人影院| 丰满少妇做爰视频| 免费观看a级毛片全部| 久久这里有精品视频免费| 好男人视频免费观看在线| 91aial.com中文字幕在线观看| h视频一区二区三区| 王馨瑶露胸无遮挡在线观看| 如日韩欧美国产精品一区二区三区| 日韩一区二区视频免费看| 午夜福利影视在线免费观看| av福利片在线| 欧美激情 高清一区二区三区| a级毛片黄视频| 日日撸夜夜添| 岛国毛片在线播放| 青春草国产在线视频| 777米奇影视久久| 亚洲人成77777在线视频|