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

    譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程

    2017-11-07 01:11:33時(shí)永興
    綜合智慧能源 2017年10期
    關(guān)鍵詞:方法

    時(shí)永興

    (華電江蘇能源有限公司句容發(fā)電廠,江蘇 鎮(zhèn)江 212400)

    譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程

    時(shí)永興

    (華電江蘇能源有限公司句容發(fā)電廠,江蘇 鎮(zhèn)江 212400)

    針對(duì)譜元方法求解高雷諾數(shù)下對(duì)流擴(kuò)散方程的穩(wěn)定性問題,采用Chebyshev譜元方法結(jié)合譜消去粘性法求解一維對(duì)流擴(kuò)散方程。利用特征分析法預(yù)測(cè)了數(shù)值方法的求解穩(wěn)定性,通過數(shù)值算例驗(yàn)證了該解法的可行性,討論了譜消去粘性參數(shù)對(duì)求解穩(wěn)定性及數(shù)值精度的影響。結(jié)果表明:和譜元方法相比,譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程的穩(wěn)定區(qū)域有了明顯的擴(kuò)大,在高雷諾數(shù)時(shí)能夠獲得具有較高精度的數(shù)值解;較大的譜消去粘性項(xiàng)有利于穩(wěn)定區(qū)域的擴(kuò)大,而在計(jì)算穩(wěn)定的條件下較小的粘性項(xiàng)有利于數(shù)值精度的提高,所以適當(dāng)?shù)卦O(shè)置粘性項(xiàng)的大小,在保證計(jì)算穩(wěn)定的同時(shí)提高數(shù)值精度。

    對(duì)流擴(kuò)散方程;譜元法;譜消去粘性法;穩(wěn)定性;高雷諾數(shù)

    0 引言

    對(duì)流擴(kuò)散方程一般用來描述濃度、溫度、渦量等物理量的轉(zhuǎn)移擴(kuò)散過程,在環(huán)境、化工、氣象、水利等多個(gè)工程領(lǐng)域中都有著廣泛的應(yīng)用,因其重要性而備受關(guān)注。由于一般情況下難以求得解析解,高精度的數(shù)值解對(duì)于準(zhǔn)確地理解物理過程中的對(duì)流擴(kuò)散現(xiàn)象就有著重要的意義。

    譜元方法[1]SEM (spectral element method)作為一種靈活的高精度數(shù)值方法,結(jié)合了譜方法的高精度和有限元法的復(fù)雜區(qū)域處理能力,非常適宜應(yīng)用于對(duì)流擴(kuò)散方程的數(shù)值求解。但當(dāng)方程中含有很高的雷諾數(shù)Re時(shí),應(yīng)用譜元方法進(jìn)行求解會(huì)出現(xiàn)傳統(tǒng)的有限差分法[2]、有限元法[3]、譜方法[4]所遇到的穩(wěn)定性問題,即對(duì)于一定的時(shí)間步長(zhǎng)和網(wǎng)格劃分存在Re數(shù)區(qū)間,只有當(dāng)計(jì)算位于該Re數(shù)區(qū)間時(shí),求解才是穩(wěn)定的。因此,擴(kuò)大穩(wěn)定區(qū)域的范圍成為譜元方法求解對(duì)流擴(kuò)散方程的關(guān)鍵問題。

    譜消去粘性法SVV(spectral vanishing viscosity)是由Tadmar[5]在利用Fourier譜方法求解雙曲型守恒方程時(shí)首先提出的,Maday[6]等將其擴(kuò)展到求解非周期問題的Legendre譜方法獲得了守恒方程穩(wěn)定的數(shù)值解;Heinrichs[7]利用譜消去粘性譜方法求解了含有邊界層的穩(wěn)態(tài)對(duì)流擴(kuò)散方程,消除了邊界層所造成的數(shù)值振蕩;Karamanos和Karniadakis[8]首次將譜消去粘性法和譜元方法相結(jié)合,成功地計(jì)算了管道湍流;Xu和Pasquetti[9]將譜消去粘性法應(yīng)用于配置點(diǎn)譜元方法,推廣了其求解高維復(fù)雜區(qū)域問題的形式,并獲得了高Re數(shù)下圓柱繞流的穩(wěn)定數(shù)值解;Kirby和Sherwin[10]利用矩陣形式分析了譜消去粘性法和譜元方法的結(jié)合方式,對(duì)三維三角形管道湍流進(jìn)行了求解。

    本文在前期譜消去粘性譜元法研究的基礎(chǔ)上,采用矩陣分析法提出求解過程的穩(wěn)定性條件,系統(tǒng)地研究譜消去粘性譜元法求解對(duì)流擴(kuò)散方程的穩(wěn)定性問題及其擴(kuò)大穩(wěn)定域的效果。討論譜消去粘性參數(shù)對(duì)穩(wěn)定性和計(jì)算精度的影響,通過數(shù)值算例對(duì)譜消去粘性譜元法求解高Re數(shù)下對(duì)流擴(kuò)散方程的有效性進(jìn)行驗(yàn)證。

    1 譜消去粘性法

    譜消去粘性法求解對(duì)流擴(kuò)散方程時(shí),在方程中加入人工粘性項(xiàng),穩(wěn)定化的一維對(duì)流擴(kuò)散方程形式為:

    00 ,

    (1)

    相應(yīng)的初始條件和邊界條件:

    u(x,0)=u0(x),0

    (2)

    u(0,t)=g1(t),u(1,t)=g2(t),t>0

    (3)

    ε,Q依賴于用式(1)求解的多項(xiàng)式的階數(shù)N,ε通常取ε≈N-1,Q在一維標(biāo)準(zhǔn)區(qū)域(-1,1)內(nèi)的定義為:

    (4)

    2 Chebyshev譜元離散

    采用譜元方法求解式(1)~式(3)即為尋找u∈H1,使得在時(shí)間區(qū)間[0,1]上滿足:

    ε(Qu,v)=(f,v),?v∈H1,

    (5)

    在邊界上滿足u(0,t)=g1(t),u(1,t)=g2(t),其中(·,·)表示定義在[0,1]上的L2內(nèi)積。

    由式(5)可得,譜消去粘性項(xiàng)的弱形式為:

    Svv(u,v)=ε(Qu,v) ,

    (6)

    為了保持其系數(shù)矩陣的對(duì)稱性,可改寫為:

    Svv(u,v)=ε(Q1/2u,Q1/2v) 。

    (7)

    由Chebyshev多項(xiàng)式的性質(zhì)可知,式(6)、式(7)所表示的譜消去粘性項(xiàng)的弱形式并不完全相同,但式(7)保證了系數(shù)矩陣的對(duì)稱性,可在計(jì)算中近似地采用。

    實(shí)際計(jì)算時(shí),將式(5)中的粘性項(xiàng)和譜消去粘性項(xiàng)相合并,則可以得到:

    (f,v),?v∈H1,

    (8)

    式中:S=I+εReQ,I為單位矩陣。

    (9)

    (10)

    在標(biāo)準(zhǔn)單元Λ內(nèi)對(duì)式(8)進(jìn)行離散,則可以得到第i個(gè)單元的譜元離散形式:

    (11)

    (12)

    初始條件為:u(x,0)=u0(x)。

    3 時(shí)間離散

    在時(shí)域上求解微分方程組一般有顯式和隱式兩種積分法。隱式積分法絕對(duì)穩(wěn)定,但每一個(gè)時(shí)間步都需要求解線性方程組,顯式積分法計(jì)算高效,需要的計(jì)算機(jī)內(nèi)存小,但一般條件穩(wěn)定。本文采用半隱式中心差分格式對(duì)時(shí)間項(xiàng)進(jìn)行離散,對(duì)流項(xiàng)采用顯式處理,粘性項(xiàng)采用隱式處理,融合了顯式方法計(jì)算簡(jiǎn)單和隱式方法穩(wěn)定性好的特點(diǎn),式(12)的全離散形式為:

    c1un+1=c2un+c3un-1+Bfn,

    (13)

    4 穩(wěn)定性分析

    4.1譜消去粘性譜元方法穩(wěn)定性條件

    將式(13)兩側(cè)左乘c1的逆矩陣得:

    (14)

    由于僅考慮第1類邊界條件,即在邊值的計(jì)算中不引入誤差,所有的誤差來自初始值的迭代過程。假定在初值的計(jì)算中引入了擾動(dòng)τ0,則擾動(dòng)τ滿足:

    τn+1=M1τn+M2τn-1,

    (15)

    Dorsselaer等[11]分析了形如式(15)的全離散格式的穩(wěn)定性條件,定義擾動(dòng)Vn+1=((τn+1)T,(τn)T)T,其滿足:

    Vn+1=EVn,n≥1 ,

    (16)

    可以驗(yàn)證當(dāng)采用Chebyshev函數(shù)為單元基函數(shù)時(shí),擾動(dòng)系數(shù)矩陣E為正規(guī)矩陣。譜消去粘性譜元法求解穩(wěn)定的充分條件為:矩陣E的譜半徑ρ(E)≤1。

    時(shí)永興:譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程

    4.2穩(wěn)定性結(jié)果

    圖1 不同時(shí)間步長(zhǎng)的穩(wěn)定性曲線

    將粘性振幅設(shè)置為ε=α/Nx,通過改變自由參數(shù)α的值,來研究ε的大小對(duì)求解穩(wěn)定性的影響。圖2給出了不同粘性振幅譜消去粘性譜元方法的求解穩(wěn)定域??梢钥闯?,對(duì)于較大的時(shí)間步長(zhǎng),ε值的增加對(duì)求解穩(wěn)定域的擴(kuò)大作用有限;對(duì)于較小的時(shí)間步長(zhǎng),求解穩(wěn)定域隨著ε值的增加不斷擴(kuò)大。說明對(duì)于較小的時(shí)間步長(zhǎng),較大的譜消去粘性項(xiàng)有利于求解穩(wěn)定域的擴(kuò)大。圖2中截?cái)嚯A數(shù)的設(shè)置為mN=Nx/2。

    圖2 不同粘性振幅的穩(wěn)定性曲線

    圖3給出了不同截?cái)嚯A數(shù)譜消去粘性譜元方法的求解穩(wěn)定域。可以看出,對(duì)于2種時(shí)間步長(zhǎng),當(dāng)截?cái)嚯A數(shù)mN從Nx減小到2時(shí),求解穩(wěn)定域基本上保持不變;當(dāng)截?cái)嚯A數(shù)mN取0或1時(shí),求解穩(wěn)定域顯著增加,同樣說明較大的譜消去粘性項(xiàng)有利于求解穩(wěn)定域的擴(kuò)大。圖3中粘性振幅的設(shè)置為ε=5.0/Nx。

    圖3 不同截?cái)嚯A數(shù)的穩(wěn)定性曲線

    5 算例及分析

    當(dāng)f=0時(shí),利用變量分離法可得式(1)的一個(gè)解析解為:

    (17)

    相應(yīng)的初始條件、邊界條件:

    u(x,0)=0, 0≤x≤1 ,

    (18)

    u(0,t)=0,u(1,t)=1 ,t>0 ,

    (19)

    該算例在邊界x=1處有一個(gè)流動(dòng)邊界層的存在,當(dāng)Re數(shù)增大時(shí),邊界層的厚度逐漸變小,速度變化劇烈,利用譜元方法進(jìn)行計(jì)算時(shí)會(huì)產(chǎn)生數(shù)值振蕩。

    本文采用ea和L2誤差來考察n時(shí)間層上數(shù)值解的準(zhǔn)確性,兩種誤差的定義為:

    (20)

    (21)

    通過數(shù)值解和解析解的比較來驗(yàn)證譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程的有效性。

    5.1譜元方法和譜消去粘性譜元方法的精度比較

    選取網(wǎng)格劃分Nm=5,Nx=7,時(shí)間步長(zhǎng)Δt=0.01,計(jì)算時(shí)間t=1.0,粘性參數(shù)ε=0.012/Nx,mN=Nx/2。由穩(wěn)定性條件可以計(jì)算出對(duì)于所選定的時(shí)間步長(zhǎng)和網(wǎng)格劃分譜元方法求解的臨界Re=557。譜元方法和譜消去粘性譜元方法在Re=570求解時(shí),不同節(jié)點(diǎn)的絕對(duì)誤差如圖4所示??梢钥闯觯诳拷吔鐚拥膮^(qū)域譜元方法的數(shù)值結(jié)果已經(jīng)開始發(fā)散,而譜消去粘性譜元方法的數(shù)值結(jié)果仍然和解析解相吻合。

    圖4 譜元方法和譜消去粘性譜元方法的絕對(duì)誤差

    譜元方法和譜消去粘性譜元方法在不同Re數(shù)的L2誤差如圖5所示。可以看出,譜元方法的L2誤差隨著Re數(shù)的增加逐步增大,特別是當(dāng)Re大于臨界值時(shí)數(shù)值結(jié)果快速發(fā)散;譜消去粘性譜元方法的L2誤差隨Re數(shù)的增加基本上保持不變,即使在穩(wěn)定區(qū)域較小的譜消去粘性項(xiàng)的增加仍有利于計(jì)算精度的提高。

    圖5 譜元方法和譜消去粘性譜元方法的L2誤差

    5.2粘性參數(shù)對(duì)譜消去粘性譜元方法的影響

    圖6給出了譜消去粘性譜元方法不同ε值的L2誤差。選取2種不同的網(wǎng)格劃分Nm=4、Nx=4,Nm=8、Nx=8,研究粘性參數(shù)的變化對(duì)兩種網(wǎng)格計(jì)算精度的影響。網(wǎng)格劃分Nm=4,Nx=4的時(shí)間步長(zhǎng)為Δt=0.05,相應(yīng)的譜元求解臨界Re=100,計(jì)算Re=150;網(wǎng)格劃分Nm=8,Nx=8的時(shí)間步長(zhǎng)為Δt=0.005,相應(yīng)的譜元求解臨界Re=874,計(jì)算Re=1 000。計(jì)算時(shí)間t=1.0,截?cái)嚯A數(shù)mN=Nx/2??梢钥闯?,在保證計(jì)算穩(wěn)定的條件下,較小的粘性振幅有利于計(jì)算精度的提高,同時(shí)對(duì)于稀疏的網(wǎng)格,保證其計(jì)算穩(wěn)定所需要的譜消去粘性項(xiàng)較大,即粘性振幅較大。

    圖6 譜消去粘性譜元方法不同ε的L2誤差

    圖7給出了譜消去粘性譜元方法不同mN值的L2誤差。選取網(wǎng)格劃分Nm=4、Nx=6,Nm=6、Nx=6,研究截?cái)嚯A數(shù)的變化對(duì)兩種網(wǎng)格計(jì)算精度的影響。網(wǎng)格劃分Nm=4、Nx=6的時(shí)間步長(zhǎng)為Δt=0.02,相應(yīng)的譜元求解臨界Re=270,計(jì)算Re=350;網(wǎng)格劃分Nm=6、Nx=6的時(shí)間步長(zhǎng)為Δt=0.01,相應(yīng)的譜元求解臨界Re=500,計(jì)算Re=600。計(jì)算時(shí)間t=1.0,粘性振幅ε=0.1/Nx??梢钥闯?,截?cái)嚯A數(shù)的變化對(duì)于計(jì)算精度的影響不是十分明顯,稀疏的網(wǎng)格需要較小的截?cái)嚯A數(shù)以保證計(jì)算的穩(wěn)定。對(duì)于Nm=4、Nx=6的網(wǎng)格其計(jì)算在mN=Nx-1時(shí)已經(jīng)發(fā)散,而對(duì)于Nm=6、Nx=6的網(wǎng)格,其所需的譜消去粘性項(xiàng)較小,在mN=Nx-1時(shí)計(jì)算誤差反而達(dá)到最小。

    圖7 譜消去粘性譜元方法不同mN的L2誤差

    6 結(jié)論

    本文采用Chebyshev譜元方法耦合譜消去粘性法求解了一維對(duì)流擴(kuò)散方程,對(duì)該方法的求解穩(wěn)定性和計(jì)算精度進(jìn)行了分析,討論了譜消去粘性參數(shù)對(duì)穩(wěn)定性及計(jì)算精度的影響,并與譜元方法的求解結(jié)果進(jìn)行了對(duì)比,可以得出如下結(jié)論。

    (1)譜元方法在求解含有高Re數(shù)的對(duì)流擴(kuò)散方程時(shí)穩(wěn)定區(qū)域十分有限,譜粘消去粘性項(xiàng)的增加有效地?cái)U(kuò)大了穩(wěn)定區(qū)域的范圍。

    (2)在高Re數(shù)時(shí)譜消去粘性項(xiàng)保證了譜元求解的數(shù)值精度,即使在穩(wěn)定區(qū)域較小的譜消去粘性項(xiàng)的增加仍然有利于數(shù)值精度的提高。

    (3)粘性參數(shù)的設(shè)置對(duì)于譜消去粘性譜元方法的穩(wěn)定區(qū)域的擴(kuò)大和計(jì)算精度的提高都有著非常大

    時(shí)永興:譜消去粘性譜元方法求解對(duì)流擴(kuò)散方程

    的影響。較大的譜消去粘性項(xiàng)有利于穩(wěn)定域的擴(kuò)大,但在保證計(jì)算穩(wěn)定的條件下,較小的粘性項(xiàng)又有利于精度的提高。

    今后的工作可以從以下兩個(gè)方面展開:應(yīng)用譜元方法求解非線性對(duì)流擴(kuò)散方程,并對(duì)求解過程的穩(wěn)定性和影響因素進(jìn)行分析;將譜消去粘性譜元方法拓展至三維非定常流動(dòng)。

    [1]PETERA A T.A spectral element method for fluid dynamics:laminar flow in a channel expansion[J].J comp phys,1984,54(3):468-488.

    [2]陶文銓.數(shù)值傳熱學(xué)[M].2版.西安:西安交通大學(xué)出版社,2001:48-186.

    [3]MOHAMED S A,MOHSMED N A,ABDEL GAWAD A F,et al.A modified diffusion coefficient technique for the convection diffusion equation[J].Appl math comput,2013,219:9317-9330.

    [4]MOFID A,PEYRET R.Stability of the Chebyshev collocation approximation to the advection-diffusion equation[J].Comput fluids,1993,22(4-5):453-465.

    [5]TADMOR E.Convergence of spectral methods for nonlinear conservation laws[J].SIAM J numer anal,1989,26(1):30-44.

    [6]MADAY Y,OULD KABER S M,TADMOR E.Legendre pseudospectral viscosity method for nonlinear conservation laws[J].SIAM J numer anal,1993,30(2):321-342.

    [7]HEINRICHS W.Spectral viscosity for convection dominated flow[J].J sci comput,1994,9(2):137-148.

    [8]KARAMANOS G S,KARNIADAKIS G E.A spectral vanishing viscosity method for large-eddy simulation[J].J comput phys,2000,163:22-50.

    [9]XU C J,PASQUETTI R.Stabilized spectral element computations of high Reynolds number incompressible flows[J].J comput phys,2004,196:680-704.

    [10]KIRBY R M,SHERWIN S J.Stabilization of spectral/hp element methods through spectral vanishing viscosity:Application to fluid mechanics modeling[J].Comput methods appl mech.engrg,2006,195:3128-3144.

    [11]VAN DORSSELAER J L M,HUNDSDORFER W.Stability estimates based on numerical ranges with an application to a spectral method[J].BIT numer math,1994,34(2):228-238.

    (本文責(zé)編:白銀雷)

    O 357.1

    A

    1674-1951(2017)10-0025-05

    2017-08-07;

    2017-09-13

    時(shí)永興(1972—),男,江蘇常州人,工程師,從事火力發(fā)電廠生產(chǎn)管理工作(E-mail:syx_jrhd@126.com)。

    猜你喜歡
    方法
    中醫(yī)特有的急救方法
    中老年保健(2021年9期)2021-08-24 03:52:04
    高中數(shù)學(xué)教學(xué)改革的方法
    化學(xué)反應(yīng)多變幻 “虛擬”方法幫大忙
    變快的方法
    兒童繪本(2020年5期)2020-04-07 17:46:30
    學(xué)習(xí)方法
    可能是方法不對(duì)
    用對(duì)方法才能瘦
    Coco薇(2016年2期)2016-03-22 02:42:52
    最有效的簡(jiǎn)單方法
    山東青年(2016年1期)2016-02-28 14:25:23
    四大方法 教你不再“坐以待病”!
    Coco薇(2015年1期)2015-08-13 02:47:34
    賺錢方法
    国产高潮美女av| 午夜福利在线观看吧| 99久久精品国产国产毛片| 成人性生交大片免费视频hd| 三级毛片av免费| 中文字幕亚洲精品专区| 国产欧美日韩精品一区二区| 亚洲经典国产精华液单| 亚洲乱码一区二区免费版| 久久婷婷人人爽人人干人人爱| 色哟哟·www| 丝袜美腿在线中文| 人妻系列 视频| 99热这里只有精品一区| 成人午夜精彩视频在线观看| 国产av码专区亚洲av| 国产三级中文精品| 亚洲欧洲国产日韩| 欧美性猛交╳xxx乱大交人| videossex国产| 国产成人aa在线观看| 欧美另类亚洲清纯唯美| av在线老鸭窝| 欧美最新免费一区二区三区| 国产三级在线视频| 国产精品日韩av在线免费观看| 欧美成人a在线观看| 狂野欧美白嫩少妇大欣赏| 在线播放国产精品三级| 国产免费又黄又爽又色| 美女脱内裤让男人舔精品视频| 国产三级在线视频| 淫秽高清视频在线观看| 亚洲欧美日韩高清专用| 老师上课跳d突然被开到最大视频| 欧美一级a爱片免费观看看| 欧美激情国产日韩精品一区| 永久网站在线| 亚洲av男天堂| 天天躁日日操中文字幕| 男人狂女人下面高潮的视频| 黄色欧美视频在线观看| 欧美区成人在线视频| 自拍偷自拍亚洲精品老妇| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 午夜福利在线观看吧| 国产乱来视频区| 久久久a久久爽久久v久久| 两性午夜刺激爽爽歪歪视频在线观看| 午夜福利在线在线| 久久人人爽人人爽人人片va| 国产精华一区二区三区| 亚洲不卡免费看| 97超碰精品成人国产| 永久网站在线| 亚洲18禁久久av| 久久6这里有精品| 秋霞在线观看毛片| 亚洲人成网站在线观看播放| 欧美一区二区精品小视频在线| 亚洲综合色惰| 欧美xxxx黑人xx丫x性爽| 亚洲一区高清亚洲精品| 欧美不卡视频在线免费观看| 久久久久久久久大av| 久久午夜福利片| 别揉我奶头 嗯啊视频| 国产黄片美女视频| 97超视频在线观看视频| 国内精品宾馆在线| 国产精品人妻久久久久久| 国产av在哪里看| 免费看美女性在线毛片视频| 内地一区二区视频在线| 成人亚洲欧美一区二区av| 久久热精品热| 激情 狠狠 欧美| 亚洲国产精品成人久久小说| 国产免费又黄又爽又色| 国产爱豆传媒在线观看| 久久久久久国产a免费观看| 国产在视频线精品| 精品99又大又爽又粗少妇毛片| 国产亚洲精品av在线| 纵有疾风起免费观看全集完整版 | 日本黄色视频三级网站网址| 搡女人真爽免费视频火全软件| 综合色丁香网| 成年免费大片在线观看| 欧美激情国产日韩精品一区| 村上凉子中文字幕在线| 最近视频中文字幕2019在线8| 国产精品久久久久久精品电影小说 | 午夜免费激情av| 乱码一卡2卡4卡精品| 欧美激情久久久久久爽电影| 成人国产麻豆网| 观看美女的网站| 亚洲欧美一区二区三区国产| 亚洲成av人片在线播放无| 永久免费av网站大全| 高清日韩中文字幕在线| 国产亚洲av嫩草精品影院| 国产综合懂色| 一个人看的www免费观看视频| 汤姆久久久久久久影院中文字幕 | av.在线天堂| 国产成人91sexporn| 91久久精品国产一区二区成人| 成人鲁丝片一二三区免费| 欧美高清成人免费视频www| 成人毛片a级毛片在线播放| 亚洲成人中文字幕在线播放| 色播亚洲综合网| 亚洲精品456在线播放app| 啦啦啦韩国在线观看视频| 永久网站在线| 国产精品日韩av在线免费观看| 国产真实乱freesex| 国产 一区 欧美 日韩| 国语自产精品视频在线第100页| 精品人妻偷拍中文字幕| 18禁在线播放成人免费| 久久99热6这里只有精品| kizo精华| 3wmmmm亚洲av在线观看| 国产免费视频播放在线视频 | 99久久成人亚洲精品观看| 大香蕉久久网| 国产精品av视频在线免费观看| 久久韩国三级中文字幕| 十八禁国产超污无遮挡网站| 免费观看在线日韩| 国产高清视频在线观看网站| 亚洲精品国产av成人精品| 又爽又黄a免费视频| 伦理电影大哥的女人| 91aial.com中文字幕在线观看| 中文字幕熟女人妻在线| 少妇熟女欧美另类| 免费一级毛片在线播放高清视频| 夜夜看夜夜爽夜夜摸| 婷婷六月久久综合丁香| 国产成人a∨麻豆精品| 99热这里只有是精品在线观看| 可以在线观看毛片的网站| 亚洲成人av在线免费| 日韩精品有码人妻一区| 欧美成人一区二区免费高清观看| 亚洲人成网站高清观看| 男女啪啪激烈高潮av片| 91久久精品国产一区二区三区| 久99久视频精品免费| videossex国产| 老师上课跳d突然被开到最大视频| 国产精品久久久久久久电影| 床上黄色一级片| 日本一二三区视频观看| 久久久欧美国产精品| 国产三级在线视频| 一个人观看的视频www高清免费观看| 精品久久久久久成人av| 国产一级毛片七仙女欲春2| 亚洲自偷自拍三级| 三级毛片av免费| 97超碰精品成人国产| 亚洲欧美精品专区久久| 亚洲精品成人久久久久久| 亚洲精品色激情综合| 1024手机看黄色片| 成人午夜精彩视频在线观看| 久久欧美精品欧美久久欧美| 搡老妇女老女人老熟妇| 老司机福利观看| 啦啦啦韩国在线观看视频| 91午夜精品亚洲一区二区三区| a级毛片免费高清观看在线播放| 国产精品久久久久久精品电影小说 | 久久久a久久爽久久v久久| 寂寞人妻少妇视频99o| 少妇的逼好多水| 久久久国产成人免费| 乱码一卡2卡4卡精品| 国产成人精品久久久久久| 天堂√8在线中文| 一个人看的www免费观看视频| 日韩制服骚丝袜av| 日韩一区二区三区影片| 亚洲av日韩在线播放| 免费看美女性在线毛片视频| 一级毛片aaaaaa免费看小| 美女大奶头视频| 久久久久九九精品影院| 国产成人freesex在线| 日韩大片免费观看网站 | 精品熟女少妇av免费看| 欧美性感艳星| 老女人水多毛片| 久久久午夜欧美精品| 国产成人a∨麻豆精品| 99久久精品一区二区三区| 少妇丰满av| 国产老妇伦熟女老妇高清| 高清视频免费观看一区二区 | 99久久精品国产国产毛片| 欧美日本视频| 亚洲国产欧洲综合997久久,| 成年av动漫网址| 欧美精品一区二区大全| 日日干狠狠操夜夜爽| 成人综合一区亚洲| 看非洲黑人一级黄片| 精品欧美国产一区二区三| 国产成人一区二区在线| 国产黄a三级三级三级人| 蜜桃亚洲精品一区二区三区| 精品一区二区三区人妻视频| 亚洲精品影视一区二区三区av| 美女被艹到高潮喷水动态| 久久精品夜夜夜夜夜久久蜜豆| 国产片特级美女逼逼视频| 亚洲欧美日韩东京热| 国产精品一区二区三区四区久久| 蜜桃久久精品国产亚洲av| 91精品伊人久久大香线蕉| 成人一区二区视频在线观看| 久久久久久久久大av| 麻豆av噜噜一区二区三区| 天堂中文最新版在线下载 | 久久精品久久久久久噜噜老黄 | 国产乱人视频| 美女cb高潮喷水在线观看| 国产亚洲精品av在线| 免费不卡的大黄色大毛片视频在线观看 | 久久精品国产99精品国产亚洲性色| 欧美日韩一区二区视频在线观看视频在线 | 成年版毛片免费区| 欧美一区二区亚洲| 国产伦精品一区二区三区四那| 草草在线视频免费看| 久久久久性生活片| 亚洲精品aⅴ在线观看| 久久6这里有精品| 啦啦啦观看免费观看视频高清| 亚洲真实伦在线观看| 我要看日韩黄色一级片| 嘟嘟电影网在线观看| 高清视频免费观看一区二区 | 男女下面进入的视频免费午夜| 六月丁香七月| 日韩亚洲欧美综合| 国产综合懂色| 久久久欧美国产精品| 边亲边吃奶的免费视频| 男女边吃奶边做爰视频| 国产美女午夜福利| 亚洲精华国产精华液的使用体验| 亚洲欧美精品专区久久| 丝袜美腿在线中文| 天堂影院成人在线观看| 欧美三级亚洲精品| 亚洲国产欧美在线一区| 99在线视频只有这里精品首页| 精品久久久久久久人妻蜜臀av| 久久久午夜欧美精品| 校园人妻丝袜中文字幕| 欧美另类亚洲清纯唯美| 插逼视频在线观看| 久久精品综合一区二区三区| 欧美成人免费av一区二区三区| 国语自产精品视频在线第100页| 黑人高潮一二区| av卡一久久| 深爱激情五月婷婷| 欧美极品一区二区三区四区| 插逼视频在线观看| 高清日韩中文字幕在线| 三级毛片av免费| 天天一区二区日本电影三级| 亚洲伊人久久精品综合 | 久久久国产成人免费| 欧美+日韩+精品| 26uuu在线亚洲综合色| av国产久精品久网站免费入址| 真实男女啪啪啪动态图| 男的添女的下面高潮视频| 国国产精品蜜臀av免费| 亚洲欧美日韩高清专用| 亚洲欧美成人精品一区二区| 中文欧美无线码| 午夜福利在线观看吧| 久久精品熟女亚洲av麻豆精品 | 国产午夜精品论理片| 午夜a级毛片| 99在线视频只有这里精品首页| 尾随美女入室| 亚洲国产精品专区欧美| 国产老妇伦熟女老妇高清| 国产伦理片在线播放av一区| 日本午夜av视频| 亚洲精品国产av成人精品| 国产激情偷乱视频一区二区| 精华霜和精华液先用哪个| 最后的刺客免费高清国语| 成人高潮视频无遮挡免费网站| 亚洲av成人av| 国产色爽女视频免费观看| 免费播放大片免费观看视频在线观看 | 午夜精品一区二区三区免费看| 午夜a级毛片| 国产精品久久久久久久久免| av播播在线观看一区| 午夜激情福利司机影院| 久久久午夜欧美精品| 亚洲欧美日韩高清专用| 亚洲精品aⅴ在线观看| 国产亚洲精品久久久com| 在线观看66精品国产| 视频中文字幕在线观看| 一个人看视频在线观看www免费| 一级毛片电影观看 | av免费在线看不卡| 别揉我奶头 嗯啊视频| eeuss影院久久| 国产69精品久久久久777片| 午夜福利高清视频| 日本一二三区视频观看| 日本-黄色视频高清免费观看| 免费搜索国产男女视频| 伦理电影大哥的女人| 校园人妻丝袜中文字幕| 欧美日韩国产亚洲二区| 淫秽高清视频在线观看| 欧美日本视频| 综合色丁香网| 人妻制服诱惑在线中文字幕| 91久久精品国产一区二区成人| 天天一区二区日本电影三级| 99国产精品一区二区蜜桃av| 欧美激情久久久久久爽电影| 国产成人a∨麻豆精品| 国产毛片a区久久久久| 精品免费久久久久久久清纯| www.av在线官网国产| 午夜亚洲福利在线播放| 看免费成人av毛片| 亚洲欧美日韩卡通动漫| 国产黄片美女视频| 最后的刺客免费高清国语| 一级黄色大片毛片| 一级毛片我不卡| 人人妻人人澡人人爽人人夜夜 | 亚洲精品国产av成人精品| 欧美日韩综合久久久久久| 亚洲aⅴ乱码一区二区在线播放| 久久久久国产网址| av卡一久久| 欧美不卡视频在线免费观看| 美女内射精品一级片tv| 伦理电影大哥的女人| 国产免费视频播放在线视频 | 欧美日韩精品成人综合77777| 桃色一区二区三区在线观看| 久久99蜜桃精品久久| 狂野欧美白嫩少妇大欣赏| 国产高潮美女av| 国产精品久久久久久精品电影小说 | 国产午夜精品久久久久久一区二区三区| 狂野欧美白嫩少妇大欣赏| 日韩av不卡免费在线播放| 久99久视频精品免费| 欧美一区二区精品小视频在线| 国产人妻一区二区三区在| 国产女主播在线喷水免费视频网站 | 最近的中文字幕免费完整| 永久免费av网站大全| 成人av在线播放网站| 午夜福利网站1000一区二区三区| 久久精品夜色国产| 国产亚洲5aaaaa淫片| 麻豆av噜噜一区二区三区| 国产一区二区三区av在线| 国产视频内射| 五月伊人婷婷丁香| 高清在线视频一区二区三区 | 亚洲人成网站在线播| 3wmmmm亚洲av在线观看| 国产av一区在线观看免费| 国产成人91sexporn| 国产大屁股一区二区在线视频| 免费电影在线观看免费观看| 精品国内亚洲2022精品成人| 国产精品日韩av在线免费观看| 赤兔流量卡办理| 淫秽高清视频在线观看| 内射极品少妇av片p| 久久6这里有精品| 人人妻人人澡人人爽人人夜夜 | 精品无人区乱码1区二区| 国产亚洲5aaaaa淫片| 女人被狂操c到高潮| 日本午夜av视频| 99热这里只有精品一区| 一本—道久久a久久精品蜜桃钙片 精品乱码久久久久久99久播 | 毛片女人毛片| 男女啪啪激烈高潮av片| 超碰97精品在线观看| 国产一区二区三区av在线| 中文字幕制服av| 久久婷婷人人爽人人干人人爱| 午夜久久久久精精品| 国产精品,欧美在线| 黄色欧美视频在线观看| 日韩强制内射视频| 搞女人的毛片| 色播亚洲综合网| 国产老妇女一区| 一本一本综合久久| 波多野结衣高清无吗| 久久99热6这里只有精品| av免费观看日本| 亚洲国产精品久久男人天堂| 亚洲一级一片aⅴ在线观看| 成人美女网站在线观看视频| a级毛片免费高清观看在线播放| 尾随美女入室| 国产成人a∨麻豆精品| 精品国产露脸久久av麻豆 | 99热6这里只有精品| 日韩欧美精品v在线| 国产精品久久电影中文字幕| 色吧在线观看| 国产午夜精品一二区理论片| 一夜夜www| 国产黄片美女视频| 国产精品久久视频播放| 看片在线看免费视频| 男女国产视频网站| 人妻制服诱惑在线中文字幕| 一本久久精品| 久久精品人妻少妇| 亚洲欧美精品综合久久99| 亚洲伊人久久精品综合 | www.色视频.com| 亚洲在线观看片| 亚洲精品国产成人久久av| 一个人看视频在线观看www免费| 一夜夜www| 色网站视频免费| 日本欧美国产在线视频| 国产av在哪里看| 午夜福利成人在线免费观看| 久久久久久久久久黄片| 天堂av国产一区二区熟女人妻| 国产视频首页在线观看| 少妇高潮的动态图| 国产精华一区二区三区| 午夜福利在线在线| 内地一区二区视频在线| АⅤ资源中文在线天堂| 欧美激情国产日韩精品一区| 性色avwww在线观看| 男的添女的下面高潮视频| 桃色一区二区三区在线观看| 菩萨蛮人人尽说江南好唐韦庄 | 日韩精品有码人妻一区| 黄色欧美视频在线观看| 免费大片18禁| 小说图片视频综合网站| 国产人妻一区二区三区在| 国产精品无大码| av线在线观看网站| 国产av不卡久久| 在线播放无遮挡| 99久久无色码亚洲精品果冻| 边亲边吃奶的免费视频| 女人被狂操c到高潮| www.色视频.com| 国产精品电影一区二区三区| 日本免费a在线| 国产精品1区2区在线观看.| 在现免费观看毛片| 亚洲av中文av极速乱| 一卡2卡三卡四卡精品乱码亚洲| 久久婷婷人人爽人人干人人爱| 联通29元200g的流量卡| 欧美一区二区亚洲| 夫妻性生交免费视频一级片| 亚洲人与动物交配视频| av国产免费在线观看| 高清在线视频一区二区三区 | av黄色大香蕉| 国产视频首页在线观看| 看十八女毛片水多多多| 午夜爱爱视频在线播放| av国产免费在线观看| 国产一区有黄有色的免费视频 | 性插视频无遮挡在线免费观看| 精品熟女少妇av免费看| 欧美丝袜亚洲另类| 五月伊人婷婷丁香| 国产色爽女视频免费观看| 乱人视频在线观看| 免费无遮挡裸体视频| 国产精品嫩草影院av在线观看| 国产精品福利在线免费观看| 中文在线观看免费www的网站| 久久久久精品久久久久真实原创| 成人一区二区视频在线观看| 成人午夜高清在线视频| 亚洲欧美精品自产自拍| 美女高潮的动态| 精品久久久久久电影网 | 日日摸夜夜添夜夜爱| 国内精品美女久久久久久| 干丝袜人妻中文字幕| 国产精品国产三级国产av玫瑰| 亚洲成av人片在线播放无| 亚洲高清免费不卡视频| 国内少妇人妻偷人精品xxx网站| 成人二区视频| 亚洲国产欧洲综合997久久,| 伦精品一区二区三区| 日韩欧美精品免费久久| 婷婷色麻豆天堂久久 | 精品久久国产蜜桃| 伦理电影大哥的女人| 人妻少妇偷人精品九色| 久久精品久久久久久噜噜老黄 | 一区二区三区四区激情视频| 国内精品美女久久久久久| 嘟嘟电影网在线观看| 日韩欧美三级三区| 免费av观看视频| 日韩欧美精品免费久久| 自拍偷自拍亚洲精品老妇| 国产亚洲av片在线观看秒播厂 | 最后的刺客免费高清国语| 汤姆久久久久久久影院中文字幕 | 亚洲人成网站在线观看播放| 性插视频无遮挡在线免费观看| 日韩一本色道免费dvd| 精品99又大又爽又粗少妇毛片| 国产精品三级大全| 三级经典国产精品| 秋霞在线观看毛片| 我的女老师完整版在线观看| 欧美+日韩+精品| 亚洲av免费高清在线观看| 亚洲成人中文字幕在线播放| 亚洲丝袜综合中文字幕| 性插视频无遮挡在线免费观看| 22中文网久久字幕| 国产在线一区二区三区精 | 久久久精品大字幕| 狠狠狠狠99中文字幕| 成人毛片60女人毛片免费| 99久国产av精品| 国产男人的电影天堂91| 精品酒店卫生间| 午夜激情欧美在线| 我的女老师完整版在线观看| 亚洲在久久综合| 麻豆成人午夜福利视频| 亚洲美女视频黄频| 久99久视频精品免费| 亚洲欧洲国产日韩| 18+在线观看网站| 亚洲高清免费不卡视频| 天天躁夜夜躁狠狠久久av| 国产精品久久电影中文字幕| 狂野欧美激情性xxxx在线观看| 在线a可以看的网站| 韩国高清视频一区二区三区| 亚洲av福利一区| 国产精品.久久久| 99热全是精品| 国产精品人妻久久久影院| 伊人久久精品亚洲午夜| 国产精品美女特级片免费视频播放器| 色综合亚洲欧美另类图片| 91精品一卡2卡3卡4卡| 毛片女人毛片| 91aial.com中文字幕在线观看| 午夜精品一区二区三区免费看| 一二三四中文在线观看免费高清| 久久欧美精品欧美久久欧美| 国产高清三级在线| 69人妻影院| 亚洲激情五月婷婷啪啪| 久久久a久久爽久久v久久| 国产高清有码在线观看视频| 亚洲精品亚洲一区二区| 国产一区有黄有色的免费视频 | 一区二区三区乱码不卡18| 国产三级在线视频| 青春草视频在线免费观看| 国产免费视频播放在线视频 | 午夜亚洲福利在线播放| 午夜福利高清视频| 老女人水多毛片| 69av精品久久久久久| 午夜福利高清视频| 中文在线观看免费www的网站| 91在线精品国自产拍蜜月| 欧美日韩一区二区视频在线观看视频在线 | 久久久久久久久大av| 日韩国内少妇激情av| 成人性生交大片免费视频hd| 国产成人免费观看mmmm| 黄色配什么色好看| 热99re8久久精品国产| 99久久中文字幕三级久久日本| 亚洲国产精品合色在线| 麻豆精品久久久久久蜜桃|