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

    應(yīng)用重心插值的Burgers方程數(shù)值解法

    2025-01-21 00:00:00滕宇航賴藝穎黃浪揚

    摘要: 利用重心插值構(gòu)造求解Burgers方程的無網(wǎng)格數(shù)值方法。在時間方向用Crank-Nicolson差分法對方程進行離散,在空間方向用重心插值切比雪夫配點法逼近函數(shù)本身及其空間導(dǎo)數(shù)。 對全離散數(shù)值格式進行相容性分析。數(shù)值算例表明:與經(jīng)典的有限差分方法比較,重心插值配點法用較少的節(jié)點就能達到較高的精度。

    關(guān)鍵詞: Burgers方程;重心插值配點法;Crank-Nicolson差分法;相容性分析

    中圖分類號: O 241.82文獻標志碼: A 文章編號: 1000-5013(2025)01-0104-09

    Numerical Solution Method of Burgers Equation Using Barycentric Interpolation

    TENG Yuhang,LAI Yiying,HUANG Langyang

    (School of Mathematical Science,Huaqiao University,Quanzhou 362021,China)

    Abstract: The meshless numerical method of Burgers equation is solved using barycentric interpolation construction. The equation is discretized using the Crank-Nicolson difference method in the time direction. The function is approximated to itself and so is its spatial derivative using barycentric interpolation Chebyshev collocation method in spatial direction. The compatibility analysis of the fully discrete numerical value scheme is perfomed. Numerical experiments show that,compared with the classical finite difference method,the barycentric interpolation collocation method can achieve higher accuracy with fewer nodes.

    Keywords: Burgers equation;barycentric interpolation collocation method;Crank-Nicolson difference method;compatibility analysis

    Burgers方程是一類描述對流與擴散相互作用的非線性偏微分方程1,是Navier-Stokes方程中忽略壓力項后的簡化形式??紤]如下Burgers方程,即

    式(1)中:Ω Rd(d=1,2)為有界凸區(qū)域,邊界分段光滑。

    Burgers方程在流體力學(xué)、氣體動力學(xué)、交通流等眾多領(lǐng)域都有著重要的作用和地位。許多學(xué)者對該問題做了大量研究,并提出了數(shù)值解法。陳蓮2建立兩種求解Burgers方程的Crank-Nicolson差分法格式,并對其進行誤差分析。Wang等3建立一種求解粘性Burgers方程的能量守恒型四階隱式緊差分格式。Zhang等4利用線性化緊致差分格式求解具有Burgers型非線性項的二維Sobolev方程。Xu等5建立兩種求解非線性Burgers方程的修正的非協(xié)調(diào)有限元全離散格式,并對其進行超收斂分析。Wang等6建立一種求解Burgers方程的多區(qū)域Galerkin方法,并分析其收斂性及穩(wěn)定性。Zhao等7利用時空連續(xù)Galerkin方法對二維Burgers方程進行數(shù)值求解,并給出先驗誤差估計。Wang等8利用弱伽遼金有限元方法求解一類時間分數(shù)階的廣義Burgers方程。

    以上求解Burgers方程的數(shù)值方法都是基于網(wǎng)格剖分建立離散格式,Chebyshev節(jié)點的重心Lagrange插值公式可以有效克服Runge現(xiàn)象。很多學(xué)者將其推廣到求解各類微分方程與積分方程,如Volterra積分方程9、Allen-Cahn方程10-13、Black-Scholes方程14、分數(shù)階微分方程15-16等。?mer17應(yīng)用重心插值配點法數(shù)值求解了二維和三維Klein-Gordon-Schr?dinger(KGS)方程?;诖耍疚膶?yīng)用重心插值的Burgers方程數(shù)值解法進行研究。

    1 重心Lagrange插值

    1.1 重心Lagrange插值函數(shù)

    設(shè)插值節(jié)點為xj,對應(yīng)的實數(shù)為yj(j=0,1,…,n),使用多項式插值,則在次數(shù)不超過n的多項式空間中,存在唯一的重心Lagrange插值函數(shù)p(x),滿足p(xj)=yj,j=0,1,…,n。重心Lagrange插值函數(shù)為

    p(x)=∑n/j=0wj/(x-xj)yj/∑n/i=0wi/(x-xi)=∑n/j=0Dj(x)yj。(2)

    式(2)中:Dj(x)為重心Lagrange插值基函數(shù);wj為重心Lagrange插值函數(shù)的重心權(quán),wj

    wj=∏n/i=0,i≠j(xj-xi-1。(3)

    由于重心Lagrange插值函數(shù)的重心權(quán)只與插值節(jié)點分布有關(guān),在利用重心插值逼近未知函數(shù)時,采用第二類 Chebyshev節(jié)點,即

    xj=cosj/nπ, j=0,1,…,n。(4)

    可以得到重心Lagrange插值函數(shù)的重心權(quán)為

    wj=(-1)δj,δj=1/2, j=0或n,1, 其他。(5)

    對式(2)進行求導(dǎo),得到函數(shù)p(x)在插值節(jié)點xi處的導(dǎo)數(shù),ps(xi)為

    ps(xi)=dsp(xi)/dxs=∑n/j=0D(s)j(xi)yj=∑n/j=0D(s)ijyj。(6)

    將其寫成矩陣形式即為

    p(s)=D(s)p。(7)

    式(7)中:D(s)為關(guān)于節(jié)點xi(i=0,1,…,n)的s階微分矩陣,D(s)=(D(s)i,jn×n,s=1,2。

    通過對重心Lagrange插值基函數(shù)求導(dǎo),分別得到一、二階微分矩陣的元素11

    D(1)i,jji/xi-xj, i≠j,D(1)i,i=-∑n/j=0,j≠iD(1)ij,D(2)i,j=2(D(1)iiD(1)i,j-D(1)i,j/xi-xj), i≠j,D(2)i,i=-∑n/j=0,j≠iD(2)i,j。(8)

    1.2 逼近性質(zhì)

    設(shè)p(x,y)是u(x,y)的重心Lagrange插值函數(shù),滿足p(xm,yn)=u(xm,yn),則定義誤差函數(shù)為

    e(x,y)=u(x,y)-p(x,y)。(9)

    引理118 如果u(x,y)∈Cm+1([a,b]×[c,d]),其中,m=max{m,n},則

    |e(x,y)|≤u(m+1)(C1elx/2mm+C2ely/2nn)。(10)

    式(10)中:e為自然對數(shù);lx代表區(qū)間[a,b]的一半;ly代表區(qū)間[c,d]的一半。

    同理可以得到

    ex(x,y)≤C*1‖?(m+1)xu‖elx/2(m-1)m-1+C2‖?(n+1)yu‖ely/2nn,exx(x,y)≤C**1‖?(m+1)xu‖elx/2(m-2)m-2+C2‖?(n+1)yu‖ely/2nn,ey(x,y)≤C1‖?(m+1)xu‖elx/2mm+C*2‖?(n+1)yu‖ely/2(n-1)n-1,ey,y(x,y)≤C1‖?(m+1)xu‖elx/2mm+C**2‖?(n+1)yu‖ely/2(n-2)n-2。(11)

    式(11)中:C1,C*1,C**1,C2,C*2,C**2為常數(shù)。

    2 二維Burgers方程的離散格式

    考慮二維Burgers方程,對時間區(qū)域進行K等分,記步長為τ=T/K,T為時間計算的終點??臻g區(qū)域采用第二類Chebyshev節(jié)點分別離散為m+1個計算節(jié)點a=x0lt;x1lt;…lt;xm=b和n+1個計算節(jié)點 c=y0lt;y1lt;…lt;yn=d,得到計算節(jié)點(xi,yj,tk),0≤i≤m,0≤j≤n,0≤k≤K。

    2.1 半離散格式及相容性分析

    記函數(shù)u(x,y,t)在節(jié)點x0,x1,…,xm上的值為u(xi,y,t)=ui(y,t),則函數(shù)u(x,y,t)在節(jié)點x0,x1,…,xm上的重心插值函數(shù)為

    u(x,y,t)=∑m/s=0Ls(x)us(y,t)。(12)

    式(12)中:Ls(x)表示x方向的重心插值基函數(shù)。

    將式(12)代入方程(1)中,并令方程在點x0,x1,…,xm上精確成立,得到常微分方程組,即

    ∑m/s=0Lsxi?us(y,t)/?t+/∑m/s=0Lsxius(y,t)∑m/s=0L′sxius(y,t)+∑m/s=0Lsxi?2us(y,t)/?y2=/

    v∑m/s=0L″sxius(y,t)+∑m/s=0Lsxi?2us(y,t)/?y2。(13)

    式(13)中:L′s(xi)=C(1)s,i為關(guān)于節(jié)點x0,x1,…,xm的一階微分矩陣C(1)中的元素,類似地,L″s(xi)=C(2)s,i為二階微分矩陣C(2)中的元素。

    記uij(t)=ui(yj,t),函數(shù)ui(y,t)在節(jié)點y0,y1,…,yn的重心插值函數(shù)為

    ui(y,t)=∑n/r=0Tr(y)ui,r(t)。(14)

    式(14)中:Tr(y)表示y方向的重心插值基函數(shù)。

    將式(14)代入式(13),微分矩陣形式為

    dU(t)/dt+diag(U(t))(C(1)In+1+Im+1D(1))U(t)=v(C(2)In+1+Im+1D(2))U(t)。(15)

    式(15)中:

    U(t)=[U0(t)…Um(t)]T;

    Ui(t)=[ui,0(t)…uin(t)]T(i=0,…,m);

    D(k)為關(guān)于節(jié)點y0,y1,…,yn的k階微分矩陣。

    為了方便進行相容性分析,算子定義為

    Γu(x,y,t):=ut+u(ux+uy)-v(uxx+uy,y)。(16)

    令u(xm,yn,t)為u(x,y,t)的數(shù)值解,則可以得到

    Γu(xm,yn,t)=0。(17)

    定理1 若u(x,y,t)∈C(m+1)([a,b]×[c,d])×C(0)(0,T],其中,m=max{m,n},假設(shè)u(xm,yn,t)有界,則誤差估計為

    u(x,y,t)-u(xm,yn,t)≤C**1?(m+1)xuelx/2(m-2)m-2+/

    C**2‖?(n+1)yu‖ely/2(n-2)n-2。(18)

    證明:Γu(x,y,t)-Γu(xm,yn,t)=ut(x,y,t)-ut(xm,yn,t)+/

    u(x,y,t)ux(x,y,t)-u(xm,yn,t)ux(xm,yn,t)+/

    u(x,y,t)uy(x,y,t)-u(xm,yn,t)uy(xm,yn,t)+/

    v(uxx(xm,yn,t)-(x,y,t))+v(uyy(xm,yn,t)-uyy(x,y,t))=/R1+R2+R3+R4+R5。(19)

    式(19)中:

    R1=ut(x,y,t)-ut(xm,yn,t)=ut(x,y,t)-ut(xm,y,t)+

    ut(xm,y,t)-ut(xm,yn,t)=et(xm,y,t)+et(xm,yn,t)。(20)

    由引理1,有R1=/et(xm,y,t)+et(xm,yn,t)≤/C1‖?(m+1)xu‖elx/2mm+C2‖?(n+1)yu‖ely/2nn。(21)

    R2的估計為

    R2=/u(x,y,t)ux(x,y,t)-u(xm,yn,t)ux(xm,yn,t)≤/ux(x,y,t)u(x,y,t)-u(xm,yn,t)+/

    u(xm,yn,t)ux(x,y,t)-ux(xm,yn,t)≤/

    C*1‖?(m+1)xu‖elx/2(m-1)m-1+C2?(n+1)yuely/2nn。(22) 同理有

    R3=u(x,y,t)uy(x,y,t)-u(xm,yn,t)uy(xm,yn,t)≤

    C1‖?(m+1)xu‖elx/2mm+C*2‖?(n+1)yu‖ely/2(n-1)n-1。 (23)

    類似地,對于R4,R5的估計分別為

    R4=vex,x(xm,y,t)+vex,x(xm,yn,t)C≤

    C**1‖?(m+1)xu‖(elx/2(m-2))m-2+C2‖?(n+1)yu‖-ely/2nn。(24)

    R5=vey,y(xm,y,t)+veyy(xm,yn,t)≤

    C1‖?(m+1)xu‖elx/2mm+C**2‖?(n+1)yu‖ely/2(n-2)n-2。 (25)

    將式(21)~(25)代入式(19),則定理得證。

    2.2 全離散格式及相容性分析

    定理2 設(shè)u(x,y,t)∈C(m+1)([a,b]×[c,d])×C(2)(0,T],假設(shè)u(xm,yn,t)有界,使用重心插值配點法求解的數(shù)值解為uh(xi,yj,tk),則如下式子成立,即

    u(x,y,t)-uh(xi,yj,tk)≤Cτ2+‖u(m+1)elx/2(m-2)m-2+ely/2(n-2)n-2。(26)

    證明:對方程時間方向采用Crank-Nicolson差分法離散得到的數(shù)值解為uk+1=u(x,y,tk+1),則

    δuk+1/2t+uk+1/2(uk+1/2x+uk+1/2y)-v(uk+1/2x,x+uk+1/2y,y)=Rk。(27)

    式(27)中:δuk+1/2t=1/τ(u(x,y,tk+1)-u(x,y,tk)),Rk=δuk+1/2t-uk+1/2t是時間方向的截斷誤差。

    由Taylor展開,可得

    Rk≤Cτ2。(28)

    空間方向利用重心Lagrange插值配點法離散,有/δuht(xi,yj,tk)+u(xi,yj,tk+1/2)(ux(xi,yj,tk+1/2)+uy(xi,yj,tk+1/2))-/

    v(uxx(xi,yj,tk+1/2)+uyy(xi,yj,tk+1/2))=Rkij。(29)

    式(29)中:εij為空間截斷誤差。

    由式(27),(29),可以得到

    /δuk+1/2t-δuht(xi,yj,tk)+uk+1/2(uk+1/2x+uk+1/2y)-/

    u(xi,yj,tk+1/2)(ux(xi,yj,tk+1/2)+uy(xi,yj,tk+1/2))-/

    v(uk+1/2x,x+uk+1/2y,y)+v(uxx(xi,yj,tk+1/2)+uy,y(xi,yj,tk+1/2))=-εi,j。(30)

    類似定理2的推導(dǎo),有

    εi,j≤C‖u(m+1)elx/2(m-2)m-2+ely/2(n-2)n-2。(31)

    結(jié)合式(28),(31),則定理得證。

    記uk(x,y)=u(x,y,tk),k=0,1,…,K,采用Crank-Nicolson差分法對時間離散,對方程中的非線性項線性化處理,得到時間半離散格式,即

    uk+1-uk/τ+uk(uk+1x+uk+1y)+uk+1(ukx+uky)/2=v(uk+1xx+uk+1y,y+ukxx+uky,y)/2。(32)

    將微分矩陣代入時間半離散格式,有Burgers方程全離散格式,即

    I+τ/2(diag(Uk)A*+diag(A*Uk))-vτ/2B*Uk+1=(I+vτ/2B*)Uk。(33)

    3 數(shù)值算例

    為驗證重心插值配點法求解Burgers方程的有效性和合理性,通過算例對重心插值配點法穩(wěn)定性進行驗證,與Crank-Nicolson差分法對比說明重心插值配點法的高精度。記Ue,Uh分別為精確解及數(shù)值解,定義在無窮范數(shù)意義下的絕對誤差和相對誤差分別為

    E=Ue-Uh, Er=Ue-Uh/Ue。

    上式中:·表示無窮范數(shù)。

    算例1 考慮d=1時的一維Burgers方程(1),取T=1,a=0,b=1,初始條件為

    當(dāng)v=0.01,固定空間節(jié)點n=16時,改變時間步長,可得到重心插值配點法的時間收斂階,如表1所示。由表1可知:隨著K的增加,利用重心插值配點法求出的方程數(shù)值解與精確解的絕對誤差的數(shù)量級從10-3減到 10-5,因此,時間收斂階為二階精度。

    當(dāng)固定時間步長τ=0.001時,改變空間節(jié)點數(shù)量,可得到重心插值配點法與有限差分法的誤差,如表2所示。由表2可知: 重心插值配點法在空間選取10個節(jié)點時,最大絕對誤差即可達到10-6量級,而有限差分法空間選取128個節(jié)點時,最大絕對誤差才能達到10-6量級,這表明重心插值配點法在空間上用較少的點可達到更高的精度。

    重心插值配點法數(shù)值解與精確解(τ=0.05,n=20),如圖1所示。由圖 1 可知:兩種解吻合。

    絕對誤差(τ=0.05,n=20),如圖2所示??臻g收斂階(τ=0.001,v=0.01),如圖3所示。由圖3可知:重心插值配點法具有指數(shù)收斂速度。不同參數(shù)下重心插值配點法的數(shù)值解,如圖4所示。

    圖4(a)為最終時間T=1,2,4,8,16的數(shù)值解圖像,圖4(b)為不同的雷諾數(shù)Re=10,20,50,100的數(shù)值解圖像。由圖4可知:隨著T,Re的增大,數(shù)值解保持穩(wěn)定,說明了數(shù)值法的穩(wěn)定性。

    重心插值配點法的數(shù)值解與精確解的等勢圖(v=0.1),如圖5所示。

    算例2 考慮如下二維Burgers方程,即

    ut+u?u/?x+?u/?y=v?2u/?x2+?2u/?y2,(x,y)∈[0,1]×[0,1],t∈[0,1]。

    初值條件和邊值條件由精確解確定,方程的精確解為

    重心插值配點法時間收斂階,如表3所示。

    當(dāng)v=0.1,空間節(jié)點數(shù)m=n,固定空間節(jié)點數(shù)n=16時,改變時間節(jié)點數(shù),重心插值配點法與有限差分法的誤差(τ=0.001),如表4所示。由表4可知:當(dāng)n取8時,重心插值配點法的誤差具有10-4量級,而有限差分法的誤差只有10-3量級;重心插值配點法n=12時誤差有10-5量級,而有限差分法在n=64時誤差只有10-5量級。

    重心插值配點法的數(shù)值解與精確解(τ=0.001,n=32),如圖6所示。由圖6可知:數(shù)值解與精確解是相符的。

    絕對誤差(τ=0.001,n=32),如圖7所示。由圖7可知:數(shù)值解與精確解是相符的。

    空間收斂階(τ=0.001),如圖8所示。由圖8可知:重心插值配點法在求解二維問題時,空間具有指數(shù)收斂。

    算例2結(jié)果表明:重心插值配點法在二維問題中,具有高精度性與有效性,且可通過更少的節(jié)點得到相對高的精度,與理論分析相符。

    5 結(jié)束語

    將Crank-Nicolson差分法與重心插值切比雪夫配點法相結(jié)合,提出了一種求解Burgers的有效數(shù)值法,對方程的非線性項進行線性化處理,避免迭代帶來大計算量,并給出了全離散法的相容性分析。通過數(shù)值算例驗證數(shù)值法的有效性及高精度,通過與Crank-Nicolson差分法進行比較,重心插值法可以用較少的節(jié)點得到較高的精度。

    參考文獻:

    [1] BURGERS J M.A mathematical model illustrating the theory of turbulence[J].Advances in Applied Mechanics,1948,1:171-199.DOI:10.1016/S0065-2156(08)70100-5.

    [2] 陳蓮.一維Burgers方程的幾種有限差分解法[D].南充:西華師范大學(xué),2020.

    [3] WANG Xuping,ZHANG Qifeng,SUN Zhizhong.The pointwise error estimates of two energy-preserving fourth-order compact schemes for viscous Burgers equation[J].Advances in Com-putational Mathematics,2021,47(2):23.DOI:10.1007/S10444-021-09848-9.

    [4] ZHANG Qifeng,QIN Yifan,SUN Zhizhong.Linearly compact scheme for 2D Sobolev equation with Burgers′ type nonlinearity[J].Numerical Algorithms,2022,91(3):1081-1114.DOI:10.1007/S11075-022-01293-Z.

    [5] XU Chao,PEU Lifang.Unconditional superconvergence analysis of two modified finite element fully discrete schemes for nonlinear Burgers′ equation[J].Applied Numerical Mathematics,2023,185:1-17.DOI:10.1016/j.apnum.2022.11.008.

    [6] WANG Chuan,WANG Tianjun.A multi-domain Galerkin method with numerical integration for the Burgers equation[J].International Journal of Computer Mathematics,2023,100(5):927-947.DOI:10.1080/00207160.2023.2171265.

    [7] ZHAO Zhihui,LI Hong.Numerical study of two-dimensional Burgers equation by using a continuous Galerkin method[J].Computers and Mathematics with Applications,2023,149(1):38-48.DOI:10.1016/J.CAMWA.2023.08.030.

    [8] WANG Haifeng,XU Da,ZHOU Jun,et al.Weak Galerkin finite element method for a class of time fractional generalized Burgers equation[J].Numerical Methods for Partial Differential Equations,2021,37(1):732-749.DOI:10.1002/num.22549.

    [9] 于孟文,張新東.重心插值配點法求解Volterra積分方程[J].新疆師范大學(xué)學(xué)報(自然科學(xué)版),2023,42(1):75-80.DOI:10.14100/j.cnki.1008-9659.2023.01.010.

    [10] DENG Yangfang,WENG Zhifeng.Barycentric interpolation collocation method based on Crank-Nicolson scheme for the Allen-Cahn equation[J].AIMS Mathematics,2021,6(4):3857-3873.DOI:10.3934/MATH.2021229.

    [11] 黃蓉,鄧楊芳,翁智峰.SAV/重心插值配點法求解Allen-Cahn方程[J].應(yīng)用數(shù)學(xué)和力學(xué),2023,44(5):573-582.DOI:10.21656/1000-0887.430149.

    [12] 鄧楊芳,姚澤豐,汪精英,等.二維Allen-Cahn方程的有限差分法/配點法求解[J].華僑大學(xué)學(xué)報(自然科學(xué)版),2020,41(5):690-694.DOI:10.11830/ISSN.1000-5013.202001001.

    [13] 翁智峰,姚澤豐,賴淑琴.重心插值配點法求解Allen-Cahn方程[J].華僑大學(xué)學(xué)報(自然科學(xué)版),2019,40(1):133-140.DOI:10.11830/ISSN.1000-5013.201806043.

    [14] 賴舒琴,華之維,翁智峰.重心插值配點法求解Black-Scholes方程[J].聊城大學(xué)學(xué)報(自然科學(xué)版),2020,33(5):1-8.DOI:10.19728/j.issn1672-6634.2020.05.001.

    [15] 黃蓉,翁智峰.時間分數(shù)階Allen-Cahn方程的重心插值配點法[J].華僑大學(xué)學(xué)報(自然科學(xué)版),2022,43(4):553-560.DOI:10.11830/ISSN.1000-5013.202104060.

    [16] LI Jin,SU Xiaoning,ZHAO Kaiyan.Barycentric interpolation collocation algorithm to solve fractional differential equations[J].Mathematics and Computers in Simulation,2023,205:340-367.DOI:10.1016/J.MATCOM.2022.10.005.

    [17] ?MER O.Application of a collocation method based on linear barycentric interpolation for solving 2D and 3D Klein-Gordon-Schr?dinger (KGS) equations numerically[J].Enginee-Ring Computations,2021,38(5):2394-2414.DOI:10.1108/EC-06-2020-0312.

    [18] YI Shichao,YAO Linquan.A steady barycentric Lagrange interpolation method for the 2D higher order time fractional teleglaph equation with nonlocal boundary condition with error analysis[J].Numerical Methods for Partial Differential Equations,2019,35(5):1694-1716.DOI:10.1002/num.22371.

    (責(zé)任編輯:陳志賢 英文審校:黃心中)

    好男人电影高清在线观看| 精品一区二区三区视频在线观看免费| 久久精品国产亚洲av涩爱 | 九九在线视频观看精品| 一个人看视频在线观看www免费 | 免费高清视频大片| 久久中文看片网| 精品乱码久久久久久99久播| eeuss影院久久| 91av网一区二区| 色在线成人网| 亚洲美女黄片视频| 日韩成人在线观看一区二区三区| 亚洲中文字幕一区二区三区有码在线看| 又紧又爽又黄一区二区| 日本成人三级电影网站| 色老头精品视频在线观看| 欧美绝顶高潮抽搐喷水| av专区在线播放| 久久久精品大字幕| 一级毛片女人18水好多| 日韩欧美在线二视频| 51国产日韩欧美| 18禁黄网站禁片午夜丰满| 18禁黄网站禁片午夜丰满| 99久久精品国产亚洲精品| 最新在线观看一区二区三区| 99热6这里只有精品| 最近最新免费中文字幕在线| 欧美zozozo另类| 少妇裸体淫交视频免费看高清| 噜噜噜噜噜久久久久久91| 哪里可以看免费的av片| 热99在线观看视频| 亚洲自拍偷在线| 欧美又色又爽又黄视频| 亚洲 国产 在线| 久久香蕉精品热| 国产午夜精品论理片| 日日摸夜夜添夜夜添小说| 色综合婷婷激情| 无限看片的www在线观看| 在线观看舔阴道视频| 欧美+日韩+精品| 欧美乱码精品一区二区三区| 亚洲成av人片在线播放无| 精品国内亚洲2022精品成人| 啦啦啦韩国在线观看视频| 日韩大尺度精品在线看网址| a级一级毛片免费在线观看| 久久久久久久久中文| 成人精品一区二区免费| 天堂av国产一区二区熟女人妻| 丁香六月欧美| 9191精品国产免费久久| 搡老熟女国产l中国老女人| 欧美乱色亚洲激情| 少妇熟女aⅴ在线视频| 久久精品国产亚洲av涩爱 | 色播亚洲综合网| 中文字幕熟女人妻在线| 偷拍熟女少妇极品色| 国产一区二区三区在线臀色熟女| 国产精品自产拍在线观看55亚洲| 欧美日韩乱码在线| 亚洲av一区综合| 国产一区二区亚洲精品在线观看| 国产精品一区二区三区四区久久| 少妇的逼好多水| 嫩草影院精品99| 中文字幕熟女人妻在线| 久久久久久久精品吃奶| 一边摸一边抽搐一进一小说| 亚洲欧美日韩高清专用| 成人性生交大片免费视频hd| 女警被强在线播放| 女警被强在线播放| 国产97色在线日韩免费| 婷婷丁香在线五月| 欧美中文综合在线视频| 日韩欧美国产在线观看| 尤物成人国产欧美一区二区三区| 中文在线观看免费www的网站| 制服人妻中文乱码| 美女被艹到高潮喷水动态| 免费一级毛片在线播放高清视频| 国产亚洲精品av在线| 国产精品一区二区三区四区免费观看 | 岛国在线免费视频观看| 欧美成人性av电影在线观看| 一区二区三区免费毛片| 日本 欧美在线| 国产欧美日韩一区二区精品| 18禁黄网站禁片午夜丰满| www国产在线视频色| 一本精品99久久精品77| 少妇的逼水好多| 国产精品亚洲美女久久久| 欧美日韩精品网址| 亚洲 欧美 日韩 在线 免费| 悠悠久久av| 少妇人妻一区二区三区视频| 国产精品久久久人人做人人爽| 男女床上黄色一级片免费看| 一个人免费在线观看的高清视频| 一夜夜www| 久久九九热精品免费| 午夜激情福利司机影院| 国产亚洲欧美98| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 狠狠狠狠99中文字幕| 一个人观看的视频www高清免费观看| 一区二区三区免费毛片| 好男人在线观看高清免费视频| 国产视频内射| 久久久国产成人精品二区| 日韩欧美国产在线观看| 一区二区三区国产精品乱码| 久久精品亚洲精品国产色婷小说| 精品不卡国产一区二区三区| 国产精品久久久久久亚洲av鲁大| 91字幕亚洲| 久久久久性生活片| 精品久久久久久久末码| 欧美乱妇无乱码| 人妻丰满熟妇av一区二区三区| 俺也久久电影网| 最近最新免费中文字幕在线| 给我免费播放毛片高清在线观看| 一本一本综合久久| 又爽又黄无遮挡网站| 日韩欧美在线二视频| 欧美一级a爱片免费观看看| 精品久久久久久久人妻蜜臀av| 亚洲成av人片在线播放无| 亚洲中文字幕日韩| 日本与韩国留学比较| 九色成人免费人妻av| 在线观看午夜福利视频| svipshipincom国产片| 网址你懂的国产日韩在线| 亚洲成人久久性| 欧美av亚洲av综合av国产av| eeuss影院久久| 一二三四社区在线视频社区8| 国产黄色小视频在线观看| 嫩草影院入口| 午夜激情福利司机影院| 久久伊人香网站| 国产探花在线观看一区二区| 少妇熟女aⅴ在线视频| 亚洲成人久久性| 亚洲国产精品sss在线观看| 丰满人妻一区二区三区视频av | 18美女黄网站色大片免费观看| 可以在线观看的亚洲视频| aaaaa片日本免费| 法律面前人人平等表现在哪些方面| 国产高清三级在线| 窝窝影院91人妻| 少妇高潮的动态图| 日韩欧美国产一区二区入口| av欧美777| 九色成人免费人妻av| a级毛片a级免费在线| 欧美zozozo另类| 成人av在线播放网站| 日本在线视频免费播放| 天堂av国产一区二区熟女人妻| 亚洲va日本ⅴa欧美va伊人久久| 国产精品久久久久久久电影 | 黄片大片在线免费观看| bbb黄色大片| 成年女人毛片免费观看观看9| 国产在视频线在精品| 成人永久免费在线观看视频| 又黄又爽又免费观看的视频| 亚洲中文字幕一区二区三区有码在线看| 一区二区三区高清视频在线| 女人高潮潮喷娇喘18禁视频| 免费电影在线观看免费观看| 18禁在线播放成人免费| 亚洲av免费高清在线观看| 国产单亲对白刺激| 日日摸夜夜添夜夜添小说| h日本视频在线播放| 精品久久久久久,| 国产精品日韩av在线免费观看| 禁无遮挡网站| 99久久综合精品五月天人人| 国产av麻豆久久久久久久| 国产老妇女一区| 久久国产乱子伦精品免费另类| 黄片小视频在线播放| 成人午夜高清在线视频| 狂野欧美白嫩少妇大欣赏| 级片在线观看| www日本在线高清视频| 成人av一区二区三区在线看| 色综合亚洲欧美另类图片| 久久久精品欧美日韩精品| 中国美女看黄片| 又爽又黄无遮挡网站| 亚洲av中文字字幕乱码综合| 老汉色∧v一级毛片| 亚洲精品成人久久久久久| 长腿黑丝高跟| 色精品久久人妻99蜜桃| 热99re8久久精品国产| 国产高清videossex| 欧美+日韩+精品| 国产高清视频在线播放一区| eeuss影院久久| 国产精品三级大全| 国产精华一区二区三区| 亚洲第一欧美日韩一区二区三区| a级一级毛片免费在线观看| 一区二区三区激情视频| 色噜噜av男人的天堂激情| 露出奶头的视频| 在线观看免费午夜福利视频| 中文字幕人妻熟人妻熟丝袜美 | 淫秽高清视频在线观看| 久久性视频一级片| 欧美中文综合在线视频| 老司机深夜福利视频在线观看| av在线蜜桃| 国产蜜桃级精品一区二区三区| 欧美性感艳星| 亚洲欧美日韩高清专用| 国产久久久一区二区三区| 在线播放无遮挡| 99视频精品全部免费 在线| 国产真人三级小视频在线观看| 亚洲国产欧美人成| 18+在线观看网站| 国产精品久久久久久久电影 | av片东京热男人的天堂| 国内精品久久久久精免费| 欧美不卡视频在线免费观看| 伊人久久大香线蕉亚洲五| 国产一区二区激情短视频| 青草久久国产| 蜜桃久久精品国产亚洲av| 欧美精品啪啪一区二区三区| 午夜福利18| 成年女人永久免费观看视频| 国产精品嫩草影院av在线观看 | 人妻夜夜爽99麻豆av| 欧美在线黄色| 叶爱在线成人免费视频播放| 熟女人妻精品中文字幕| 欧美极品一区二区三区四区| 两性午夜刺激爽爽歪歪视频在线观看| 国内少妇人妻偷人精品xxx网站| 免费av毛片视频| 免费观看人在逋| 国产中年淑女户外野战色| 人妻丰满熟妇av一区二区三区| 午夜两性在线视频| 日韩欧美一区二区三区在线观看| 亚洲欧美日韩卡通动漫| 欧美中文日本在线观看视频| 国内精品久久久久精免费| 亚洲专区国产一区二区| 国产野战对白在线观看| 亚洲内射少妇av| 精品一区二区三区视频在线 | 久久亚洲精品不卡| 国产黄色小视频在线观看| 欧美黄色淫秽网站| 欧美高清成人免费视频www| 97碰自拍视频| 国产精品爽爽va在线观看网站| 成人性生交大片免费视频hd| 尤物成人国产欧美一区二区三区| 白带黄色成豆腐渣| www日本黄色视频网| 亚洲aⅴ乱码一区二区在线播放| 成熟少妇高潮喷水视频| 免费一级毛片在线播放高清视频| 久久香蕉精品热| 国产欧美日韩一区二区三| 久久亚洲真实| 欧美国产日韩亚洲一区| av天堂中文字幕网| 久久久成人免费电影| 色播亚洲综合网| 国产69精品久久久久777片| 久久久久免费精品人妻一区二区| 日本五十路高清| 18禁在线播放成人免费| 精品午夜福利视频在线观看一区| 婷婷亚洲欧美| 一进一出好大好爽视频| 男女视频在线观看网站免费| 日韩欧美在线二视频| 日韩国内少妇激情av| av天堂在线播放| 国产成人福利小说| 免费在线观看亚洲国产| 亚洲av中文字字幕乱码综合| 日韩欧美 国产精品| 少妇的逼好多水| 免费在线观看成人毛片| 香蕉久久夜色| 欧美av亚洲av综合av国产av| 成人午夜高清在线视频| 成熟少妇高潮喷水视频| 少妇熟女aⅴ在线视频| 男女做爰动态图高潮gif福利片| 亚洲av熟女| 欧美日韩亚洲国产一区二区在线观看| 国产毛片a区久久久久| 亚洲自拍偷在线| 免费看美女性在线毛片视频| x7x7x7水蜜桃| 国产伦一二天堂av在线观看| 91av网一区二区| 天天躁日日操中文字幕| 久久久久精品国产欧美久久久| 黄色成人免费大全| 久久久成人免费电影| 成人无遮挡网站| 色尼玛亚洲综合影院| 全区人妻精品视频| 精品国产美女av久久久久小说| 日本一二三区视频观看| 成人欧美大片| 日本与韩国留学比较| 91久久精品国产一区二区成人 | 亚洲av二区三区四区| 丰满的人妻完整版| 一本一本综合久久| 搡老熟女国产l中国老女人| 国产真实伦视频高清在线观看 | 国产极品精品免费视频能看的| 丰满人妻一区二区三区视频av | 国产高清videossex| 三级国产精品欧美在线观看| 亚洲内射少妇av| 午夜激情福利司机影院| 特级一级黄色大片| 少妇人妻一区二区三区视频| 欧美成人a在线观看| 国产aⅴ精品一区二区三区波| 两个人的视频大全免费| 色综合站精品国产| 亚洲欧美激情综合另类| 757午夜福利合集在线观看| 99久久精品热视频| 久9热在线精品视频| 亚洲在线自拍视频| 好男人电影高清在线观看| 亚洲性夜色夜夜综合| 国产97色在线日韩免费| 伊人久久精品亚洲午夜| 亚洲精品日韩av片在线观看 | 首页视频小说图片口味搜索| 成年版毛片免费区| 99久久无色码亚洲精品果冻| 99热这里只有是精品50| 日韩人妻高清精品专区| 观看美女的网站| 99热这里只有精品一区| 亚洲人成网站在线播放欧美日韩| 久久中文看片网| 精品熟女少妇八av免费久了| 少妇的逼水好多| 亚洲精品成人久久久久久| 在线观看午夜福利视频| 亚洲av成人精品一区久久| 成人国产一区最新在线观看| 免费看美女性在线毛片视频| 99久久综合精品五月天人人| 色老头精品视频在线观看| 久久99热这里只有精品18| 亚洲成av人片免费观看| netflix在线观看网站| 精华霜和精华液先用哪个| 亚洲黑人精品在线| 真实男女啪啪啪动态图| 日本黄色视频三级网站网址| 一个人免费在线观看的高清视频| 国产亚洲欧美98| 日本黄色视频三级网站网址| 天堂√8在线中文| 黄色片一级片一级黄色片| 成人一区二区视频在线观看| 97碰自拍视频| 国产高清有码在线观看视频| 中文字幕人妻熟人妻熟丝袜美 | 99riav亚洲国产免费| 人妻丰满熟妇av一区二区三区| 国产精品久久久久久亚洲av鲁大| 国产黄色小视频在线观看| 99热这里只有精品一区| 久久精品91蜜桃| 听说在线观看完整版免费高清| 亚洲av不卡在线观看| 久久国产精品影院| www.999成人在线观看| 国产成人欧美在线观看| 亚洲精品在线美女| 乱人视频在线观看| 国产精品亚洲av一区麻豆| 亚洲av二区三区四区| 国产一区二区亚洲精品在线观看| 亚洲av不卡在线观看| 亚洲一区二区三区色噜噜| 一个人免费在线观看电影| 级片在线观看| 人人妻,人人澡人人爽秒播| 十八禁人妻一区二区| 99久久99久久久精品蜜桃| 久久久久久久久大av| 国产亚洲欧美98| 网址你懂的国产日韩在线| 有码 亚洲区| 蜜桃亚洲精品一区二区三区| 欧美日韩黄片免| а√天堂www在线а√下载| 国产成人a区在线观看| 国产成人系列免费观看| 国产亚洲精品久久久com| 欧美日韩瑟瑟在线播放| 久久香蕉国产精品| 日本三级黄在线观看| 悠悠久久av| 成人特级av手机在线观看| 精品一区二区三区视频在线 | 久久中文看片网| 男人和女人高潮做爰伦理| 特级一级黄色大片| 免费电影在线观看免费观看| 蜜桃亚洲精品一区二区三区| 悠悠久久av| 老鸭窝网址在线观看| 亚洲午夜理论影院| 美女黄网站色视频| 国产单亲对白刺激| 亚洲一区二区三区色噜噜| av视频在线观看入口| 色视频www国产| 国产精品电影一区二区三区| 亚洲成人中文字幕在线播放| 真人一进一出gif抽搐免费| 欧美激情久久久久久爽电影| 亚洲精品粉嫩美女一区| 久久伊人香网站| 国产男靠女视频免费网站| 男插女下体视频免费在线播放| 国内精品一区二区在线观看| 一个人免费在线观看的高清视频| 中文在线观看免费www的网站| 亚洲国产精品999在线| 亚洲成人久久爱视频| 夜夜看夜夜爽夜夜摸| 欧美在线一区亚洲| 99久久99久久久精品蜜桃| 国产成年人精品一区二区| 精品午夜福利视频在线观看一区| 日本黄色视频三级网站网址| 欧美zozozo另类| 免费av观看视频| 久久久久久久久久黄片| 成年女人毛片免费观看观看9| 精品人妻偷拍中文字幕| 午夜福利在线观看吧| 一个人看视频在线观看www免费 | 一本精品99久久精品77| 午夜福利18| 日韩有码中文字幕| 亚洲国产高清在线一区二区三| 亚洲国产精品成人综合色| 成人永久免费在线观看视频| 免费高清视频大片| 老师上课跳d突然被开到最大视频 久久午夜综合久久蜜桃 | 国产欧美日韩精品一区二区| 国产精品香港三级国产av潘金莲| x7x7x7水蜜桃| 国产高清视频在线播放一区| 亚洲精品久久国产高清桃花| 在线十欧美十亚洲十日本专区| 精品久久久久久久毛片微露脸| 亚洲欧美日韩高清在线视频| 久久婷婷人人爽人人干人人爱| 性欧美人与动物交配| 色哟哟哟哟哟哟| 亚洲狠狠婷婷综合久久图片| 亚洲男人的天堂狠狠| 国产三级中文精品| 欧美性猛交╳xxx乱大交人| 日韩欧美国产在线观看| 中文字幕人成人乱码亚洲影| 99热6这里只有精品| 国产精品综合久久久久久久免费| 久久久久国内视频| 亚洲av二区三区四区| 亚洲成人免费电影在线观看| 亚洲av熟女| 色老头精品视频在线观看| 狠狠狠狠99中文字幕| 成人高潮视频无遮挡免费网站| 性色av乱码一区二区三区2| 老司机深夜福利视频在线观看| 超碰av人人做人人爽久久 | 宅男免费午夜| 欧美日本视频| 亚洲欧美一区二区三区黑人| 国内精品美女久久久久久| 国产精品亚洲av一区麻豆| 宅男免费午夜| 亚洲一区二区三区不卡视频| 亚洲五月婷婷丁香| 亚洲美女黄片视频| xxx96com| 午夜亚洲福利在线播放| 精品一区二区三区视频在线观看免费| 亚洲一区高清亚洲精品| 岛国视频午夜一区免费看| 久久久久免费精品人妻一区二区| 国产成人影院久久av| 国产av在哪里看| 精品人妻一区二区三区麻豆 | 免费在线观看日本一区| 91在线精品国自产拍蜜月 | 性欧美人与动物交配| 精品乱码久久久久久99久播| 国产伦在线观看视频一区| 黄片大片在线免费观看| 亚洲自拍偷在线| 黑人欧美特级aaaaaa片| 一夜夜www| 麻豆久久精品国产亚洲av| 欧美一级毛片孕妇| 亚洲熟妇熟女久久| 欧美不卡视频在线免费观看| 一进一出抽搐动态| 婷婷六月久久综合丁香| 99久久99久久久精品蜜桃| 一边摸一边抽搐一进一小说| 亚洲精品一区av在线观看| 久久人妻av系列| 国内毛片毛片毛片毛片毛片| 国产极品精品免费视频能看的| 国产单亲对白刺激| 黄片小视频在线播放| 中文资源天堂在线| 亚洲成av人片免费观看| 欧美在线黄色| 天堂av国产一区二区熟女人妻| 嫩草影院入口| 日本五十路高清| 久久国产精品人妻蜜桃| 51午夜福利影视在线观看| 1000部很黄的大片| 欧洲精品卡2卡3卡4卡5卡区| 九九热线精品视视频播放| 欧美区成人在线视频| 国产69精品久久久久777片| 日本精品一区二区三区蜜桃| 757午夜福利合集在线观看| 精品熟女少妇八av免费久了| 九九久久精品国产亚洲av麻豆| 99久久99久久久精品蜜桃| 男人的好看免费观看在线视频| 少妇人妻一区二区三区视频| 免费人成在线观看视频色| 老鸭窝网址在线观看| 人人妻,人人澡人人爽秒播| 久久人妻av系列| 91在线精品国自产拍蜜月 | 一卡2卡三卡四卡精品乱码亚洲| 亚洲国产欧美人成| 女同久久另类99精品国产91| 女人被狂操c到高潮| 哪里可以看免费的av片| 亚洲精品一卡2卡三卡4卡5卡| 国产日本99.免费观看| 美女高潮喷水抽搐中文字幕| 欧美精品啪啪一区二区三区| 18禁黄网站禁片午夜丰满| 久久久久久久久大av| 一个人看视频在线观看www免费 | 国产一区在线观看成人免费| 波多野结衣巨乳人妻| 俺也久久电影网| 五月伊人婷婷丁香| 亚洲中文日韩欧美视频| 日本一二三区视频观看| 午夜福利欧美成人| 中文字幕av在线有码专区| 男人的好看免费观看在线视频| 亚洲精品成人久久久久久| 一区二区三区国产精品乱码| 俺也久久电影网| 国产真人三级小视频在线观看| 中文字幕av在线有码专区| 禁无遮挡网站| 久久精品夜夜夜夜夜久久蜜豆| 国产亚洲精品av在线| 1024手机看黄色片| 日韩欧美一区二区三区在线观看| 国产伦一二天堂av在线观看| 成人av在线播放网站| 久久久久久久久大av| 成人高潮视频无遮挡免费网站| 少妇人妻精品综合一区二区 | 国内毛片毛片毛片毛片毛片| 成人av一区二区三区在线看| 99久久成人亚洲精品观看| 一夜夜www| 波多野结衣高清无吗| 高潮久久久久久久久久久不卡|