999精品在线视频,手机成人午夜在线视频,久久不卡国产精品无码,中日无码在线观看,成人av手机在线观看,日韩精品亚洲一区中文字幕,亚洲av无码人妻,四虎国产在线观看 ?

無閥納米泵中水流的反常堵塞*

2024-05-13 07:41:30李偉健周曉艷陸杭軍
物理學(xué)報 2024年9期
關(guān)鍵詞:振動區(qū)域

李偉健 周曉艷 陸杭軍

(浙江師范大學(xué),凝聚態(tài)物理研究所,金華 321004)

對于顆粒物質(zhì),在錐形通道的窄口處容易發(fā)生堵塞現(xiàn)象,實驗中通常利用機械振動進行疏通.而對于無孔不入的水卻不同,即使在納米尺度的碳納米管中仍然可以快速通透.本文利用分子動力學(xué)模擬研究了由錐形碳納米管與石墨烯平面構(gòu)成的無閥納米泵,發(fā)現(xiàn)水輸運在一定條件下也會出現(xiàn)反常堵塞.與顆粒物質(zhì)截然不同的是振動方法無法恢復(fù)水流,相反,它促使水在納米水泵的通道窄口處發(fā)生堵塞.通過分析堵塞區(qū)水的密度分布、氫鍵壽命、水分子的結(jié)構(gòu)特征,揭示了反常堵塞是由腔體中振動膜的高頻振動引發(fā)水的結(jié)構(gòu)相變造成的.

1 引言

無閥微泵是一種小型化、高精度和可控的微流控器件,由于其在化學(xué)、生物醫(yī)藥、信息技術(shù)、能源等眾多領(lǐng)域具有廣泛的應(yīng)用前景而受到越來越多的關(guān)注[1-5].近年來,納米科學(xué)和納米技術(shù)的迅猛發(fā)展不斷將流控器件的尺寸推向納米量級.2001 年,Hummer等[6]率先從理論上預(yù)言了水分子能夠快速通透剛好能容納一個水分子的(6,6)碳納米管.Ghosh等[7]在2003 年從實驗上證明了水可以在納米尺度下的碳管中流動.從此開啟了對碳納米管中受限水性質(zhì)及其應(yīng)用的研究.系列研究表明,受限在納米孔道里的流體呈現(xiàn)出許多不同于宏觀狀態(tài)下的性質(zhì)[8-12].Holt等[9]發(fā)現(xiàn),通過直徑為1—2 nm的碳納米管中的水流量比連續(xù)流體力學(xué)模型計算的值高3 個數(shù)量級以上,進一步研究還發(fā)現(xiàn)當(dāng)碳管直徑小于一定值后,碳管內(nèi)水的摩擦力會消失[13].利用水在碳納米管中幾乎無摩擦地流動的特性,各種基于碳納米管的無閥納米泵的設(shè)計被提出[14-21].水流在納米泵中表現(xiàn)出豐富的物理現(xiàn)象[22-28].然而,在納米尺度下,宏觀流體力學(xué)規(guī)律的適用性將受到挑戰(zhàn),尺度效應(yīng)變得十分明顯,因此納米流體器件中的流體行為及其背后的機理仍需進一步深入研究.

錐形通道常被作為無閥泵的重要部分,由于其不對稱的形狀,在整流振蕩流體方面起著至關(guān)重要的作用[29-31].而且,不對稱的錐形通道在生物系統(tǒng)中也很常見[32-34],主要負責(zé)控制生命活動所需的物質(zhì)跨膜輸運.由于生物通道結(jié)構(gòu)和化學(xué)組分的復(fù)雜性,其微觀機理仍未徹底弄清楚.研究在納米尺度下不對稱性孔道中水的輸運性質(zhì)對理解跨膜輸運過程及設(shè)計高效可控的納米泵都有重要的意義.

在錐形通道的孔口,由于顆粒與孔壁之間的摩擦以及顆粒內(nèi)部之間的相互作用,在一定條件下會形成阻塞拱門,顆粒物質(zhì)(如糧食、沙子或白砂糖等)會發(fā)生堵塞[35,36].外部機械振動通常是恢復(fù)流動的一種常用方法[37].而水與顆粒物質(zhì)不同,理論及實驗研究表明水分子可以快速通透納米級通道[9,38].甚至對于直徑僅為0.81 nm 的(6,6)碳納米管,水分子也能快速通透[6].本文利用分子動力學(xué)模擬發(fā)現(xiàn)在無閥納米泵的錐形碳納米管中出現(xiàn)了反常的水分子堵塞現(xiàn)象.值得注意的是,納米泵中所采用的開口錐形碳納米管最窄區(qū)域的直徑為1.2 nm,大于(6,6)碳納米管的直徑,更遠大于一個水分子的直徑(約0.28 nm).在正常情況下,水分子應(yīng)該可以快速通過這樣的納米孔道.模擬結(jié)果分析表明,在外部高頻振動驅(qū)動下,水在納米泵窄口處發(fā)生了無序到有序的轉(zhuǎn)變.其主要原因在于振動膜的高頻振動大幅延長了錐形碳納米管窄口區(qū)域內(nèi)水分子間氫鍵的平均壽命,促使受限在該區(qū)域的水分子形成了更加有序的氫鍵網(wǎng)絡(luò).

2 模擬方法

2.1 納米泵模型與參數(shù)

模擬系統(tǒng)如圖1 所示,無閥納米泵由左右兩個錐形碳納米管和中間的泵體腔構(gòu)成.中間的泵體腔,由上下2 片(5 nm × 5 nm)單層石墨烯及左右2 片相同尺寸的碳平面組成,x方向上采用周期性邊界條件.其中,底部的石墨烯固定,而頂部的石墨烯作為振動膜,除了邊緣一圈碳原子固定外其余的碳原子并不固定.左邊的碳平面中間開有直徑為1.5 nm 的孔,與左邊的錐形碳納米管相連.右邊的碳平面中間開有直徑為3.5 nm 的孔,與右側(cè)的錐形碳納米管相連.作為通道的兩個錐形碳納米管尺寸相同,左端管口直徑為3.2 nm,右端管口直徑為1.2 nm,管長為3.2 nm.在納米泵的兩端通過2 個中間開口的碳平面分別連接左右兩側(cè)的水浴.為了納米泵的穩(wěn)固,除了泵體頂部的振動膜,其余組件在模擬過程中被固定.整個體系里的水分子數(shù)為9935 個.位于泵體頂部石墨烯(振動膜)中心的標(biāo)記碳原子(橘色小球)為振動原子,用來提供驅(qū)動力,振動原子沿著y軸運動,其運動方程為y=y0+Acos(ωt+φ),其中y0是振動碳原子的初始位置,A是振幅,ω是角頻率,φ是初始相位.頻率f為ω/(2π) .因此,當(dāng)振動原子周期性振蕩時,其所在的石墨烯平面上的其他碳原子由于共價鍵的作用會協(xié)同振動.為了研究小幅微擾產(chǎn)生的影響,而非力學(xué)形變產(chǎn)生的排空效應(yīng),本文中選取振幅A為0.18 nm[11].黃色標(biāo)記的部分代表左右兩個錐形納米管窄口端厚度為0.4 nm 的區(qū)域,分別記為LP (left pore)和RP (right pore).

圖1 模擬系統(tǒng)圖,由左右兩個錐形碳納米管及中間泵體組成的無閥納米泵,通過兩個碳平面(灰色小球)連接左右水浴,紅色小球和白色小球分別表示水分子中的氧原子和氫原子;橘色小球代表振動原子;左右錐形碳納米管的窄口端用黃色標(biāo)記,分別記為LP (left pore)和RP (right pore)Fig.1.Simulation system.A valve-less nanopump consists of carbon nanocones on both sides,a fluid cavity in the center,connected by carbon planes (gray spheres) to two reservoirs.Red and white spheres represent oxygen and hydrogen atoms in water molecules,respectively.The orange sphere represents the vibrating atom.The narrow ends of the carbon nanocones are marked in yellow and labeled as LP (left pore) and RP (right pore),respectively.

對于上述系統(tǒng),利用分子動力學(xué)軟件 Gromacs4.0.7[39],在正則系綜下進行模擬.模擬盒子大小為5.1 nm × 5.0 nm × 18.2 nm.采用 V-rescale方法使系統(tǒng)溫度保持在 300 K,溫度耦合系數(shù)τ=0.1 ps[40].系統(tǒng)中的水分子采用TIP4P 模型.本文中石墨烯均不帶電,石墨烯平面的總勢能包括碳原子間的范德瓦耳斯相互作用及成鍵相互作用,可以用以下勢能函數(shù)來表示:其中所有參數(shù)來自文獻[7],σCC=0.34 nm,εCC=0.3612 kJ/mol,鍵長r0=0.142 nm,鍵角θ0=120°,勁度系數(shù)分別為kb=393 960 kJ/(mol·nm2),kθ=527 kJ/(mol·(°)2)和kξ=52.718 kJ/(mol·(°)2) .碳與水分子的范德瓦耳斯相互作用參數(shù)σCO=0.3275 nm,εCO=0.4802 kJ/mol.截斷半徑選取為1.4 nm.對于水-水之間的長程靜電相互作用,利用PME (particle-mesh ewald)方法進行統(tǒng)計,截斷半徑為 1.4 nm.系統(tǒng)在3 個方向均采用周期性邊界條件.運動方程采用Leapfrog 方法進行積分,時間步長為 1.0 fs.

2.2 流量的統(tǒng)計

本文設(shè)計了一個由2個開口的碳納米錐(錐形碳納米管)和1個泵體構(gòu)成的納米泵(圖1),泵體中膜的振動使系統(tǒng)遠離平衡態(tài),而左右2個不對稱的錐形碳納米管則破壞了空間對稱性,根據(jù)棘輪理論,這個系統(tǒng)可以形成水的定向流.定向流的大小由平均凈水通量來描述,其定義為

其中,t是總模擬時間,Nright(t) 指在t時間內(nèi)從納米泵的左端進入后從右端出去的水分子總數(shù),而Nleft(t)則表示從納米泵的右端進入后從左端出去的水分子總數(shù).本文選取頻率范圍0—1667 GHz內(nèi)的24 個頻率進行模擬研究,每個頻率對應(yīng)的模擬時間不少于20 ns,最后5 ns 的數(shù)據(jù)用來分析統(tǒng)計.

2.3 氫鍵的平均壽命

水分子之間的氫鍵定義本文采用幾何方法: 當(dāng)兩個水分子的氧原子間距小于3.5 ?,同時O—H···O 角小于30°時,則形成一個氫鍵.氫鍵的平均壽命由氫鍵維持時間的概率分布積分得到,設(shè)壽命為t的氫鍵的概率分布函數(shù)為P(t),氫鍵平均壽命計算公式如下:

且P(t)滿足歸一化條件,=1 .具體地,在模擬過程中,每一個狀態(tài)氫鍵的統(tǒng)計時間為5 ns,保存步長為1 fs,對于每個時間步,檢查系統(tǒng)中的所有可能氫鍵,記錄氫鍵形成和斷裂的事件,計算出它們的壽命,通過直方圖進行統(tǒng)計分析,最后根據(jù)壽命分布通過加權(quán)平均算出氫鍵的平均壽命.

3 結(jié)果與討論

本文設(shè)計的納米泵可利用外部振動實現(xiàn)定向流.圖2 顯示了納米泵中的水流量隨振動頻率變化的規(guī)律.顯然,當(dāng)系統(tǒng)處于平衡態(tài)時,納米泵中的水流量為0.在低頻范圍(0—300 GHz)內(nèi),水通量在零附近小幅波動.當(dāng)振動頻率繼續(xù)增大,水分子定向流動明顯增強,在625 GHz 時,納米泵中流量達到峰值,約為83 ns-1,相當(dāng)于生物水通道中水流量(3.9 ± 0.6) ns-1的20倍[41].有趣的是再進一步提高振動頻率時,水流量反而迅速減小,并在1000 GHz 時減小至0,而且,在振動頻率大于1000 GHz 的范圍內(nèi)流量始終為零.通過分析模擬數(shù)據(jù),發(fā)現(xiàn)從左側(cè)(右側(cè))水浴通過納米泵到達右側(cè)(左側(cè))水浴的水分子數(shù)Nright(t)和Nleft(t) 竟然為0,這意味著納米水泵中的錐形納米管被水堵塞了!值得注意的是,本研究采用的錐形納米管小口直徑為1.2 nm,遠大于一個水分子的尺度(~ 0.28 nm),比水分子可以快速通透的(6,6) 碳納米管的直徑也要大0.4 nm 左右,但是水流卻被堵住.通常,在宏觀世界里,當(dāng)顆粒物質(zhì)在通過孔口被堵住時,通常采用機械振動來恢復(fù)流動狀態(tài)[37].但是,在本文中無閥納米泵中的水分子卻因為高頻振動而導(dǎo)致從流動狀態(tài)轉(zhuǎn)變?yōu)槎氯麪顟B(tài).

圖2 水流量隨振動頻率的變化Fig.2.Average water flux as a function of the vibration frequency.

為了更深入地理解納米泵從流動到堵塞的轉(zhuǎn)變機制,在最易引起堵塞的錐形納米管的窄口區(qū)域里研究了水分子概率密度分布在不同頻率下呈現(xiàn)的特征.所選窄口區(qū)域的厚度為0.4 nm,為了討論方便,將左右兩個納米管的窄口區(qū)分別標(biāo)記為左孔LP 和右孔RP,如圖1 所示.分析結(jié)果顯示水分子概率密度分布特征在LP 和RP 兩個區(qū)域類似.圖3 展示了水分子在左孔LP 區(qū)域的概率密度分布,其中藍色代表低密度,紅色代表高密度.圖3(a)顯示,在f=384 GHz 的情況下,水的密度分布幾乎是按角度均勻分布的,這表明水分子可以自由地通過錐形碳納米管的窄口,即水處于流動狀態(tài).在f=625 GHz的情況下也可以看到類似的水分子密度分布,見圖3(b).當(dāng)振動頻率進一步升高至833 GHz時,水分子的概率密度分布發(fā)生了顯著的變化.由圖3(c)可觀察到水分子的密度分布不再均勻而是出現(xiàn)4 個高密度區(qū)域,這種非均勻的分布特征表明在LP 區(qū)中的水分子出現(xiàn)了有序的結(jié)構(gòu).而這種有序結(jié)構(gòu)使得納米泵中的水流量迅速減小(見圖2).當(dāng)f=1250 GHz 時,水分子的分布仍然呈現(xiàn)出非均勻性,形成5 個較規(guī)則的高密度區(qū)域,相對于f=833 GHz時的形態(tài),其高密度區(qū)域的數(shù)量增加.相應(yīng)地,水流量為零,通道已完全被堵塞.

圖3 在錐形納米管窄口端0.4 nm 厚的區(qū)域內(nèi)水分子的概率密度分布圖(a) f=384 GHz;(b) f=625 GHz;(c) f=833 GHz;(d) f=1250 GHzFig.3.Probability density distributions of water molecules in LP for different vibrational frequency: (a) f=384 GHz;(b) f=625 GHz;(c) f=833 GHz;(d) f=1250 GHz.

由上述分析知,當(dāng)振動頻率達到一定的值后,水分子在錐形碳納米管窄口區(qū)域的分布會出現(xiàn)局部高密度,而且高密度區(qū)域的數(shù)目會隨頻率的升高而增加.為了對水分子的分布有更清晰的認(rèn)識,進一步分析了LP 和RP 區(qū)域中的平均水分子個數(shù)隨頻率的變化.每個分析區(qū)域的厚度為0.4 nm,大約可容納2 層水分子.為了便于討論,稱右側(cè)的一層水為第1 水分子層(記為L1和R1),左側(cè)相鄰的另一層水為第2 水分子層(記為L2和R2).如圖4所示,在L1中,當(dāng)振動頻率不超過625 GHz 時,平均水分子數(shù)幾乎保持不變,〈N〉約為4.4 個.當(dāng)振動頻率從625 GHz 增加到1250 GHz 時,〈N〉 從4.4 個增加到5.4 個,進一步升高振動頻率時,該水層的分子數(shù)不再增加.第2 水層L2的平均分子數(shù)隨頻率變化也經(jīng)歷一個類似的過程,從低頻到高頻水分子數(shù)增幅約為1.5 個.LP 區(qū)域右邊的RP 里,第1 水層和第2 水層的平均水分子數(shù)同樣也出現(xiàn)了階梯型變化,相應(yīng)的水分子數(shù)的增幅分別為1 個和2 個.第2 層水的平均水分子數(shù)增幅比第1 層水大是由于第2 層水所處的納米管的橫截面積更大所致.可見,錐形碳納米管窄口區(qū)域的水分子數(shù)增加引起了該區(qū)域里水分子概率密度的空間分布的突變,即從無序的均勻分布到有序的規(guī)則分布.對比圖2,LP 和RP 區(qū)域中的平均水分子個數(shù)隨頻率增加的區(qū)間(625—1250 GHz)與納米泵中水流開始減小并完全堵塞的區(qū)間是一致的.

圖4 LP 和RP 區(qū)域的平均水分子數(shù),兩個區(qū)域厚度為0.4 nm,分別包含2 層水分子,L1 和R1 分別表示兩個區(qū)域中右側(cè)的(管徑最小)水分子層;L2 和R2 則表示左側(cè)的水分子層Fig.4.Average number of water molecules in each water layer in LP and RP.L1 and R1 represent the right-side (minimum diameter) water molecule layers in the two regions,while L2 and R2 represent the left-side water molecule layers.

本文計算了LP 區(qū)域內(nèi)氫鍵的平均壽命.同時還選取了泵體腔內(nèi)氫鍵的平均壽命作為比較.氫鍵壽命的定義和詳細計算方法見模擬方法部分的介紹.圖5(a)顯示了LP 區(qū)域內(nèi)水分子之間形成氫鍵的平均壽命隨頻率的變化規(guī)律.中間腔室里氫鍵的平均壽命在0.3 ps 附近小幅波動.在f <600 GHz的范圍內(nèi),LP 區(qū)域里氫鍵的平均壽命保持在0.27 ps 左右,而在這個頻率范圍,納米管中的水流量隨頻率的增大明顯增大(見圖2),可見管中水的流速對氫鍵壽命并無明顯的影響,而此時,水分子空間概率分布處于一種無序的結(jié)構(gòu),如圖5(a)插圖I 所示.當(dāng)振動頻率從600 GHz 增至1000 GHz,氫鍵平均壽命急劇延長.在f=625 GHz 的情況下,氫鍵平均壽命約為0.4 ps.當(dāng)頻率為833 GHz 時,氫鍵平均壽命增加到約0.51 ps.氫鍵壽命的延長表明水分子之間形成更穩(wěn)定的氫鍵網(wǎng)絡(luò).當(dāng)頻率大于1000 GHz 時,氫鍵壽命在0.68 ps 附近小幅變化,意味著氫鍵網(wǎng)絡(luò)基本達到穩(wěn)定狀態(tài),水分子形成了有序的結(jié)構(gòu).圖5(a)插圖II 顯示了該頻率范圍內(nèi)的一個典型的結(jié)構(gòu)快照圖,可以清晰地看到水分子在管口LP 區(qū)域形成了一個有序的五邊形結(jié)構(gòu).為了進一步研究LP 區(qū)域水分子的有序結(jié)構(gòu),詳細統(tǒng)計了水分子方位角的概率分布.水分子方位角的定義為,以管口橫截面的中心為參考點,氧原子與參考點的連線與水平軸x軸的夾角(見圖1).概率分布如圖5(b)所示.當(dāng)頻率為f=1250 GHz時,分布圖中出現(xiàn)了5 個峰值,其中角度位于第一、二、三象限中的4 個峰值非常尖銳,說明對應(yīng)的水分子的空間方位非常穩(wěn)定;在第四象限的峰值相對平緩,說明對應(yīng)的水分子在約20°角度范圍內(nèi)運動.相比于頻率250 GHz 時的分布情況,水分子的分布有序度大大增強.與圖3 中水分子概率密度分布從空間均勻到非均勻分布的轉(zhuǎn)變相一致.

圖5 (a)氫鍵的平均壽命隨頻率的變化,插圖分別為f=250 GHz和f=1250 GHz 對應(yīng)的 LP 里第1 層水的模擬快照(青色小球表示碳原子,紅色和白色小球分別表示水分子中的氧原子和氫原子);(b) LP 區(qū)域第1 層水分子在橫截面上的方位角概率分布圖Fig.5.(a) Lifetime of hydrogen bonds as a function of the vibration frequency.Two typical snapshots are shown in insets.The cyan spheres represent carbon atoms,with red and white spheres representing oxygen and hydrogen atoms in water molecules.One is for the low frequency (f=250 GHz) and the other is for the high frequency (f=1250 GHz).(b) Probability distribution of the azimuthal angle of water molecules in the first layer in the LP region on the cross-section for f=250 GHz and f=1250 GHz.

4 結(jié)論

本文通過分子動力學(xué)模擬,發(fā)現(xiàn)在無閥納米泵的錐形孔道中水流堵塞的現(xiàn)象,盡管孔道相對于水分子的尺寸足夠大且并無其他雜質(zhì)分子.這易讓人聯(lián)想到宏觀世界里常見的漏斗中顆粒物質(zhì)堵塞現(xiàn)象.然而,與漏斗中顆粒物堵塞不同的是,具有高頻振動的外部擾動并不能恢復(fù)水流.相反,在高頻振動的幫助下,錐形通道孔口處的水表現(xiàn)出從無序到有序的結(jié)構(gòu)相變.引起這種結(jié)構(gòu)相變的原因是,在高頻振動的推動下水分子快速擠入錐形孔道里,水分子間形成了緊密的氫鍵網(wǎng)絡(luò),而且,隨著振動頻率的增大,氫鍵平均壽命顯著延長,最終在孔口處形成了有序的多邊形結(jié)構(gòu).錐形通道常被作為無閥泵的重要部分,由于其不對稱的形狀,在整流振蕩流體方面起著至關(guān)重要的作用.這項工作系統(tǒng)研究了振動頻率對基于錐形納米管構(gòu)成的無閥納米泵中水輸運的影響,發(fā)現(xiàn)了納米級無閥泵中水流堵塞的新現(xiàn)象,并揭示了堵塞背后的物理機制.

數(shù)據(jù)可用性聲明

本篇論文的關(guān)聯(lián)數(shù)據(jù)可在科學(xué)數(shù)據(jù)銀行https://www.doi.org/10.57760/sciencedb.j00213.00014 中訪問獲取.

猜你喜歡
振動區(qū)域
振動的思考
噴水推進高速艇尾部振動響應(yīng)分析
永久基本農(nóng)田集中區(qū)域“禁廢”
分割區(qū)域
This “Singing Highway”plays music
振動攪拌 震動創(chuàng)新
中國公路(2017年18期)2018-01-23 03:00:38
中立型Emden-Fowler微分方程的振動性
關(guān)于四色猜想
分區(qū)域
基于嚴(yán)重區(qū)域的多PCC點暫降頻次估計
電測與儀表(2015年5期)2015-04-09 11:30:52
主站蜘蛛池模板: 亚洲国产成人综合精品2020| 一区二区影院| 中文字幕无码电影| 91po国产在线精品免费观看| 亚洲国模精品一区| 国产a网站| 中文字幕久久亚洲一区| 亚洲国产天堂在线观看| 国产69精品久久久久孕妇大杂乱 | 国产jizzjizz视频| 国产在线视频导航| 男人的天堂久久精品激情| 九九九九热精品视频| 欧美亚洲一区二区三区在线| 亚洲国产欧美国产综合久久| 国产成人精品免费av| 亚洲第一视频网| 欧美a在线看| 成年网址网站在线观看| 国产精品林美惠子在线观看| 免费一级无码在线网站| 黄色一级视频欧美| 国产欧美亚洲精品第3页在线| 欧美性天天| 国产av一码二码三码无码| 成年免费在线观看| 小说区 亚洲 自拍 另类| 在线观看精品国产入口| 亚洲大尺码专区影院| 国产欧美又粗又猛又爽老| 好吊妞欧美视频免费| 国产福利影院在线观看| 久99久热只有精品国产15| 亚洲性一区| 精品久久高清| 久久公开视频| 亚洲妓女综合网995久久| 亚洲无码在线午夜电影| 欧美日韩成人| 全免费a级毛片免费看不卡| 久久亚洲综合伊人| 91麻豆精品国产高清在线| 国产91丝袜| 欧洲在线免费视频| 欧美成人午夜影院| 伊人色天堂| 日韩精品毛片人妻AV不卡| 亚洲欧美激情另类| 国产在线视频自拍| 亚洲色婷婷一区二区| 精品久久国产综合精麻豆| 国产又大又粗又猛又爽的视频| 日韩 欧美 小说 综合网 另类| 国产超碰在线观看| 国产香蕉国产精品偷在线观看| 蜜桃视频一区二区三区| 亚洲高清日韩heyzo| 丝袜久久剧情精品国产| 欧美精品一区二区三区中文字幕| 思思热精品在线8| 国产精品永久在线| 毛片在线看网站| 国产高潮流白浆视频| 永久毛片在线播| 欧美a在线视频| 丰满人妻被猛烈进入无码| 国产拍揄自揄精品视频网站| 香蕉国产精品视频| 欧美午夜性视频| 亚洲a级毛片| 高潮毛片免费观看| www中文字幕在线观看| 狠狠躁天天躁夜夜躁婷婷| 亚洲日韩欧美在线观看| 国产精品人莉莉成在线播放| 亚洲午夜综合网| 91po国产在线精品免费观看| 国产办公室秘书无码精品| 亚洲第一色网站| 国产女人在线视频| 国产www网站| 亚洲国产系列|