亚洲免费av电影一区二区三区,日韩爱爱视频,51精品视频一区二区三区,91视频爱爱,日韩欧美在线播放视频,中文字幕少妇AV,亚洲电影中文字幕,久久久久亚洲av成人网址,久久综合视频网站,国产在线不卡免费播放

        ?

        Klein-Gordon-Schr?dinger方程的幾種差分格式及比較

        2024-01-25 03:30:12林周瑾汪佳玲霍昱安
        關(guān)鍵詞:能量守恒內(nèi)積將式

        林周瑾 汪佳玲 霍昱安

        摘要: 探究在特定的初值和邊界條件下一維Klein-Gordon-Schr?dinger方程的幾種差分格式并進(jìn)行比較。利用經(jīng)典的向前差分算子、中心差分算子、Crank-Nicolson方法和緊差分算子分別為Klein-Gordon-Schr?dinger方程構(gòu)造向前Euler式、Crank-Nicolson格式及緊差分格式。結(jié)果表明:Crank-Nicolson格式及緊差分格式能夠精確地保持離散電荷和能量守恒。數(shù)值實(shí)驗(yàn)驗(yàn)證了理論結(jié)果的正確性。

        關(guān)鍵詞:Klein-Gordon-Schr?dinger方程; 向前Euler格式; Crank-Nicolson格式; 緊差分格式; 電荷守恒; 能量守恒

        中圖分類號(hào): O 241.82文獻(xiàn)標(biāo)志碼: A?? 文章編號(hào): 1000-5013(2024)01-0108-13

        Several Difference Schemes and Comparisons for ?Klein-Gordon-Schr?dinger Equation

        LIN Zhoujin, WANG Jialing, HUO Yu′an

        (School of Mathematics and Statistics, Nanjing University of Information Science and Technology, Nanjing 210044, China)

        Abstract: Several difference schemes of one-dimensional Klein-Gordon-Schr?dinger equation under specific initial value and boundary conditions are investigated and contrasted. The classical forward difference operator, central difference operator, Crank-Nicolson method and compact difference operator are used to construct forward Euler scheme, Crank-Nicolson scheme and compact difference scheme respectively. Results show that Crank-Nicolson scheme and the compact difference scheme can accurately conserve the discrete charge and energy conservation. The correctness of the theoretical result has been verified by numerical experiments.

        Keywords: Klein-Gordon-Schr?dinger equation; forward Euler scheme; Crank-Nicolson scheme; compact difference scheme; charge conservation; energy conservation

        Klein-Gordon-Schr?dinger(KGS)方程是薛定諤方程的狹義相對(duì)論形式,該系統(tǒng)于1970年被Yukawa首次提出。1975年,由Fukuda和Tsutsumi提出了帶有Yukawa作用的KGS系統(tǒng)模型[1],被用來(lái)描述量子場(chǎng)理論中守恒復(fù)標(biāo)量核子場(chǎng)與實(shí)標(biāo)量介子場(chǎng)之間相互作用,是相對(duì)論量子力學(xué)和量子場(chǎng)論中的最基本方程。隨著學(xué)術(shù)科研的發(fā)展與科學(xué)技術(shù)的創(chuàng)新,KGS方程的研究越來(lái)越受到國(guó)內(nèi)外學(xué)者的重視。在過(guò)去的二十年中,許多學(xué)者們針對(duì)KGS方程的解析解和數(shù)值解進(jìn)行了一系列的研究。

        在數(shù)學(xué)方面,F(xiàn)ukuda等[1]討論了三維空間中耦合的KGS方程的初邊值問(wèn)題,建立了初邊值問(wèn)題整體解的存在唯一性定理。Baillon等[2]討論了耦合的KGS方程的柯西問(wèn)題,并且證明了KGS方程柯西問(wèn)題的唯一整體解的存在性。Darwish等[3]設(shè)計(jì)了一種代數(shù)方法來(lái)統(tǒng)一構(gòu)造KGS耦合方程的一系列顯式精確解。Wang等[4]用雅可比橢圓函數(shù)展開法的推廣得到KGS方程的周期波解。文獻(xiàn)[5-10]也在數(shù)學(xué)上對(duì)KGS方程展開研究。

        然而,該方程的解析解很難得到,大多數(shù)情形只能靠數(shù)值方法進(jìn)行求解。因此,對(duì)于如何得到能夠長(zhǎng)時(shí)間地保持系統(tǒng)解的行為的KGS方程的數(shù)值解就顯得尤為重要。在數(shù)值方面,學(xué)者們利用許多不同的數(shù)值方法對(duì)KGS方程進(jìn)行了數(shù)值計(jì)算[11-17]。Wang[11]提出一個(gè)緊差分格式來(lái)計(jì)算具有齊次Dirichlet邊界條件的KGS方程。通過(guò)連接合適的辛Runge-Kutta-type方法和辛Runge-Kutta-Nystr?m-type方法,Hong等[12]提出了KGS方程的顯式多辛格式,并證明用該方法構(gòu)造的方法是多辛的,可在適當(dāng)?shù)倪吔鐥l件下精確地保持離散電荷守恒定律。Wang等[13]提出用傅里葉譜方法求解具有周期邊界條件的空間分?jǐn)?shù)階KGS方程,并且表明該格式可以保持離散電荷和能量守恒。

        基于此,本文在一定的初值和邊值條件下,利用不同的差分格式求解一維KGS方程并進(jìn)行比較。

        1 數(shù)值格式的構(gòu)造

        在區(qū)域Ω=[a,b]×[0,T]上考慮一維KGS方程,即

        選取初值條件

        φ(x,0)=φ0(x),? u(x,0)=u0(x),? ut(x,0)=u1(x),? x∈[a,b](3)

        和Dirichlet零邊界條件

        φ(a,t)=φ(b,t)=0,? u(a,t)=u(b,t)=0,? t∈(0,T]。(4)

        式(3)中:φ0(x)是給定的具有足夠光滑性的復(fù)值函數(shù);u0(x)和u1(x)是兩個(gè)給定的具有足夠光滑性的實(shí)值函數(shù),這3個(gè)函數(shù)充當(dāng)求解過(guò)程中的初始解。

        式(1)~(4)具有電荷守恒律和能量守恒律,即

        定義空間

        V0h={v|v={vj|0≤j≤J}∈Vh,v0=vJ=0}

        和三對(duì)角矩陣

        其中:矩陣A根據(jù)二階中心差分算子可得,矩陣B為對(duì)角占優(yōu)矩陣,因此是可逆矩陣。

        設(shè)u,v∈Vh,定義離散內(nèi)積和離散范數(shù),即

        2 幾種差分格式

        2.1 向前Euler格式

        令Φnj=φ(xj,tn),Unj=u(xj,tn)。在節(jié)點(diǎn)(xj,tn)處考慮KGS方程(1)~(2),有

        式(7),(8)中:0≤n≤N-1;1≤j≤J-1。

        由向前差分算子及二階中心差分算子,有

        將式(9)~(12)代入式(7),(8),得到

        結(jié)合式(3),(4),可得

        忽略式(13),(14)的小量項(xiàng),則有

        并用φnj,unj分別代替Φnj,Unj,得到差分格式為

        式(15)~(18)即為KGS方程的向前Euler格式。稱R(1)j,n和R(2)j,n為差分格式(15)和差分格式(16)的局部截?cái)嗾`差。記

        則可知截?cái)嗾`差R(1)j,n,R(2)j,n滿足

        |R(1)j,n|≤c1(τ+h2),? 0≤n≤N-1, 1≤j≤J-1,|R(2)j,n|≤c2(τ2+h2),? 1≤n≤N-1, 1≤j≤J-1。

        其中:c1,c2是與h和τ無(wú)關(guān)的常數(shù)。

        注1 向前Euler格式(15)~(18)是一個(gè)非線性顯性格式,并且該格式下φ的數(shù)值解在時(shí)間方向和空間方向上分別具有1階和2階精度,u的數(shù)值解在時(shí)間方向和空間方向上都具有2階精度。

        2.2 Crank-Nicolson格式

        其中:0≤n≤N-1;1≤j≤J-1。

        應(yīng)用公式

        可得到

        再利用式(9)及

        可以得到

        在(xj,tn+1/2)處考慮方程(2),即

        式(24)中:1≤n≤N-1;1≤j≤J-1。

        結(jié)合式(11),(12),(21)及

        可將式(24)改寫為

        略去式(23)和式(28)的小量項(xiàng),則有

        結(jié)合初值條件(3)和邊值條件(4),并用φnj,unj分別代替Φnj,Unj,得到Crank-Nicolson差分格式為

        式(29)~(32)即為KGS方程的Crank-Nicolson格式。稱R(3)j,n和R(4)j,n為差分格式(23)和差分格式(28)的局部截?cái)嗾`差。記

        則可知截?cái)嗾`差R(3)j,n,R(4)j,n滿足

        |R(3)j,n|≤c3(τ2+h2),? 0≤n≤N-1, 1≤j≤J-1,|R(4)j,n|≤c4(τ2+h2),? 1≤n≤N-1, 1≤j≤J-1。

        其中:c3,c4是與h和τ無(wú)關(guān)的常數(shù)。

        注2 與向前Euler格式不同,Crank-Nicolson格式(29)~(32)是一個(gè)非線性隱性格式,并且該格式下φ和u的數(shù)值解在時(shí)間方向和空間方向上都具有2階精度。

        2.3 緊差分格式

        在點(diǎn)(xj,tn+1/2)處考慮方程(1),有

        其中:0≤n≤N-1;1≤j≤J-1。

        結(jié)合式(19)~(22),有

        式(33)中:0≤n≤N-1;1≤j≤J-1。

        式(33)兩邊同時(shí)左乘緊差分算子Ah,可以得到

        由于有

        所以有

        將式(35)代入式(34),有

        在點(diǎn)(xj,tn+1/2)處考慮方程(2),有

        式(37)中:1≤n≤N-1;1≤j≤J-1。

        將式(25)~(27)代入式(37),可以得到

        將式(38)兩邊同時(shí)左乘緊差分算子Ah,并利用式(11)及

        可得

        略去式(36)和式(39)中的小量項(xiàng),則有

        結(jié)合初值條件(3)和邊值條件(4),并且用φnj,unj分別代替Φnj,Unj,得到KGS方程的緊差分格式為

        式(40)~(43)即為KGS方程的緊差分格式。稱R(5)j,n和R(6)j,n為差分格式(36)和差分格式(39)的局部截?cái)嗾`差。記

        則可知截?cái)嗾`差R(5)j,n,R(6)j,n分別滿足

        |R(5)j,n|≤c5(τ2+h4),? 0≤n≤N-1, 1≤j≤J-1,

        |R(6)j,n|≤c6(τ2+h4),? 1≤n≤N-1, 1≤j≤J-1。

        其中:c5,c6是與h和τ無(wú)關(guān)的常數(shù)。

        注3 緊差分格式(40)~(43)也是一個(gè)非線性隱性格式,并且該格式下的φ和u的數(shù)值解在時(shí)間方向和空間方向上分別具有2階和4階精度。

        3 守恒性

        引理1[18] 對(duì)于任意的u,v∈V0h,有〈δ2xu,v〉=-〈δ+xu,δ+xv〉。

        引理2[18] 對(duì)于任意的u∈Vh,n=0,1,…,N-1,則有

        Re〈-B-1A(un+1+un),(un+1-un)〉=|||δxun+1|||2-|||δxun|||2,

        Im〈B-1A(un+1+un),(un+1-un)〉=0。

        其中:Re和Im分別表示取函數(shù)的實(shí)部和虛部。

        定理1 Crank-Nicolson格式(29)~(32)能夠精確地保持離散的電荷和能量守恒,即

        Qn=φn2≡Q0,? n=1,2,…,N,(44)

        En=un2+δ+tun2+δ+xun2+δ+xφn2-2〈un,|φn|2〉=E0,? n=1,2,…,N-1。(45)

        證明:式(29),(30)可以表示為

        將式(46)與φn+1/2作內(nèi)積,并取虛部,有

        由引理1可知

        〈δ2xφn+1/2,φn+1/2〉=-〈δ+xφn+1/2,δ+xφn+1/2〉=-δ+xφn+1/22∈R。

        因此,有

        故有

        φn+12=φn2。(48)

        因此,式(44)成立。

        將式(46)與φn+1-φn做內(nèi)積,并取實(shí)部,有

        對(duì)上式進(jìn)行逐項(xiàng)分析,即

        整理可以得到

        -δ+xφn+12+δ+xφn2+〈un+1+un,|φn+1|2〉-〈un+1+un,|φn|2〉=0。(49)

        將式(47)與un+1-un做內(nèi)積,有

        〈δ2tun+1/2,un+1-un-1〉-〈δ2xun+1/2,un+1-un〉+〈un+1/2,un+1-un〉-12〈|φn|2+|φn+1|2,un+1-un〉=0。

        對(duì)上式進(jìn)行逐項(xiàng)分析,可得到

        整理得到

        δ+tun+12-δ+tun2+δ+xun+12-δ+xun2+un+12-un2-〈|φn|2+|φn+1|2,un+1-un〉=0。(50)

        用式(50)減去式(49)得到

        un+12+δ+tun+12+δ+xun+12+δ+xφn+12-2〈un+1,|φn+1|2〉=un2+δ+tun2+δ+xun2+δ+xφn2-2〈un,|φn|2〉,

        因此,式(45)成立。

        定理2 緊差分格式(40)~(43)能夠精確保持離散電荷和能量守恒,即

        Qn=φn2≡Q0,? n=0,1,…,N,(51)

        En=un2+δ+tun2+|||δxun|||2+|||δxφn|||2-2〈un,|φn|2〉≡E0,? n=0,1,…,N-1。(52)

        證明:利用前面定義的矩陣A和B,式(40),(41)可以表示為

        將式(53)與φn+1/2作內(nèi)積,并取虛部,則有

        對(duì)上式進(jìn)行逐項(xiàng)分析,有

        又〈un+1/2φn+1/2,φn+1/2〉∈R,所以有

        由此可知,式(51)成立。

        將式(53)與δ+tφn作內(nèi)積,并取實(shí)部,有

        逐項(xiàng)分析,有

        則式(56)可表示為

        |||δxφn+1|||2-|||δxφn|||2-〈un+1+un,|φn+1|2〉+〈un+1+un,|φn|2〉=0。(57)

        將式(54)與δ+tun作內(nèi)積,并取實(shí)部,有

        分析式(58)的每一項(xiàng),可得

        則式(58)可以寫成

        δ+tun+12-δ+tun2+|||δxun+1|||2-|||δxun|||2+un+12-un2-〈|φn|2+|φn+1|2,un+1-un〉=0。(59)

        結(jié)合式(57)與式(59),有

        un+12+δ+tun+12+|||δxun+1|||2+|||δxφn+1|||2-2〈un+1,|φn+1|2〉=un2+δ+tun2+|||δxun|||2+|||δxφn|||2-2〈un,|φn|2〉。

        由此可知,式(52)成立。

        4 數(shù)值實(shí)驗(yàn)

        通過(guò)數(shù)值實(shí)驗(yàn)驗(yàn)證前面的理論結(jié)果。根據(jù)文獻(xiàn)[19]可以得到KGS方程的解析解,即

        式(60),(61)中:v為孤立波的傳播速度;x0為初始相位。對(duì)于固定的t,當(dāng)x→∞時(shí),φ(x,t)和u(x,t)迅速衰減到0。因此,在數(shù)值上可以在有限區(qū)域(a,b)中求解KGS方程。其中,-a,b1,邊界條件為零邊界。

        4.1 數(shù)值解

        考慮初值條件

        φ0(x)=φ(x,0,v,0),? u0(x)=u(x,0,v,0),? u1(x)=ut(x,0,v,0)。

        計(jì)算主要在區(qū)間[-20,20]中進(jìn)行,選取空間步長(zhǎng)h為0.2,時(shí)間步長(zhǎng)τ為0.001 s,傳播速度v為0.1。向前Euler格式、Crank-Nicolson格式和緊差分格式在數(shù)值運(yùn)算時(shí)間(T)分別為1,16 s時(shí)得到的數(shù)值解,如圖1~3所示。

        由圖1~3可知:當(dāng)T=16時(shí),向前Euler格式的數(shù)值解出現(xiàn)了一些輕微的振蕩,Crank-Nicolson格式和緊差分格式的數(shù)值解較為光滑。這表明相較于其他兩種穩(wěn)定的隱式格式,作為顯式格式的向前Euler格式相對(duì)不穩(wěn)定。

        當(dāng)T分別為1,16 s時(shí),分別運(yùn)用向前Euler格式、Crank-Nicolson格式和緊差分格式求解KGS方程時(shí)的CPU運(yùn)行時(shí)間(tCPU),結(jié)果如表1所示。由表1可知:顯式的向前Euler格式的計(jì)算速度明顯優(yōu)于隱式的Crank-Nicolson格式和緊差分格式,這是因?yàn)橄蚯癊uler格式在計(jì)算過(guò)程中沒有迭代。

        4.2 電荷守恒與能量守恒

        分別定義離散電荷誤差error Q和能量誤差error E為

        式(62)中:Qn和En分別表示第n步的電荷值和能量值。計(jì)算在區(qū)間[-20,20]中進(jìn)行,選取空間步長(zhǎng)h為0.2,時(shí)間步長(zhǎng)τ為0.001 s,傳播速度v為0.1。分別繪制T=10 s時(shí)Crank-Nicolson格式和緊差分格式的電荷、能量值及其離散誤差,如圖4,5所示。

        由圖4可知:電荷誤差和能量誤差分別在10-13和10-12左右,表明Crank-Nicolson格式的電荷和能量是守恒的。由圖5可知:緊差分格式的守恒量誤差分別在10-11和10-12左右,表明緊差分格式能夠很好地保持離散電荷和能量守恒。

        5 結(jié)束語(yǔ)

        利用經(jīng)典的差分算子為一維KGS方程分別構(gòu)造向前Euler格式、Crank-Nicolson格式和緊差分格式。利用相關(guān)理論知識(shí)討論了3種格式的精度,詳細(xì)證明了Crank-Nicolson格式和緊差分格式能夠精確保持離散電荷守恒及能量守恒。數(shù)值實(shí)驗(yàn)結(jié)果表明,與Crank-Nicolson格式和緊差分格式相比,向前Euler格式長(zhǎng)時(shí)間計(jì)算的穩(wěn)定性稍差,但是其計(jì)算效率更高。另外,在數(shù)值上,Crank-Nicolson格式和緊差分格式能夠精確地保持離散的電荷和能量守恒,驗(yàn)證了理論結(jié)果的正確性。

        通過(guò)對(duì)3種格式的比較,可以看出它們?cè)谇蠼釱GS方程時(shí)的優(yōu)缺點(diǎn),為不同工程應(yīng)用提供合適的選擇。

        參考文獻(xiàn):

        [1] FUKUDA I,TSUTSUMI M.On coupled Klein-Gordon-Schr?dinger equations Ⅱ[J].Journal of Mathematical Analysis and Applications,1978,66(2):358-378.DOI:10.1016/0022-247X(78)90239-1.

        [2] BAILLON J B,CHADAM J M.The cauchy problem for the coupled Schr?dinger-Klein-Gordon equations[J].North-Holland Mathematics Studies,1978,30:37-44.DOI:10.1016/S0304-0208(08)70857-0.

        [3] DARWISH A,F(xiàn)AN Engui.A series of new explicit exact solutions for the coupled Klein-Gordon-Schr?dinger equations[J].Chaos,Solitons and Fractals,2004,20(3):609-617.DOI:10.1016/S0960-0779(03)00419-3.

        [4] WANG Mingliang,ZHOU Yubin.The periodic wave solutions for the Klein-Gordon-Schr?dinger equations[J].Physics Letters A,2003,318(1/2):84-92.DOI:10.1016/j.physleta.2003.07.026.

        [5] WANG Baoxiang.Classical global solutions for nonlinear Klein-Gordon-Schr?dinger equations[J].Mathematical Methods in the Applied Sciences,1997,20(7):599-616.DOI:10.1002/(SICI)1099-1476(19970510)20:7<599::AID-MMA866>3.0.CO;2-7.

        [6] CAVALCANTI M,CAVALCANTI V.Global existence and uniform decay for the coupled Klein-Gordon-Schr?din ger equations[J].Nonlinear Differential Equations and Applications,2000,7(3):285-307.DOI:10.1007/PL0000 1426.

        [7] HAYASHI N,WAHL W.On the global strong solutions of coupled Klein-Gordon-Schr?dinger equations[J].Journal of the Mathematical Society of Japan,1987,39(3):489-497.DOI:10.2969/jmsj/03930489.

        [8] OHTA M.Stability of stationary states for the coupled Klein-Gordon-Schr?dinger equations[J].Nonlinear Analysis: Theory,Methods and Applications,1996,27(4):455-461.DOI:10.1016/0362-546X(95)00017-P.

        [9] FUKUDA I,TSUTSUMI M.On the Yukawa-coupled Klein-Gordon-Schr?dinger equations in three space dimensions[J].Proceedings of the Japan Academy Series A Mathematical Sciences,1975,51(6):402-405.DOI:10.3792/pja/1195518563.

        [10] XIA Jingna,HAN Shuxia,WANG Mingliang.The exact solitary wave solutions for the Klein-Gordon-Schr?dinger equations[J].Applied Mathematics and Mechanics,2002,23(1):58-64.DOI:10.1007/BF02437730.

        [11] WANG Tingchun.Optimal point-wise error estimate of a compact difference scheme for the Klein-Gordon-Schr? dinger equation[J].Journal of Mathematical Analysis and Applications,2014,412(1):155-167.DOI:10.1016/j.jmaa.2013.10.038.

        [12] HONG Jialin,JIANG Shanshan,LI Chun.Explicit multi-symplectic methods for Klein-Gordon-Schr?dinger equations[J].Journal of Computational Physics,2009,228(9):3517-3532.DOI:10.1016/j.jcp.2009.02.006.

        [13] WANG Junjie,XIAO Aiguo.Conservative Fourier spectral method and numerical investigation of space fractional Klein-Gordon-Schr?dinger equations[J].Applied Mathematics and Computation,2019,350:348-365.DOI:10.1016/j.amc.2018.12.046.

        [14] WANG Tingchun,ZHAO Xiaofei,JIANG Jiaping.Unconditional and optimal H2-error estimates of two linear and conservative finite difference schemes for the Klein-Gordon-Schr?dinger equation in high dimensions[J].Advances in Computational Mathematics,2018,44:477-503.DOI:10.1007/s10444-017-9557-5.

        [15] ZHANG Luming.Convergence of a conservative difference scheme for a class of Klein-Gordon-Schr?dinger equations in one space dimension[J].Applied Mathematics and Computation,2005,163(1):343-355.DOI:10.1016/j.amc.2004.02.010.

        [16] KONG Linghua,WANG Lan,JIANG Shanshan,et al.Multi-symplectic Fourier pseudo-spectral integrators for Klei n-Gordon-Schr?dinger equations[J].Science China Mathematics,2013,56:915-932.DOI:10.1007/s11425-013-4575-3.

        [17] XIANG Xinmin.Spectral method for solving the system of equations of Schr?dinger-Klein-Gordon field[J].Journal of Computational and Applied Mathematics,1988,21(2):161-171.DOI:10.1016/0377-0427(88)90265-8.

        [18] WANG Tingchun,ZHAO Xiaofei.Unconditional L∞-convergence of two compact conservative finite difference schemes for the nonlinear Schr?dinger equation in multi-dimensions[J].Calcolo,2018,55(3):34-59.DOI:10.1007/s10092-018-0277-0.

        [19] WANG Jialing,LIANG Dong,WANG Yushun.Analysis of a conservative high-order compact finite difference scheme for the Klein-Gordon-Schr?dinger equation[J].Journal of Computational and Applied Mathematics,2019,358:84-96.DOI:10.1016/j.cam.2019.02.018.

        猜你喜歡
        能量守恒內(nèi)積將式
        壓縮sCO2儲(chǔ)能系統(tǒng)動(dòng)態(tài)特性仿真
        AKNS方程的三線性型及周期孤立波解
        因子von Neumann代數(shù)上非線性*-Lie導(dǎo)子的刻畫
        單自由度系統(tǒng)
        運(yùn)用能量守恒,巧解物理難題
        動(dòng)量能量守恒齊用難題不難求解完勝
        電磁場(chǎng)能量守恒研究
        考試周刊(2017年3期)2017-02-13 00:17:15
        基于矩陣的內(nèi)積函數(shù)加密
        關(guān)于矩陣的Frobenius內(nèi)積的一個(gè)推廣
        阻尼系統(tǒng)的特征
        91在线观看国产自拍| 区二区欧美性插b在线视频网站| 日本欧美在线播放| 插入中文字幕在线一区二区三区| 亚洲最新精品一区二区| 绝顶潮喷绝叫在线观看| 国产乱子伦精品免费无码专区| 手机AV片在线| 免费一区二区在线观看视频在线| 24小时日本在线视频资源| 欧洲-级毛片内射| 国产丰满乱子伦无码专| 五月激情四射开心久久久| 97人伦色伦成人免费视频| 亚洲综合无码一区二区| 亚洲av永久无码精品成人| 国产在线精品成人一区二区三区| 五月丁香六月综合缴清无码| 特级毛片a级毛片免费播放| 无码91 亚洲| 久久一区二区国产精品| 久久99精品久久水蜜桃| 99爱这里只有精品| 区三区久久精品水蜜桃av| 国产日本精品一二三四区| 日本无码人妻波多野结衣| 亚洲VA中文字幕欧美VA丝袜| 亚洲av色香蕉一区二区三区软件| 成人无码av免费网站| 国产欧美日韩一区二区三区在线| 久久99精品久久久久九色| 精品高清一区二区三区人妖| 亚洲国产成人久久三区| 亚洲美腿丝袜综合一区| 蜜桃av一区在线观看| 久久久国产精品| 中国xxx农村性视频| 亚洲无码图| 91精品久久久老熟女91精品| 国产成人a人亚洲精品无码| 亚洲中文字幕在线一区二区三区|