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

        ?

        一維非線性功能梯度材料的熱整流反轉(zhuǎn)*

        2015-12-19 05:28:38黃桂芳詹斯琦黃維清
        關(guān)鍵詞:熱流量熱傳導(dǎo)聲子

        黃桂芳,詹斯琦,黃維清

        (湖南大學(xué) 物理與微電子科學(xué)學(xué)院,湖南 長(zhǎng)沙 410082)

        近年來,一維非線性晶格體系的熱傳導(dǎo)性質(zhì)引起了廣泛的關(guān)注.人們?cè)噲D利用解析或非平衡分子動(dòng)力學(xué)方法從微觀動(dòng)力學(xué)機(jī)理的角度去解釋宏觀的熱傳導(dǎo)現(xiàn)象[1-6].盡管至今沒能得到一個(gè)完善的關(guān)于低維材料的熱傳導(dǎo)理論,但在經(jīng)典一維系統(tǒng)的熱傳導(dǎo)研究過程中,發(fā)現(xiàn)許多新奇的現(xiàn)象,如熱整流效應(yīng),界面熱阻,負(fù)微分熱阻等.這些現(xiàn)象的存在使得熱量的載體——聲子能像電子那樣得到控制與管理,因此開辟了熱二極管[7],熱晶體管[8],熱邏輯門[9]等熱學(xué)器件的研究領(lǐng)域.其中熱二極管的研究最為廣泛,人們?cè)噲D建立各種模型來實(shí)現(xiàn)材料的熱整流性質(zhì)[10-12].經(jīng)過大量的理論研究后,總結(jié)出熱二極管模型需要滿足兩個(gè)必要條件:系統(tǒng)的非線性(簡(jiǎn)諧)及不對(duì)稱結(jié)構(gòu).此外,實(shí)驗(yàn)方面也取得了不錯(cuò)的成績(jī),借助納米科技的迅速發(fā)展,Chang等人[13]首次在實(shí)驗(yàn)中成功制備出碳納米管固體熱二極管.隨后不同形態(tài)及材料構(gòu)成的熱二極管也相繼在實(shí)驗(yàn)中成功制備[14-15],這使熱二極管的實(shí)際應(yīng)用成為可能.

        功能梯度材料在各領(lǐng)域中得到廣泛的應(yīng)用,它是一種在結(jié)構(gòu)或性能上連續(xù)/準(zhǔn)連續(xù)變化的非均質(zhì)復(fù)合材料.對(duì)于功能梯度材料的研究,目前已逐步擴(kuò)大到機(jī)械、電子、聲學(xué)、光學(xué)、核工程、生物等領(lǐng)域.相比這些性質(zhì),人們對(duì)它的熱學(xué)性質(zhì)的研究并不夠多.文獻(xiàn)[16]最先將功能梯度材料概念引入一維非線性熱傳導(dǎo)領(lǐng)域,并提出一維線性質(zhì)量梯度下的Fermi-Pasta-Ulam(FPU)鏈,通過計(jì)算實(shí)現(xiàn)了熱整流效應(yīng).隨后線性質(zhì)量梯度下的Frenkel-Kontorova(FK)鏈及指數(shù)型質(zhì)量梯度下的FPU 鏈也相繼被提出用以構(gòu)建熱二極管模型[17-18].在實(shí)驗(yàn)方面,Chang等人[13]也是通過采用質(zhì)量梯度沉積制備出了不對(duì)稱的無定型C6H16Pt固體熱二極管.由此可見,功能梯度材料在低維非線性熱傳導(dǎo)領(lǐng)域極具應(yīng)用價(jià)值.

        本文利用非平衡分子動(dòng)力學(xué)模擬方法,系統(tǒng)地研究一維非線性功能梯度材料中不對(duì)稱熱傳導(dǎo)性質(zhì).在以往的研究中,改變系統(tǒng)的熱整流方向需要進(jìn)行結(jié)構(gòu)重建.然而,本文的模型能通過調(diào)節(jié)熱浴溫差進(jìn)行熱整流方向的控制,即系統(tǒng)出現(xiàn)了熱整流反轉(zhuǎn)現(xiàn)象.熱整流的反轉(zhuǎn)是由負(fù)微分熱阻引起的,并且它還與系統(tǒng)的體系大小及界面耦合強(qiáng)度緊密相關(guān).鑒于該模型能在原地控制熱整流方向,文章進(jìn)一步討論了其實(shí)際應(yīng)用的可能性.

        1 模型和計(jì)算方法

        圖1為系統(tǒng)結(jié)構(gòu)示意圖,該系統(tǒng)由兩段質(zhì)量連續(xù)變化的FK鏈組成,并由一彈性系數(shù)為kint的簡(jiǎn)諧彈簧連接,我們稱該彈簧為界面彈簧.系統(tǒng)通過引入線性遞減分布的質(zhì)量梯度來實(shí)現(xiàn)不對(duì)稱性(功能梯度),為方便觀察,質(zhì)量大小由晶格的大小表示.最左端與最右端晶格分別與溫度為TL和TR的熱浴接觸.

        系統(tǒng)的哈密頓量為:

        式中:N為系統(tǒng)總粒子數(shù);HLFK和HRFK分別為左右兩段FK 鏈的哈密頓量.

        圖1 系統(tǒng)結(jié)構(gòu)示意圖Fig.1 Schematic picture of the system

        式中:pi和xi分別為第i個(gè)晶格的動(dòng)量和位置.取k=a=1,分別表示晶格間的簡(jiǎn)諧作用強(qiáng)度及晶格常數(shù).第i個(gè)粒子的質(zhì)量為:

        固定Mmax=20,Mmin=1,Mmax為最左端晶格的質(zhì)量,Mmin為最右端晶格的質(zhì)量.基底勢(shì)強(qiáng)度為V=5.

        模擬計(jì)算中,采用固定邊界條件,并使兩端的晶格與NOSE-HOOVER 熱浴耦合,保持溫度分別為TL和TR.為簡(jiǎn)單起見,定義系統(tǒng)的平均溫度為T0=(TL+TR)/2,晶格鏈兩端的溫度差為ΔT=TL-TR=2T0Δ.當(dāng)Δ>0時(shí),最左端晶格與溫度較高的熱浴耦合,最右端的晶格與溫度較低的熱浴耦合;當(dāng)Δ<0時(shí),情況則相反.系統(tǒng)的局域化溫度定義為格點(diǎn)i處晶格的平均動(dòng)能:Ti=m〈˙xi2〉.系統(tǒng)的局域熱流定義為Ji=k〈˙xi(xi-xi-1)〉,(〈〉均表示對(duì)時(shí)間的平均).計(jì)算中,使用改進(jìn)后的龍格庫(kù)塔算法對(duì)晶格的運(yùn)動(dòng)方程進(jìn)行積分,設(shè)置步長(zhǎng)為0.1,并有足夠長(zhǎng)的計(jì)算時(shí)間(107)使系統(tǒng)最終達(dá)到非平衡穩(wěn)態(tài).當(dāng)系統(tǒng)的熱傳導(dǎo)過程達(dá)到非平衡穩(wěn)態(tài)以后,系統(tǒng)中流過不同格點(diǎn)的局域熱流量相同且為一恒定值.我們定義從左到右流過系統(tǒng)的熱流量的絕對(duì)值為J+,當(dāng)熱浴反轉(zhuǎn),從右到左的熱流量的絕對(duì)值則為J-,r=J+/J- 為系統(tǒng)的熱整流效率.

        2 結(jié)果與討論

        2.1 不對(duì)稱熱傳導(dǎo)

        本文模型中,當(dāng)界面耦合強(qiáng)度為kint=1.0時(shí),系統(tǒng)為一完整的FK 鏈.線性質(zhì)量梯度下的FK 鏈已被證明具有熱整流效應(yīng)[17].因此,我們先研究界面耦合強(qiáng)度較小時(shí)系統(tǒng)的熱傳導(dǎo)性質(zhì).圖2為當(dāng)N=100,kint=0.05時(shí),系統(tǒng)在不同平均溫度下的溫度分布圖.由圖2 可見,系統(tǒng)的溫度分布是不對(duì)稱的,并且由于界面耦合彈簧的存在,所有情況中,系統(tǒng)都存在由界面熱阻(ITR)引起的溫度躍變.而且熱浴反轉(zhuǎn)后,界面熱阻是不對(duì)稱的,正是因?yàn)榻缑鏌嶙璧牟粚?duì)稱性導(dǎo)致了系統(tǒng)的不對(duì)稱熱傳導(dǎo)現(xiàn)象.此外,該現(xiàn)象還與系統(tǒng)的平均溫度有關(guān),溫度過高或過低都不利于實(shí)現(xiàn)不對(duì)稱熱傳導(dǎo),因此我們固定系統(tǒng)的平均溫度為T0=0.1.

        圖2 不同平均溫度下系統(tǒng)的歸一化溫度分布圖Fig.2 Normalized temperature profile of the system for different average temperature

        2.2 熱整流反轉(zhuǎn)

        質(zhì)量梯度下的FK 鏈,沿晶格質(zhì)量增加方向的熱流量要大于沿質(zhì)量減小方向的熱流量[17].然而,當(dāng)FK 鏈中引入一彈性系數(shù)較小的界面彈簧后,這種單向的不對(duì)稱熱傳導(dǎo)現(xiàn)象將被打破.圖3 為當(dāng)kint=0.05時(shí),3種不同體系大小下熱流量與熱浴溫差因子的關(guān)系圖.由圖3(a)(N=100)可以看到,沿晶格質(zhì)量減小方向的熱流量J+隨|Δ|的增加幾乎呈線性增長(zhǎng),而沿質(zhì)量增加方向的熱流在時(shí)達(dá)到最大值后隨著|Δ|的增加而減小,導(dǎo)致了時(shí)J+<J-到J+>J-的熱整流反轉(zhuǎn)現(xiàn)象.因此系統(tǒng)僅在條件下,沿著晶格質(zhì)量增加方向的熱流量大于沿質(zhì)量減小方向的熱流量,而當(dāng)|Δ|>0.4時(shí),則轉(zhuǎn)為沿質(zhì)量減小方向的熱流量大于沿質(zhì)量增加方向的熱流量.可見,熱整流方向的改變僅僅是通過控制熱浴兩端溫差實(shí)現(xiàn)的.一般而言,熱流量會(huì)隨著熱浴溫差的增加而增大,然而在圖3(a)中,當(dāng)|Δ|>0.5時(shí)出現(xiàn)了負(fù)導(dǎo)熱現(xiàn)象,正因?yàn)樗某霈F(xiàn)直接導(dǎo)致了熱整流的反轉(zhuǎn).并且負(fù)導(dǎo)熱現(xiàn)象出現(xiàn)的區(qū)域隨著體系的增大逐步移向更高的熱浴溫差區(qū).因此,當(dāng)體系大小為N=1 000時(shí),從J+<J-到J+>J-的熱整流的反轉(zhuǎn)現(xiàn)象出現(xiàn)在|Δ|≈0.7,如圖3(b)所示.隨著體系進(jìn)一步增大,負(fù)導(dǎo)熱現(xiàn)象最終消失,熱整流反轉(zhuǎn)也隨之消失,因此,當(dāng)N=2 000時(shí),系統(tǒng)中沿著晶格質(zhì)量增加方向的熱流量大于沿質(zhì)量減小方向的熱流量,即J+<J- .為了進(jìn)一步研究界面耦合強(qiáng)度及體系大小對(duì)熱傳導(dǎo)的影響,圖4給出了系統(tǒng)在以上3種體系下熱流量J±與界面耦合強(qiáng)度kint的關(guān)系,固定熱浴溫差因子.所有的情況下,隨著界面耦合強(qiáng)度的加強(qiáng),J±迅速增加,并且在一定范圍內(nèi)呈現(xiàn)J±~k2int的關(guān)系,這種關(guān)系普遍存在于聲子通過一較弱耦合界面的情況中.隨著kint進(jìn)一步增加,J±將逐漸達(dá)到飽和,二者的關(guān)系也趨近穩(wěn)定.特別是,在圖4(a)和(b)中還存在J+>J-到J+<J-的交叉現(xiàn)象,即當(dāng)界面耦合強(qiáng)度較小時(shí),系統(tǒng)沿晶格質(zhì)量減小方向的熱流量大于沿質(zhì)量增加方向的熱流量,隨著耦合強(qiáng)度的增加,變?yōu)檠刭|(zhì)量增加方向的熱流量大于沿質(zhì)量減小方向的熱流量.而當(dāng)體系較大時(shí)(圖4(c)),系統(tǒng)的熱整流方向不受界面耦合強(qiáng)度的影響,沿晶格質(zhì)量增加方向的熱流量始終大于沿質(zhì)量減小方向的熱流量,即J+<J-.可見,無論體系多大,當(dāng)界面彈簧彈性系數(shù)較大時(shí),均有J+<J-,這是因?yàn)楫?dāng)系統(tǒng)處于強(qiáng)耦合狀態(tài)時(shí),兩段FK 鏈組成一個(gè)整體,成為一條完整的FK 鏈,與文獻(xiàn)[17]中得出的熱整流方向一致,且系統(tǒng)的熱整流方向不隨體系的改變而改變.而當(dāng)系統(tǒng)處于弱耦合狀態(tài)時(shí),系統(tǒng)的熱整流方向會(huì)隨著熱浴溫差,體系大小的改變而改變,出現(xiàn)熱整流反轉(zhuǎn)現(xiàn)象.

        2.3 熱整流機(jī)理

        針對(duì)以上不對(duì)稱熱傳導(dǎo)現(xiàn)象,我們根據(jù)有效聲子理論研究了其熱整流機(jī)理:當(dāng)FK 鏈中的晶格處于高溫極限時(shí),其動(dòng)能足以克服基底勢(shì)能的束縛,因此基底勢(shì)可忽略不計(jì)并且系統(tǒng)可近似地看做為簡(jiǎn)諧晶格鏈,晶格中聲子的振動(dòng)頻率范圍為0<ω/2π<;當(dāng)晶格處于低溫極限時(shí)則被困在基底勢(shì)壘中,通過線性化運(yùn)動(dòng)方程可得到晶格中聲子的振動(dòng)頻率范圍為

        因此,當(dāng)系統(tǒng)處于強(qiáng)耦合狀態(tài)時(shí),系統(tǒng)回歸為一條完整的FK 鏈,不存在界面熱阻,我們主要研究熱浴兩端晶格的聲子頻率:當(dāng)系統(tǒng)左端質(zhì)量較大的晶格與高溫?zé)嵩〗佑|時(shí),晶格處于高溫極限中,聲子的振動(dòng)頻率為0<ω/2π<0.071.而與低溫?zé)嵩〗佑|的右端質(zhì)量較小的晶格中,聲子的頻率為0.35<ω/2π<0.47.可見左右兩端晶格中聲子的頻率不重合,不能匹配,這就意味著聲子不能較順利地從左往右輸運(yùn)熱量,因此系統(tǒng)中的熱流量較小.

        圖3 不同體系下熱流量J 與|Δ|的關(guān)系圖Fig.3 Dependence of the heat flux Jon|Δ|for different system size

        圖4 不同體系下熱流量J 與kint的關(guān)系圖Fig.4 Dependence of the heat flux J on kintfor different system size

        當(dāng)系統(tǒng)左端晶格與低溫?zé)嵩〗佑|時(shí),其聲子頻率為0.079<ω/2π<0.106.而右端晶格處于高溫極限,其聲子頻率為0<ω/2π<0.308.可見左右兩端晶格中聲子的頻率能較好地匹配,更有利于聲子的熱輸運(yùn),因此,J+<J-,這與我們的結(jié)果吻合.然而,當(dāng)系統(tǒng)處于弱耦合狀態(tài)時(shí),以上結(jié)論無法解釋熱整流的反轉(zhuǎn).但可以很直觀地看出,負(fù)導(dǎo)熱現(xiàn)象直接導(dǎo)致了熱整流反轉(zhuǎn),而負(fù)導(dǎo)熱現(xiàn)象被認(rèn)為是由負(fù)微分熱阻所引起,因此需要考慮界面效應(yīng).正是由于彈性系數(shù)較小的界面彈簧存在,當(dāng)系統(tǒng)右端質(zhì)量較小粒子與低溫?zé)嵩〗佑|時(shí),導(dǎo)致系統(tǒng)出現(xiàn)了負(fù)微分熱阻,即隨著系統(tǒng)右端的溫度進(jìn)一步減小,界面兩端晶格中聲子頻率呈現(xiàn)由匹配到分離的現(xiàn)象.于是J-隨著熱浴溫差的增加而減小,熱整流方向由最初的J+<J-變?yōu)镴+>J-,因此熱整流反轉(zhuǎn)了.由此可見,只有界面耦合強(qiáng)度較小情況下,才能通過調(diào)節(jié)熱浴溫差控制熱整流的方向.

        此外,熱整流反轉(zhuǎn)除了與耦合強(qiáng)度有關(guān)外,還與體系的大小密切相關(guān).由圖3和圖4可知,只有在體系較小的條件下才能出現(xiàn)熱整流反轉(zhuǎn),當(dāng)體系大于一定的范圍后,無論界面耦合強(qiáng)度如何,系統(tǒng)都呈現(xiàn)出和強(qiáng)耦合狀態(tài)下一樣的熱整流性質(zhì),沒有熱整流反轉(zhuǎn)現(xiàn)象.因此,圖3(c)中熱整流反轉(zhuǎn)的消失可理解如下:界面彈簧對(duì)系統(tǒng)的影響隨著體系的增大而減小,在較大的體系中,界面彈簧對(duì)系統(tǒng)的作用可以忽略不計(jì),系統(tǒng)的兩部分仍可被看做強(qiáng)耦合在一起的整體.因此,可以總結(jié)得出,熱整流反轉(zhuǎn)的根本原因是界面彈簧的引入,并且只有在特定的體系下,當(dāng)系統(tǒng)為弱耦合狀態(tài)時(shí),才能通過調(diào)節(jié)熱浴溫差控制熱整流方向.

        2.4 熱整流效率

        值得一提的是,對(duì)于特定的體系(如N=100),系統(tǒng)的熱整流效率r=J+/J- 可從r→0.94變?yōu)閞→2.1,與目前實(shí)驗(yàn)中制備的固體熱二極管的熱整流效率較為一致.這意味著,該模型在熱整流方向控制領(lǐng)域具有潛在的應(yīng)用價(jià)值.此外,由于負(fù)微分熱阻的存在以及功能梯度材料的廣泛應(yīng)用,該模型也是構(gòu)建熱晶體管,熱邏輯門等熱學(xué)器件不錯(cuò)的候選者.

        3 結(jié) 論

        本文研究了質(zhì)量梯度下一維非線性功能梯度材料的熱整流性質(zhì).該系統(tǒng)在特定的界面耦合強(qiáng)度及特定的體系大小范圍內(nèi),可以僅通過控制熱浴溫差實(shí)現(xiàn)由J+<J-到J+>J-的熱整流反轉(zhuǎn).引起熱整流反轉(zhuǎn)的根本原因是系統(tǒng)中存在的界面彈簧.鑒于該系統(tǒng)無需進(jìn)行結(jié)構(gòu)重建便能在原地進(jìn)行熱整流方向的控制,討論了它在實(shí)驗(yàn)制備及應(yīng)用方面的可能性.

        [1]LEPRI S,LIVI R,POLITI A.Heat conduction in chains of nonlinear oscillators[J].Physical Review Letters,1997,78(10):1896-1899.

        [2]HU Bam-bi,LI Bao-wen,ZHAO Hong.Heat conduction in one-dimensional chains[J].Physical Review E,1998,57(3):2992-2995.

        [3]LEPRI S,LIVI R,POLITI A.Thermal conduction in classical low dimensional lattices[J].Physics Reports-Review Section of Physics Letters,2003,377(1):1-80.

        [4]PROSEN T,CAMPBELL D K.Normal and anomalous heat transport in one-dimensional classical lattices [J].Chaos,2005,15(1):015117.

        [5]GAUL C,BUETTNER H.Quantum mechanical heat transport in disordered harmonic chains[J].Physical Review E,2007,76(1):011111.

        [6]LI Nian-bei,REN Jie,WANG Lei,etal.Colloquium:Phononics:Manipulating heat flow with electronic analogs and beyond[J].Reviews of Modern Physics,2012,84(3):1045-1066.

        [7]TERRANEO M,PEYRARD M,CASATI G.Controlling the energy flow in nonlinear lattices:A model for a thermal rectifier[J].Physical Review Letters,2002,88(9):094302.

        [8]LI Bao-wen,WANG Lei,CASATI G.Negative differential thermal resistance and thermal transistor[J].Applied Physics Letters,2006,88(14):143501.

        [9]WANG Lei,LI Bao-wen.Thermal logic gates:Computation with phonons[J].Physical Review Letters,2007,99(17):177208.

        [10]LI Bao-wen,WANG Lei,CASATI G.Thermal diode:Rectification of heat flux[J].Physical Review Letters,2004,93(18):184301.

        [11]YANG Nuo,ZHANG Gang,LI Bao-wen.Thermal rectification in asymmetric graphene ribbons[J].Applied Physics Letters,2009,95(3):033107.

        [12]YANG Nuo,ZHANG Gang,LI Bao-wen.Carbon nanocone:A promising thermal rectifier[J].Applied Physics Letters,2008,93(24):243111.

        [13]CHANG C W,OKAWA D,MAJUMDAR A,etal.Solidstate thermal rectifier[J].Science,2006,314(5802):1121-1124.

        [14]KOBAYASHI W,TERAOKA Y,TERASAKI I,An oxide thermal rectifier[J].Applied Physics Letters,2009,95(17):171905.

        [15]TIAN He,XIE Dan,YANG Yi,etal.A novel solid-state thermal rectifier based on reduced graphene oxide[J].Scientific Reports,2012,2:523.

        [16]YANG Nuo,LI Nian-bei,WANG Lei,etal.Thermal rectification and negative differential thermal resistance in lattices with mass gradient[J].Physical Review B,2007,76(2):020301.

        [17]HU Tao,HU Ke,TANG Yi.Anharmonic on-site potential implies asymmetry thermal transport in 1D mass-graded harmonic lattice[J].Physica B-Condensed Matter,2010,405(21):4407-4412.

        [18]SHAH TEJAL N,GAJJAR P N.Study of thermal conductivity and thermal rectification in exponential mass graded lattices[J].Physics Letters A,2012,376(8):438-441.

        猜你喜歡
        熱流量熱傳導(dǎo)聲子
        基于泵閥聯(lián)合調(diào)節(jié)的供暖系統(tǒng)節(jié)能優(yōu)化運(yùn)行
        加熱型織物系統(tǒng)的熱傳遞性能
        半無限板類聲子晶體帶隙仿真的PWE/NS-FEM方法
        一類三維逆時(shí)熱傳導(dǎo)問題的數(shù)值求解
        納米表面聲子 首次實(shí)現(xiàn)三維成像
        聲子晶體覆蓋層吸聲機(jī)理研究
        熱傳導(dǎo)方程解的部分Schauder估計(jì)
        一類非線性反向熱傳導(dǎo)問題的Fourier正則化方法
        基于聲子晶體理論的導(dǎo)線防舞方法及數(shù)值驗(yàn)證
        低溫貯箱連接支撐結(jié)構(gòu)優(yōu)化設(shè)計(jì)
        載人航天(2016年2期)2016-05-24 07:49:22
        成人影院免费观看在线播放视频| 麻豆亚洲一区| 久久人妻无码一区二区| 97精品超碰一区二区三区| 天天天天躁天天爱天天碰| 成人做爰高潮尖叫声免费观看| 极品 在线 视频 大陆 国产| 久久亚洲精精品中文字幕早川悠里 | 亚洲国产精品sss在线观看av| 亚洲欧美另类自拍| 国产精品福利小视频| 最新国产一区二区三区| 亚洲天堂亚洲天堂亚洲色图| 夜夜春亚洲嫩草影院| 激情亚洲一区国产精品| 日韩亚洲制服丝袜中文字幕| 国产av在线观看91| 国产一区二区三区日韩精品| 国产 麻豆 日韩 欧美 久久 | 狠狠噜天天噜日日噜| 中文字幕亚洲好看有码| 亚洲国产一区二区,毛片| 在线观看 国产一区二区三区| 欧洲女人与公拘交酡视频| 亚洲成av人片天堂网| 国内精品大秀视频日韩精品| 91在线观看国产自拍| 午夜少妇高潮在线观看| 妺妺窝人体色www聚色窝| 处破痛哭a√18成年片免费| 精品无码一区二区三区小说| 亚洲丰满熟女乱一区二区三区| 久久久精品国产免大香伊| 久久精品人人爽人人爽| 国产视频嗯啊啊啊| 丝袜美腿人妻第一版主| 国自产精品手机在线观看视频| XXXXBBBB欧美| 在线女同免费观看网站| 91九色老熟女免费资源| 午夜三级a三级三点|