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

        ?

        三氨基三硝基苯基高聚物粘結炸藥熱力學性質(zhì)的理論計算研究*

        2019-06-04 05:31:44范航何冠松楊志劍聶福德陳鵬萬
        物理學報 2019年10期
        關鍵詞:界面振動

        范航 何冠松 楊志劍 聶福德 陳鵬萬

        1)(北京理工大學機電學院,北京 100081)

        2)(中國工程物理研究院化工材料研究所,綿陽 621900)

        1 引 言

        高聚物粘結炸藥(PBX)是一類在武器行業(yè)有著重要應用的顆粒復合含能材料,具有機械感度低、加工性能良好、力學性能優(yōu)異等特點,在現(xiàn)代高性能戰(zhàn)斗部系統(tǒng)中得到廣泛應用[1,2].在PBX炸藥中,以三氨基三硝基苯(TATB)為基的PBX在武器中有著廣泛的應用前景.TATB的機械感度和熱感度相比于其他含能材料有著明顯的優(yōu)勢,而且其熔點高,耐熱性能好,具有非線性光學性質(zhì)[3-5],是唯一通過美國能源部鈍感高能炸藥11項安全鑒定的高能炸藥.美國已率先將TATB基PBX用于核戰(zhàn)斗部,從而顯著提高了核戰(zhàn)斗部的安全性.在武器服役的過程中,外界復雜的力-熱環(huán)境,不僅要求PBX部件具有良好的結構完整性,還要具備較強的環(huán)境適應性,因此對PBX的力學性能提出了設計要求[6].另一方面,在復雜的熱-力環(huán)境作用下,熱傳導及由于熱作用引起的內(nèi)應力都可能導致PBX產(chǎn)生內(nèi)部損傷,嚴重的可能導致結構破壞.但是由于PBX炸藥是一類高聚物混合材料,主要由高能單質(zhì)炸藥、高聚物黏結劑和相關添加劑組成,在界面包覆處不可避免地存在缺陷,這些缺陷的存在不僅使界面成為制約PBX熱傳導的重要因素,而且這些PBX內(nèi)部缺陷在受外部刺激時易產(chǎn)生熱點,還會帶來嚴重的安全性問題[7,8].因此對高聚物粘結炸藥的熱力學性質(zhì),尤其是炸藥界面熱傳導的理論研究,對于炸藥感度的理解和PBX的設計都有著重要意義.

        國內(nèi)外僅有少數(shù)研究團隊開展了含能材料熱力學性質(zhì)的理論計算研究,且多集中于單質(zhì)含能材料.Gee等[9]和Bedrov等[10]分別開發(fā)了針對TATB分子動力學(MD)計算的力場參數(shù);Stevens等[11]和Rykounov[12]研究了TATB的狀態(tài)方程;Kroonblawd和Sewell[13,14]利用分子動力學方法系統(tǒng)研究了TATB的各向異性熱導;Fan等[15]采用德拜聲子理論方法研究了TATB熱導與壓力和溫度的關系;蔣文燦等[16]、Wu等[17]和Liu等[18]也基于第一性原理方法研究了TATB的聲子振動模式.目前,單質(zhì)TATB的狀態(tài)方程、熱容、熱膨脹、熱導率等性質(zhì)目前已較為完善,但是關于PBX等復合含能材料的熱導率、界面熱阻的研究則鮮有報道.Long和Chen[19]利用分子動力學方法最早研究了奧克托今(HMX)基PBX的界面熱導率.該方法基于線性響應理論和德拜聲子理論,認為聲子在界面處的散射率等同于彈性波的散射率,以此得到了HMX和氟聚物(F2312,F2313),TATB等界面處波的散射率.界面熱導在低頻區(qū)正比于振動頻率的平方,在高頻區(qū)因分子振動模的耦合而變得無規(guī)律.Long和Chen[20]采用第一性原理方法,通過計算TATB和石墨的聲子譜,得到TATB和石墨界面的聲子折射率和界面熱導率,并界定出五個主導TATB/石墨界面熱導率的入射聲子模.總體來看,目前針對含能材料的熱力學性質(zhì)的理論計算研究已有初步進展,開展了部分含能材料的聲子譜研究,這對于今后從晶格動力學角度開展熱導率計算有著很好的助益.

        含能材料的熱力學理論計算主要有第一性原理和分子動力學計算兩種途徑.但是,受制于含能材料結構的復雜性,目前用于炸藥分子的力場參數(shù)的可靠性尚未得到有效檢驗.第一性原理計算是目前能準確針對復雜體系含能材料進行聲子研究的有效方法,但需要注意的是該方法仍在計算能力方面存在著一定的制約.本文對于TATB和PVDF的聲子色散問題將采用第一性原理計算方法進行.而在熱傳導問題上,分子動力學方法通過牛頓力學方程描述體系中原子的相互作用,相較于第一性原理的計算更為簡單和直觀,且計算的結果與實驗符合較好.本文用分子動力學方法進行在TATB熱導率的計算.第一性原理和分子動力學方法可以從微觀尺度上解釋其結構與性能之間的關系.這不僅有助于含能材料的配方設計,還可縮短研制周期并提高研制過程中的安全性.

        本文的研究基于晶格動力學理論[21],結合第一性原理計算和分子動力學方法系統(tǒng)探討了TATB基PBX炸藥的熱力學性質(zhì).首先采用凍聲子法計算得到TATB和PVDF的聲子色散曲線和熱容等熱力學性質(zhì);根據(jù)散射錯配模型,詳細討論了TATB與高聚物界面的聲子透射率和界面熱導隨溫度的變化情況;最后,結合平衡分子動力學獲得的熱導率數(shù)據(jù),對TATB基PBX炸藥的熱導率做了簡要分析.研究結果為PBX混合炸藥的熱力學參數(shù)研究方法的構建,以及對理解PBX中熱傳導的本質(zhì)可以提供理論上的解釋.

        2 理論方法與計算模型

        2.1 第一性原理計算模型

        第一性原理計算的TATB和PVDF結構見圖1.其中TATB(C6H6N6O6)晶體屬三斜晶系,單胞中包含兩個TATB分子,共48個原子.TATB分子的構型與石墨層狀結構類似.PVDF是一種高分子聚合物材料,目前實驗已發(fā)現(xiàn)α,β,γ,δ四種晶相的九種構型形式[22,23],其中β相單胞中有兩條全反式分子鏈,分子鏈中相鄰的CH2和CF2具有相同朝向.目前針對β相的研究最為廣泛,為便于結果比較,在此我們選取β相構型為代表構型,進行本文的聲子譜計算和界面熱導率研究.針對上述兩種晶體,構建了2×2×2的TATB超胞和3×3×3的PVDF超胞分別用于聲子譜的計算.

        圖1 TATB和PVDF晶體結構圖Fig.1.The crystal structure of TATB (a),and PVDF (b).

        TATB和PVDF的聲子譜由凍聲子法[24]計算得到,結構優(yōu)化和力的計算利用VASP[25]軟件進行.計算中,我們選取投影綴加平面波(PAW)贗勢[26,27]結合Perdew-Burke-Ernzerhof (PBE)交換關聯(lián)泛函[28,29]的方法.其平面波截斷能為550 eV,能量的收斂標準為10—6eV,結構優(yōu)化過程中力的收斂標準為10—4eV/?,K點的設置為4×4×4.二階力常數(shù)和聲子譜的處理均通過phonopy[30]軟件進行.

        2.2 簡諧晶格動力學

        聲子是晶格振動的簡正模式的量子化描述,在晶格的動力學行為和熱力學性質(zhì)中起著重要的作用,晶體中的聲子色散曲線可以依據(jù)簡諧晶格動力學理論[31,32]計算.晶格動力學理論認為,原子能夠在格點附近發(fā)生振動.對于包含N個原子的晶體體系,其勢能函數(shù)可以在平衡位置附近展開成泰勒級數(shù)的形式:

        式中α,β,γ是笛卡爾坐標的三個方向;Φ0,Φα(lκ),Φαβ(lκ,l′κ′),Φαβγ(lκ,l′κ′,l′′κ′′),··· 分別對應0階、一階、二階、三階等力常數(shù).對于固定體積的晶體體系,我們對體系內(nèi)的原子施加一個微小的位移,原子的振動問題可由二階項的簡諧近似得以解決.通過構建體系的動力學矩陣D(q):

        在動力學矩陣做對角化處理后,即可得到簡諧近似中原子的振動頻率.

        2.3 界面熱傳導模型

        雖然基于傅里葉定律的傳熱理論宏觀上可以使用,但界面的研究多涉及微觀的原子或分子尺度,此時界面聲子的散射影響加大,必須從微觀角度重新考慮界面?zhèn)鳠釂栴}.首先定義界面處的凈熱流密度和界面熱阻[33],其計算公式分別為

        在單位時間、單位面積內(nèi)物體1通過界面?zhèn)鬏數(shù)轿矬w2的能量為

        式中,j為聲子模,N1為物體1的聲子密度,ω為聲子頻率,θ為聲子入射角,? 為約化普朗克常數(shù),ωd為聲子截斷頻率,c1,j為物體1中聲子模為j的聲速,a1→2(θ,j,ω)為聲子透射率.通過(6)式可以發(fā)現(xiàn),界面聲子透射率a1→2(θ,j,ω)是界面熱輸運的關鍵參數(shù).根據(jù)Swartz和Pohl提出的散射失配模型(DMM)[33],我們認為所有聲子在入射到界面處時發(fā)生隨機散射,并且與入射聲子的初始狀態(tài)無關.聲子通過界面的概率與該聲子對應的態(tài)密度成正比,且與入射角θ無關.聲子透射率α1→2(ω)的計算公式有

        結合(6)式得到界面熱導率

        目前,DMM模型已在界面熱導率的計算上廣泛應用,計算結果也較接近于實驗值.不過需要指出的是DMM模型預測的數(shù)值是界面熱阻的下限,其真實數(shù)值會根據(jù)界面的粗糙程度有不同程度的偏高.

        2.4 平衡態(tài)分子動力學

        平衡態(tài)分子動力學(EMD)計算熱導率的工作,主要是基于格林-久保公式進行[34].格林-久保公式建立了非平衡過程的輸運系數(shù)與平衡態(tài)中相應物理量的漲落之間的聯(lián)系.輸運系數(shù)等于自相關函數(shù)對于時間的積分.因此熱導率的計算公式[35,36]為

        式中,kB表示玻爾茲曼常數(shù),T和V表示體系溫度和體積,Jμ表示在μ方向上的熱流,t′ 表示關聯(lián)時間.

        本文的計算使用由Gee等[9]提出的TATB全原子力場進行分子動力學計算.上述力場已在研究中廣為應用[37,38],其模擬結果與實驗匹配較好.力場詳細參數(shù)見文獻[9].我們構建了3×4×5的超胞用于分子動力學模擬.平衡態(tài)分子動力學的計算主要利用LAMMPS(large-scale atomic molecular massively parallel simulator)模擬程序包[39].在分子動力學的計算中,計算時間步長設定為0.5 fs,力的收斂標準為10—6kcal/mol·?.首先利用分子力學對建立的模型進行構型優(yōu)化,確保體系能量最小化.隨后在300 K溫度及大氣壓條件下進行等壓等溫(NPT)分子動力學模擬,使體系得到充分的弛豫,模擬時長為100 ps.最后使體系在微正則(NVE)系綜下繼續(xù)模擬,總時長為1000 ps,關聯(lián)時長為100 ps.

        3 計算結果與討論

        3.1 聲子譜計算

        聲子色散關系是固體的一個重要物理量,可以反映出晶格振動模式與動量、能量的關系.聲子色散關系中,Γ點處的晶格振動與材料的紅外和Raman特性有關.基于凍聲子法,我們得到了TATB,PVDF的聲子色散關系和聲子振動模式結果(見圖2(a)和圖2(b)以及圖3(a)和圖(b)).TATB晶體的單胞有48個原子,共有144個振動模式,其中包含3個聲學支和141個光學支.

        對比表1和圖2(c)中理論計算結果和實驗數(shù)據(jù),本文采用廣義梯度近似(GGA)交換關聯(lián)泛函的計算結果與其他研究組采用不同交換關聯(lián)泛函的計算結果符合較好.交換關聯(lián)泛函及考慮范德瓦耳斯修正等因素對計算結果的影響不大.在高頻區(qū),尤其是大于3000 cm—1時,GGA的結果會略小于Jiang(局域梯度近似(LDA))的計算結果,更接近于實驗值.TATB晶格比熱的計算數(shù)值與其他計算結果相符.其中在溫度293 K時計算的比熱(498.3 J/mol·K)與實驗值(495.7 J/mol·K)[41]的誤差為0.52%.結合圖3和表1分析TATB的晶格振動情況,可以發(fā)現(xiàn)在低頻區(qū)( < 300 cm—1)尤其是頻率低于100 cm—1的區(qū)域,TATB晶體中硝基(—NO2)的振動是主要的振動模式.隨著振動頻率提高,在300—1600 cm—1區(qū)間,苯環(huán)的振動貢獻有著明顯的提升,同時硝基(—NO2)和氨基(—NH2)的振動也有著一定程度的貢獻.在高頻區(qū)(> 3000 cm—1),TATB晶格振動的貢獻絕大部分來自于氨基(-NH2)的振動.伴隨著晶格的振動,固體內(nèi)的傳熱問題隨之而生.在此,我們基于振動模式的結果,簡要分析材料熱導率這一表征傳熱能力的重要的物理參數(shù).熱導率是聲子弛豫時間、聲子群速度和聲子比熱的函數(shù).聲子群速度和聲子比熱可由簡諧晶格動力學獲得,而聲子弛豫時間的獲得,則需要通過非簡諧晶格動力學方法求得三階等高階力常數(shù).TATB等復雜材料的高階力常數(shù)計算仍較為困難,目前仍暫未獲得.我們在此僅對TATB和PVDF的聲子比熱進行分析,以初步探究TATB和PVDF的導熱情況.晶格動力學中晶格熱容與聲子態(tài)密度g(ω)有關,公式如下:

        圖2 (a),(b)TATB和PVDF聲子色散曲線;(c)TATB晶格熱容隨溫度的變化;(d),(e)TATB和PVDF第一布里淵區(qū)高對稱點積分路徑Fig.2.The phonon dispersion relation of TATB (a),and PVDF (b);(c)the heat capacity of TATB as a function of temperature;the Brillouin zone and high symmetry points of TATB (d),and PVDF (e).

        圖3 (a),(b)TATB和PVDF的Γ點聲子模分析;(c),(d)TATB和PVDF的加權聲子態(tài)密度和聲子態(tài)密度(小圖)Fig.3.Decomposition of the gamma-point eigenmodes of TATB (a),and PVDF (b);the weighted phonon density of states of TATB (c),and PVDF (d).

        定義聲子的加權態(tài)密度為

        根據(jù)圖3(c)和圖3(d),我們分析了TATB和PVDF在300 K溫度時的聲子加權態(tài)密度情況.可以看到,盡管TATB和PVDF在分別在100 THz和90 THz附近有著較大的態(tài)密度,但是經(jīng)過權重分析后,二者的高頻聲子比例明顯減少,低頻聲子在導熱中占據(jù)著主導貢獻.

        另外需要注意到的是,我們的計算結果中TATB在部分高對稱點附近存在著一定的虛頻,該現(xiàn)象在蔣文燦等[16]和Wu等[17]的結果同樣也有出現(xiàn).虛頻的產(chǎn)生,有計算中未考慮非簡諧效應和溫度效應的原因,凍聲子計算超胞的結構較小等.本文計算中,TATB的2×2×2超胞有384個原子.受制于結構的復雜性,目前計算資源無法高效針對更大超胞做出計算.與已有實驗和理論數(shù)據(jù)對比,本文的聲子譜計算結果具有可靠性,可用于后續(xù)的界面聲子散射研究.

        3.2 界面熱導

        界面熱導率的計算結果如圖4所示.從圖4(a)可以看出,界面熱導隨溫度的升高而相應提高,并在高溫時趨于平穩(wěn).在低溫區(qū)間,界面兩側的晶格中只有少部分聲子模式被激發(fā),隨著溫度的升高,更多的聲子模式隨之被激發(fā),進而參與界面處的傳熱,當溫度接近或達到德拜溫度時,被激發(fā)的聲子模式數(shù)量接近飽和,聲子模對界面熱導的貢獻不再增加.進一步分析,根據(jù)聲子透射率公式(7)式,可以發(fā)現(xiàn)界面的熱傳導主要發(fā)生在界面兩側振動頻率重疊的地方,且與重疊程度大小有關.圖4(b)展示了聲子透射率和累積界面熱導率隨頻率的變化關系,TATB-PVDF界面的熱傳導集中在聲子頻率小于45 THz的范圍內(nèi).其中在0—5 THz區(qū)間,累積界面熱導率已經(jīng)達到總熱導率的50%,隨后的聲子貢獻程度逐漸下降,在大于45 THz的高頻區(qū)間,聲子對界面熱傳導幾乎無貢獻.界面熱導在聲子態(tài)密度重疊程度最大處達到最大值.另外,聲子從TATB側到PVDF側的透射率,明顯大于反向的情況.從圖4(b)中可以看出,聲子透射率在2,10和20 THz附近出現(xiàn)了拐點,這與圖3(c)TATB的加權聲子態(tài)密度的變化趨勢基本一致.DMM模型認為,界面的聲子透射率與入射聲子的聲子模、聲子群速度和入射角無關,聲子通過界面的概率與該聲子對應的態(tài)密度成正比,因此我們并沒有分析各聲子模對聲子透射率的影響.TATB聲子主導著TATB-PVDF界面的聲子透射,也進一步主導著界面的熱傳導.

        表1 TATB部分Raman活性模計算數(shù)值與文獻結果的比較Table 1.Comparison of the Raman modes of TATB crystal obtained in the present and previous calculations with experimental results.

        圖4 (a)TATB-PVDF界面熱導率隨溫度變化情況;(b)TATB-PVDF界面聲子透射率(黑)和累積熱導率(紅)Fig.4.(a)The interface thermal conductance of TATBPVDF as a function of temperature;(b)phonon transmission of TATB-PVDF interface as a function of frequency.

        以上的研究表明,PBX炸藥各組分的聲子譜計算,可以作為指導炸藥配方設計的一種手段.通過計算各組分聲子譜,我們可以根據(jù)振動頻率匹配程度,初步篩選出界面熱阻小的配方.另外,以上研究也預示著通過高分子功能化、接枝手段,改變界面兩側振動頻率的情況,可以獲得更高界面熱導率的配方.目前在實驗方面,我們尚無有效測量手段獲知炸藥界面的熱導數(shù)據(jù),而且PBX炸藥界面的理論研究也不多見,我們較完整地研究了PBX炸藥的界面熱傳導問題.

        3.3 熱導率

        我們通過EMD方法獲得TATB的熱導率為0.58 W/mK,結果見圖5(a),計算細節(jié)參考2.4節(jié).EMD的交換關聯(lián)時長為100 ps,本文TATB的計算共重復10組,并對最后的熱導率結果做平均.可以看到交換關聯(lián)時間在25 ps逐漸收斂.由于TATB的低熱導率,其聲子平均自由程也僅有0.54 nm,遠小于硅、鍺等中高熱導率材料的自由程.因此,我們用于EMD模擬的參數(shù),包括模型尺寸和交換關聯(lián)時間均較合理.其300 K熱導率計算結果也與實驗值0.54 W/mK[42]符合較好.此外我們選取PVDF的熱導率實驗值0.13 W/mK,和3.2節(jié)中獲得的界面熱導率進行PBX熱導率計算.

        圖5 (a)EMD方法計算TATB熱導率結果;(b)TATB基PBX炸藥熱導率隨顆粒尺寸的變化Fig.5.(a)EMD calculation of TATB thermal conductivity;(b)thermal conductivity of TATB-based PBXs as a function of particle size.

        對于復合材料的熱導率,通常我們考慮采用串并聯(lián)的基本模型預測熱導.結合3.2節(jié)中DMM模型所獲得的界面熱導率,我們基于串聯(lián)模型計算得到了TATB-PVDF高聚物粘結炸藥的熱導率.該體系的熱導率由三部分貢獻組成,分別是主炸藥TATB、黏結劑PVDF和TATB-PVDF界面的熱導率.當我們從某一方向?qū)Τ叽鐬長的炸藥體系施加一個熱流q,則該體系的溫度變化 ΔT可由式(13)得到:

        式中,ΔTc,ΔTa,ΔTI分別為主炸藥、黏結劑和界面的溫度差;kc,ka,kI為對應三者的熱導率;d表示主炸藥的粒徑尺寸;δd為包覆厚度.

        PBX炸藥的熱導由(14)式得出,

        將(13)式代入(14)式,并認為包覆的厚度δd極小可忽略不計,最終化簡得到PBX的熱導率

        圖5展示了TATB-PVDF高聚物粘結炸藥的熱導率隨顆粒尺寸的變化情況.當顆粒尺寸大于100 nm時,界面熱阻對于PBX熱導率的影響有限.此外,由于TATB的平均聲子自由程僅為5 ?左右,界面層的厚度遠大于聲子自由程,因此TATB-PVDF界面的熱傳導將不會受到聲子彈道輸運的限制,界面熱傳導滿足宏觀的傅里葉擴散理論.考慮到PBX的實際結構,其粘結劑包覆的炸藥顆粒是一個諸多炸藥晶體的集合體,顆粒間的接觸熱阻、以及顆粒尺寸也都會影響PBX的熱導.目前暫無有效實驗數(shù)據(jù)進行比較,本研究采用的的串聯(lián)傳熱模型,雖預測的是復合材料熱導率的下限,但是考慮到PBX炸藥中主成分的質(zhì)量分數(shù)超過90%且該部分熱導率明顯高于粘結劑組分的事實,該模型仍然適用.在分析顆粒尺寸對熱導率影響的問題上,串聯(lián)模型則可以給出不受界面熱阻影響的顆粒臨界最小尺寸,對PBX設計有著指導意義.復合含能材料有效熱導率模型研究,以及利用分子動力學方法模擬真實結構PBX熱導的工作將在接下來的工作中開展.

        4 結 論

        本文基于晶格動力學,結合第一性原理和分子動力學方法,系統(tǒng)研究了TATB基高聚物黏結炸藥的聲子色散關系、界面透射系數(shù)和熱導率.

        1)計算了TATB和PVDF的聲子色散關系,計算結果與實驗符合較好.TATB在低頻區(qū)主要以硝基(-NO2)的振動為主,中頻區(qū)苯環(huán)的振動貢獻有著明顯的提升,而在高頻區(qū)TATB晶格振動的貢獻絕大部分來自于氨基(-NH2)的振動.低頻聲子在導熱中有著明顯的主導貢獻.

        2)利用DMM模型和聲子色散關系,研究了TATB與PVDF界面的熱傳導.熱導率隨溫度升高而上升,并且在高溫情況下接近不變.通過提高界面兩側的振動頻率的重合程度,可以改善界面熱導.DMM模型僅考慮了界面二聲子的彈性散射情況,多聲子的非彈性散射修正的界面熱傳導將在接下來的研究進行.

        3)基于串聯(lián)模型,分析了界面熱導率和顆粒尺寸對PBX熱導率的關系.當顆粒尺寸大于100 nm時,界面熱阻對于PBX熱導率的影響有限.現(xiàn)實情況中,由于PBX炸藥中主炸藥顆粒并非以單晶狀態(tài)存在,其粘結劑包覆的炸藥顆粒是一個諸多炸藥小顆粒的集合體,顆粒間的接觸熱阻、以及顆粒尺寸都會影響PBX的熱導.

        鑒于炸藥分子的復雜結構,PBX的大規(guī)模第一性原理計算目前仍難以實現(xiàn).目前計算結果的準確程度仍有提高的空間.本文構建了PBX炸藥界面熱導的研究方法,對炸藥配方研究、PBX設計有著指導作用.

        猜你喜歡
        界面振動
        振動的思考
        科學大眾(2023年17期)2023-10-26 07:39:14
        噴水推進高速艇尾部振動響應分析
        國企黨委前置研究的“四個界面”
        當代陜西(2020年13期)2020-08-24 08:22:02
        This “Singing Highway”plays music
        振動攪拌 震動創(chuàng)新
        中國公路(2017年18期)2018-01-23 03:00:38
        中立型Emden-Fowler微分方程的振動性
        基于FANUC PICTURE的虛擬軸坐標顯示界面開發(fā)方法研究
        空間界面
        金秋(2017年4期)2017-06-07 08:22:16
        電子顯微打開材料界面世界之門
        人機交互界面發(fā)展趨勢研究
        久久亚洲av成人无码电影| 日韩中文字幕网站| 免费在线av一区二区| 日本中文一区二区在线| 大地资源中文第3页| 丁香五月缴情综合网| 日本韩国三级aⅴ在线观看| 国产丝袜美腿一区二区三区| 人妻中文字幕在线网站| 思思久久96热在精品国产| 欧美久久久久中文字幕| 丰满少妇又爽又紧又丰满动态视频| 极品尤物人妻堕落沉沦| 亚洲精品美女久久久久久久| 在线视频这里只有精品| 日本高清一区二区在线观看| 国内精品亚洲成av人片| 亚洲国产精华液网站w| 真正免费一级毛片在线播放| 日本一区二区三区的免费视频观看 | 欧美中文字幕在线看| 中文乱码字幕人妻熟女人妻| 图片小说视频一区二区| 欧美两根一起进3p做受视频| 手机看片国产日韩| 久久开心婷婷综合中文| 二区三区亚洲精品国产| 好看的日韩精品视频在线| 久久久久国产综合av天堂| 成人在线激情网| 免费啪啪av人妻一区二区 | 浪货趴办公桌~h揉秘书电影| 久久av无码精品人妻出轨| 国产视频不卡在线| 亚洲国产综合精品中久| 北条麻妃国产九九九精品视频 | 最新国产不卡在线视频| 成人爽a毛片在线视频| 亚洲AV无码专区国产H小说| 99久久精品人妻一区| 欧美老妇多毛xxxxx极瑞视频|