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

基于多物理場(chǎng)耦合的毛細(xì)水高度研究

2021-06-17 11:14:22鄧改革何建國(guó)康寧波
水土保持研究 2021年4期
關(guān)鍵詞:模型研究

鄧改革, 何建國(guó), 康寧波

(1.寧夏大學(xué) 土木與水利工程學(xué)院, 銀川 750021; 2.中國(guó)礦業(yè)大學(xué)銀川學(xué)院 信息工程學(xué)院, 銀川 750021; 3.寧夏大學(xué) 農(nóng)學(xué)院, 銀川 750021)

毛細(xì)作用是指浸潤(rùn)液體在細(xì)管里升高和不浸潤(rùn)液體在細(xì)管里降低的現(xiàn)象,是液體的慣性力、黏性力、毛細(xì)力與重力等共同作用的結(jié)果[1]。土壤毛細(xì)水是靠毛細(xì)吸引力而保持于土壤毛細(xì)孔隙中的水。土壤毛細(xì)水上升高度對(duì)包括路基凍壞以及農(nóng)田土壤次生鹽漬化等在內(nèi)的眾多過(guò)程具有重要的影響,目前國(guó)內(nèi)外學(xué)者在毛細(xì)水上升高度方面已經(jīng)開(kāi)展了大量的研究工作但尚不深入。

于丹等[2]進(jìn)行了不同干密度重塑黃土的毛細(xì)上升速率和最大高度的研究,發(fā)現(xiàn)毛細(xì)最大上升高度與干密度之間非線性關(guān)系,而是呈現(xiàn)近似拋物線的函數(shù)關(guān)系。李先瑞等[3]針對(duì)改性黃土進(jìn)行了垂直方向毛細(xì)水上升作用的研究,發(fā)現(xiàn)石灰能夠有效降低土體內(nèi)的含水率,且隨著含量的增加,吸水作用越明顯,而水泥對(duì)于土體內(nèi)部結(jié)構(gòu)的改性作用更大,提升土體強(qiáng)度和遇水穩(wěn)定性,阻礙毛細(xì)水上升作用顯著。肖紅宇等[4]進(jìn)行了基于黏性土分形特征的毛細(xì)水上升高度研究,推導(dǎo)出土體內(nèi)毛細(xì)水上升高度預(yù)測(cè)公式,并將公式預(yù)測(cè)結(jié)果與實(shí)測(cè)結(jié)果進(jìn)行了對(duì)比分析,發(fā)現(xiàn)誤差在工程允許范圍之內(nèi),具有一定的實(shí)際應(yīng)用價(jià)值。落宇杰等[5]進(jìn)行的壓實(shí)黃土狀粉土毛細(xì)特性試驗(yàn)發(fā)現(xiàn),毛細(xì)水上升的速率隨毛細(xì)上升高度的增加而逐漸衰減。袁玉卿等[6]通過(guò)室內(nèi)試驗(yàn)研究了豫東黃泛區(qū)粉砂土毛細(xì)水上升規(guī)律及控制技術(shù),研究發(fā)現(xiàn)隨時(shí)間的延長(zhǎng)毛細(xì)水上升高度逐漸穩(wěn)定,毛細(xì)水上升速度與壓實(shí)度呈反比,級(jí)配碎石、水泥穩(wěn)定土、纖維水泥穩(wěn)定土能有效阻隔毛細(xì)水上升。

胡明鑒等[1]進(jìn)行了砂性土條件下,不同粒徑及潛水礦化度組合的毛細(xì)水上升試驗(yàn),研究了不同粒徑及潛水礦化度下對(duì)毛細(xì)水上升高度及上升速度的規(guī)律。王聰[7]、栗現(xiàn)文[8]等通過(guò)觀測(cè)毛細(xì)水上升速度和高度,研究不同濃度鹽溶液和鹽漬土對(duì)毛細(xì)水上升的影響,并對(duì)這些影響因素進(jìn)行機(jī)理分析。

魏樣[9]、童玲[10]等分別研究了石油及柴油污染對(duì)土壤毛細(xì)水上升特性的影響。章求才等[11]研究了溫度和氣壓對(duì)某金屬礦山尾礦壩中毛細(xì)水上升規(guī)律的影響,發(fā)現(xiàn)溫度對(duì)毛細(xì)水上升的影響較氣壓大;在毛細(xì)水上升初期,溫度、氣壓對(duì)毛細(xì)水上升規(guī)律的影響不大,但在上升后期其影響逐漸明顯。

Masoodi等[12]采用達(dá)西定律模擬了液體向由纖維素和高吸水性纖維制成的紙狀膨脹多孔介質(zhì)中的芯吸現(xiàn)象。提出了新的毛細(xì)模型,將多孔介質(zhì)中達(dá)西定律與質(zhì)量守恒方程相結(jié)合,通過(guò)向質(zhì)量守恒方程添加附加的匯或源項(xiàng),以說(shuō)明基質(zhì)的膨脹和液體吸收。新模型預(yù)測(cè)芯吸速率較修正的沃什伯恩方程預(yù)測(cè)芯吸速率更好。

Aghajani等[13]提出了一種改善土壤中毛細(xì)水上升高度準(zhǔn)確性的方法,該方法重新考慮了傳統(tǒng)方法中被忽略的毛細(xì)區(qū)域中水分和吸力的變化,試驗(yàn)結(jié)果表明,改進(jìn)的方法比以前的毛細(xì)管上升解決方案更加準(zhǔn)確和通用。

Fili等[14]進(jìn)行了采用多相格子玻爾茲曼方法(LBM)檢驗(yàn)土壤持水曲線中毛細(xì)現(xiàn)象的研究。LBM模型首先針對(duì)基準(zhǔn)問(wèn)題進(jìn)行驗(yàn)證,然后用于模擬靜態(tài)粒子陣列,以此來(lái)研究初始流體密度分布對(duì)系統(tǒng)毛細(xì)管響應(yīng)的影響。結(jié)果表明,濕潤(rùn)鋒的統(tǒng)一性直接影響土壤骨架毛細(xì)響應(yīng)峰值的大小。

Zhou等[15]研究了非飽和土的抗剪強(qiáng)度與毛管持水曲線的關(guān)系,提出了一種新的、考慮吸附水毛細(xì)凝結(jié)的非飽和土壤保水模型。模型中,將土壤的飽和度分為基于毛細(xì)管水的飽和度和基于吸附水的飽和度。通過(guò)對(duì)部分飽和的雙圓柱系統(tǒng)的分析,提出了一種新的非飽和土抗剪強(qiáng)度準(zhǔn)則。研究結(jié)果將大大改進(jìn)對(duì)試驗(yàn)數(shù)據(jù)預(yù)測(cè)的準(zhǔn)確性。

以上對(duì)毛細(xì)水高度的研究多側(cè)重于室內(nèi)土柱試驗(yàn),研究?jī)?nèi)容多集中于不同土質(zhì)結(jié)構(gòu)、礦化水、污染物以及外界溫度壓力等方面,取得了一定的研究成果。而基于多物理場(chǎng)耦合的毛細(xì)水上升高度機(jī)理方面研究尚未見(jiàn)公開(kāi)報(bào)道。本文使用多孔介質(zhì)中多相流理論,通過(guò)多物理場(chǎng)耦合軟件COMSOL Multiphysics建立模型,耦合多孔介質(zhì)界面中達(dá)西定律和相輸運(yùn),最后通過(guò)實(shí)驗(yàn)室土柱試驗(yàn)對(duì)模擬結(jié)果進(jìn)行驗(yàn)證。研究結(jié)果將對(duì)弄清毛細(xì)水上升機(jī)理提供新的思路和方法。

1 研究方法

1.1 數(shù)值模擬研究

1.1.1 理論模型 在經(jīng)典的多孔介質(zhì)模型中,構(gòu)成其基體的尺寸及其之間的孔隙在毛細(xì)吸水過(guò)程中保持恒定。目前計(jì)算毛細(xì)水上升高度的模型主要有兩種,一種基于Lucas-Washburn方程,另一種基于達(dá)西定律[12]。這兩種方法在計(jì)算毛細(xì)上升高度方面均存在缺陷,致使理論計(jì)算值與實(shí)際結(jié)果之間偏差較大。下面介紹多物理場(chǎng)耦合多孔介質(zhì)中的多相流計(jì)算模型,該多物理場(chǎng)耦合了多孔介質(zhì)界面中的達(dá)西定律和相輸運(yùn)。多孔介質(zhì)界面中的相輸運(yùn)遵循潤(rùn)濕或非潤(rùn)濕流體體積分?jǐn)?shù)的獨(dú)立方程:

(1)

式中:εp為孔隙率;ρi為第i相的流體密度(kg/m3);?為梯度算子;si為相的體積分?jǐn)?shù);κ為滲透率(m2);κri為相對(duì)滲透率(給定流體的飽和度的函數(shù));μi為流體的動(dòng)態(tài)黏度(Pa/s);pi為壓力(Pa);?pi為壓力差(Pa);g為重力加速度(m/s2)。

由于兩相的體積分?jǐn)?shù)之和為1,因此可以從以下公式計(jì)算剩余的體積分?jǐn)?shù):

s1=1-s2

(2)

式中:s1為非潤(rùn)濕相體積分?jǐn)?shù);s2為潤(rùn)濕相體積分?jǐn)?shù)。

毛管壓力pc作為潤(rùn)濕相sw(模型中為s2)的飽和度和入口毛細(xì)管壓力pec的函數(shù)進(jìn)行計(jì)算。通過(guò)使用Brooks-Corey模型[16-18],毛管壓力由下式得出:

(3)

式中:λp為孔隙分布指數(shù);pec為入口毛細(xì)管壓力。基于Brooks-Corey模型,濕相和非濕相的相對(duì)滲透率由下式給出:

(4)

(5)

式中:κrs為濕相相對(duì)滲透率;κrsn為非濕相相對(duì)滲透率。

通過(guò)達(dá)西定律接口將達(dá)西定律與連續(xù)性方程相結(jié)合:

(6)

式中:ρ為流體密度;μ為流體動(dòng)力黏度;κ為滲透率(m2)。

1.1.2 模擬參數(shù)確定 以1,2,3號(hào)管土柱為數(shù)值模擬的研究對(duì)象,建立模型水力邊界條件:地下水位線在砂土柱的最底端并保持穩(wěn)定,砂土柱頂端設(shè)置為自由邊界[6]。由于豎管的對(duì)稱結(jié)構(gòu),同時(shí)為了減小網(wǎng)格數(shù)量以提高計(jì)算效率,建立豎管的二維模型,所建立模型寬度為41 mm,高度為1 000 mm,并對(duì)所建立模型進(jìn)行網(wǎng)格劃分,建模相關(guān)參數(shù)見(jiàn)表1—2。計(jì)算時(shí)間為1 600 min,時(shí)間步長(zhǎng)為5 min,通過(guò)COMSOL Multiphysics多物理場(chǎng)仿真軟件,利用有限差分法對(duì)上述方程進(jìn)行離散化求解。

表1 水、空氣和多孔材料的材料參數(shù)

表2 不同編號(hào)砂土多孔介質(zhì)特性參數(shù)

1.2 土柱試驗(yàn)研究

本文采用豎管法研究毛細(xì)水上升高度。試驗(yàn)用砂土編號(hào)為#1,#2,#3,按照相關(guān)規(guī)范進(jìn)行物理指標(biāo)試驗(yàn),測(cè)定相應(yīng)參數(shù),結(jié)果見(jiàn)表3。

表3 試驗(yàn)用砂土物理參數(shù)

試驗(yàn)采用自制的毛細(xì)管水上升高度裝置(圖1)。裝置由儀器支架、3支有機(jī)玻璃管、有機(jī)玻璃水箱、標(biāo)尺、微型循環(huán)水泵及相連接的橡膠管組成。有機(jī)玻璃管長(zhǎng)度為100 cm,內(nèi)徑41 mm,外徑45 mm,頂端敞口可以與大氣連通。將編號(hào)為#1,#2,#3的砂土分別裝入編號(hào)為1,2,3的有機(jī)玻璃管中。三根有豎管放置在有機(jī)玻璃水箱托盤(pán)中,托盤(pán)高出水箱液面一定距離,通過(guò)微型循環(huán)水泵將水箱中水補(bǔ)充至托盤(pán)中,使托盤(pán)中液面保持恒定。

圖1 毛細(xì)管水上升高度裝置

毛細(xì)水上升高度主要試驗(yàn)步驟如下:將篩分配制好的一定質(zhì)量的砂樣,通過(guò)漏斗將土樣裝入有機(jī)玻璃管中,使其具有相同的壓實(shí)度。打開(kāi)微型循環(huán)水泵并開(kāi)始計(jì)時(shí),觀測(cè)濕潤(rùn)鋒的位置,記錄毛細(xì)水上升高度與時(shí)間。前期毛細(xì)水上升較快,每隔10 min記錄濕潤(rùn)鋒所在位置,之后隨著毛細(xì)水上升速度變慢,讀數(shù)時(shí)間間隔相應(yīng)加大。

1.3 現(xiàn)場(chǎng)試驗(yàn)研究

在實(shí)訓(xùn)基地內(nèi)部選擇合適的場(chǎng)地,在場(chǎng)地上方提前30 d布置遮擋物防止外界雨水滲入,然后開(kāi)挖長(zhǎng)寬高分別為1 m×1 m×0.8 m的土坑,通過(guò)水泵向坑底部供水,不斷調(diào)整水泵的流量使水面保持恒定,通過(guò)鋼卷尺讀取土壤濕潤(rùn)鋒位置的變化見(jiàn)圖2。

圖2 毛細(xì)水上升高度現(xiàn)場(chǎng)試驗(yàn)

2 結(jié)果與分析

2.1 數(shù)值模擬結(jié)果及分析

從圖3可知在模擬1 600 min后1,2,3號(hào)管濕潤(rùn)鋒位置高度分別為16.8,23.1,37.5 cm,該模擬結(jié)果與后面實(shí)驗(yàn)室室內(nèi)土柱試驗(yàn)結(jié)果較為接近。由圖4可知,在模擬初期毛細(xì)水高度隨時(shí)間增加較快,在模擬后期毛細(xì)水高度增加緩慢并逐漸趨于穩(wěn)定。在毛細(xì)水上升高度不高的情況下,由于基質(zhì)吸力做功遠(yuǎn)大于水分重力勢(shì),所以毛細(xì)水上升速度較快;但隨著高度的上升,重力勢(shì)增加明顯,基質(zhì)吸力做的功相對(duì)減少,上升的速度就逐漸變慢[6]。

圖3 不同編號(hào)毛細(xì)管1 600 min毛細(xì)水高度

圖4 3號(hào)管7個(gè)不同時(shí)刻毛細(xì)水高度情況

2.2 土柱試驗(yàn)結(jié)果及分析

由圖5可知,毛細(xì)水上升試驗(yàn)開(kāi)始的前120 min毛細(xì)水上升速度較快。試驗(yàn)開(kāi)始供水的第2小時(shí),1,2,3號(hào)管毛細(xì)水分別上升至14.2,18,28 cm,這表明在相同壓實(shí)度條件下不同粒徑的土柱毛細(xì)水上升高度不同。3號(hào)管(粒徑范圍:0.15~0.28 mm)毛細(xì)水上升高度高于1號(hào)管(粒徑范圍:0.36~0.45 mm)及2號(hào)管(粒徑范圍:0.28~0.36 mm),說(shuō)明砂土粒徑越小毛細(xì)水上升的高度越大。隨后毛細(xì)水高度繼續(xù)增大,但是其上升速度逐漸減小。至510 min,1,2,3號(hào)管毛細(xì)水分別上升至14.6,19.6,34 cm。

圖5 毛細(xì)水前510 min上升高度

由圖6可知,在隨后的試驗(yàn)時(shí)間內(nèi),試驗(yàn)土柱內(nèi)毛細(xì)水高度持續(xù)增加,3號(hào)管毛細(xì)水高度最大,2號(hào)管次之,1號(hào)管最小。至1 600 min 1,2,3號(hào)管毛細(xì)水分別上升至15.6,22.0,37.5 cm。毛細(xì)高度隨時(shí)間變化曲線斜率代表毛細(xì)水上升速度(圖7),前120 min曲線斜率較大,表明毛細(xì)水上升速度較大,隨后曲線斜率逐漸減小,毛細(xì)水上升速度也隨之逐漸減小但并未停止,毛細(xì)水高度增加緩慢并逐漸趨于穩(wěn)定。

圖6 毛細(xì)水1 600 min上升高度

圖7 毛細(xì)水上升速度隨時(shí)間變化(0~1 600 min)

利用Origin軟件,對(duì)圖7數(shù)據(jù)的曲線進(jìn)行函數(shù)擬合,所得關(guān)系式見(jiàn)表4。上升高度關(guān)系式兩邊對(duì)時(shí)間求導(dǎo)即可得毛細(xì)水上升速度表達(dá)式(表4)。

表4 毛細(xì)水上升高度擬合關(guān)系式

由表4可知,毛細(xì)水上升高度符合冪函數(shù)性質(zhì),為增函數(shù),其變化特征是初期高度增加快,到一定程度后增加緩慢。而毛細(xì)水上升速度屬于指數(shù)減函數(shù),毛細(xì)水上升速度隨時(shí)間減小。

將數(shù)值模擬結(jié)果與室內(nèi)試驗(yàn)結(jié)果進(jìn)行對(duì)比(表5),在誤差范圍內(nèi)數(shù)值模擬與室內(nèi)試驗(yàn)結(jié)果吻合較好。

表5 模擬高度與試驗(yàn)高度對(duì)比

2.3 現(xiàn)場(chǎng)試驗(yàn)結(jié)果及分析

由于現(xiàn)場(chǎng)土坑中砂土參數(shù)與3號(hào)管中砂土參數(shù)較為接近,因此圖8試驗(yàn)數(shù)據(jù)與3號(hào)管試驗(yàn)數(shù)據(jù)較為接近。通過(guò)分析數(shù)據(jù)可知毛細(xì)高度在30 cm以內(nèi),現(xiàn)場(chǎng)試驗(yàn)毛細(xì)高度小于實(shí)驗(yàn)室內(nèi)土柱試驗(yàn)高度,其試驗(yàn)誤差小于10%,當(dāng)毛細(xì)高度超過(guò)30 cm后,試驗(yàn)誤差逐漸增大至15%左右。現(xiàn)場(chǎng)觀測(cè)發(fā)現(xiàn),土坑中30 cm以上高度區(qū)域砂土粒徑增大,且砂土中石塊等雜物增多,對(duì)毛細(xì)水上升高度產(chǎn)生了影響。將試驗(yàn)誤差考慮在內(nèi),可以認(rèn)為現(xiàn)場(chǎng)毛細(xì)水上升規(guī)律與室內(nèi)土柱試驗(yàn)以及多物理場(chǎng)仿真結(jié)果相吻合。

圖8 毛細(xì)水1 600 min上升高度現(xiàn)場(chǎng)試驗(yàn)

3 討論與結(jié)論

相同壓實(shí)度不同粒徑砂土柱內(nèi),毛細(xì)水上升高度與時(shí)長(zhǎng)近似呈冪函數(shù)關(guān)系,為增函數(shù),其變化特征是初期高度增加快,到一定程度后增加緩慢;毛細(xì)水上升速度與時(shí)長(zhǎng)近似呈指數(shù)函數(shù)關(guān)系,為減函數(shù),毛細(xì)水上升速度表現(xiàn)為初始階段較大但隨后迅速減小以致各組間差異不明顯[8]。

在本試驗(yàn)中,相同壓實(shí)度不同粒徑砂土柱毛細(xì)水上升高度由大到小依次為:3號(hào)管(粒徑范圍0.36~0.45 mm),2號(hào)管(粒徑范圍0.28~0.36 mm),1號(hào)管(粒徑范圍0.15~0.28 mm),毛細(xì)水上升高度與砂土粒徑呈負(fù)相關(guān),總體表現(xiàn)為粒徑越小,毛細(xì)水上升高度越大。在較小粒徑條件下,黏性力的作用較強(qiáng),在小粒徑范圍內(nèi)起主要作用。相反,在大粒徑范圍下,重力作用表現(xiàn)更加明顯,尤其是到毛細(xì)運(yùn)移的后期,重力的影響增大,作用不可忽略[19]。砂柱粒徑越小其內(nèi)部毛細(xì)孔隙越小,毛細(xì)管力越大,毛細(xì)潤(rùn)濕峰上移迅速,同樣時(shí)間間隔內(nèi),上升距離較高,且在試驗(yàn)結(jié)束時(shí)的高度也遠(yuǎn)遠(yuǎn)高于大粒徑的高度[19]。這一結(jié)論和胡明鑒[1]、董斌[19]等的研究結(jié)論具有相同之處。

本文通過(guò)耦合多物理場(chǎng)多孔介質(zhì)中達(dá)西定律和相輸運(yùn),較好地模擬了毛細(xì)水上升高度,并且在試驗(yàn)誤差范圍內(nèi)理論值與室內(nèi)試驗(yàn)結(jié)果吻合較好,該研究方法在其他參考文獻(xiàn)中尚未見(jiàn)公開(kāi)報(bào)道。

壓實(shí)度作為土柱試驗(yàn)的參數(shù)之一對(duì)毛細(xì)水高度具有一定的影響,本文僅研究了相同壓實(shí)度條件下不同粒徑的毛細(xì)水高度,而基于多物理場(chǎng)耦合條件下,不同壓實(shí)度對(duì)毛細(xì)水高度的影響尚有待今后進(jìn)一步研究。

試驗(yàn)用砂土作為一種多孔介質(zhì),其內(nèi)部毛細(xì)流動(dòng)現(xiàn)象是由慣性力、黏性力、毛細(xì)力與重力共同作用產(chǎn)生的,但在不同毛細(xì)高度各項(xiàng)作用力對(duì)毛細(xì)上升的貢獻(xiàn)并不相同,這也是進(jìn)一步理解毛細(xì)流動(dòng)現(xiàn)象機(jī)理的關(guān)鍵所在[19]。

由于現(xiàn)場(chǎng)砂土粒徑及土中雜物的影響使得現(xiàn)場(chǎng)試驗(yàn)與實(shí)驗(yàn)室試驗(yàn)及數(shù)值模擬結(jié)果存在一定的誤差,但是將試驗(yàn)誤差考慮在內(nèi),可以認(rèn)為現(xiàn)場(chǎng)毛細(xì)水上升規(guī)律與室內(nèi)土柱試驗(yàn)以及多物理場(chǎng)仿真結(jié)果相吻合。該結(jié)論應(yīng)用于工程領(lǐng)域時(shí),應(yīng)該充分考慮現(xiàn)場(chǎng)因素對(duì)毛細(xì)水高度的影響,以現(xiàn)場(chǎng)試驗(yàn)作為補(bǔ)充,盡量減小試驗(yàn)誤差。

猜你喜歡
模型研究
一半模型
FMS與YBT相關(guān)性的實(shí)證研究
2020年國(guó)內(nèi)翻譯研究述評(píng)
遼代千人邑研究述論
重要模型『一線三等角』
重尾非線性自回歸模型自加權(quán)M-估計(jì)的漸近分布
視錯(cuò)覺(jué)在平面設(shè)計(jì)中的應(yīng)用與研究
科技傳播(2019年22期)2020-01-14 03:06:54
EMA伺服控制系統(tǒng)研究
新版C-NCAP側(cè)面碰撞假人損傷研究
3D打印中的模型分割與打包
主站蜘蛛池模板: 国产91熟女高潮一区二区| 91探花在线观看国产最新| 国产va在线| 全色黄大色大片免费久久老太| 国产丝袜91| 91av成人日本不卡三区| 欧美国产在线看| www.日韩三级| 国产成人亚洲日韩欧美电影| 99精品视频九九精品| 亚洲国产中文欧美在线人成大黄瓜 | 98精品全国免费观看视频| 亚国产欧美在线人成| 视频在线观看一区二区| 天堂成人av| a欧美在线| 国产成人凹凸视频在线| 蝌蚪国产精品视频第一页| 国产网站免费看| 高潮爽到爆的喷水女主播视频| 精品人妻AV区| 亚洲有码在线播放| 国产一级无码不卡视频| 无码内射中文字幕岛国片| 免费一级成人毛片| 免费看a级毛片| 日韩毛片免费观看| 久久国产精品娇妻素人| 国产在线观看第二页| 2021国产乱人伦在线播放| 亚洲国产成人自拍| 亚洲精品日产精品乱码不卡| 人妻少妇久久久久久97人妻| 精品久久久久久成人AV| 国产无码精品在线播放| 国产高颜值露脸在线观看| 亚洲成肉网| 国产精品99r8在线观看| 人妻精品久久无码区| 久久精品无码中文字幕| 不卡视频国产| 91成人在线观看| 亚洲国产天堂久久综合226114| 国产午夜精品鲁丝片| 欧美亚洲一二三区| 亚洲综合精品香蕉久久网| 欧美高清视频一区二区三区| 亚洲欧美日韩动漫| 一本一道波多野结衣av黑人在线| 日韩黄色在线| 四虎亚洲国产成人久久精品| 欧美日韩中文字幕在线| 真实国产乱子伦高清| 国产福利2021最新在线观看| 午夜精品国产自在| 亚洲第一页在线观看| 九九热免费在线视频| 亚洲网综合| 无码人妻热线精品视频| 欧美视频在线播放观看免费福利资源| 久久久久国产精品熟女影院| 日韩欧美国产中文| av在线手机播放| 亚国产欧美在线人成| 日韩中文无码av超清| 国产91特黄特色A级毛片| 99视频在线免费| 国产三级韩国三级理| 中文字幕欧美日韩| 99视频免费观看| 99热这里只有精品国产99| 色哟哟国产精品一区二区| 亚洲精品成人福利在线电影| 久久国产高清视频| 国产在线精彩视频二区| 婷婷伊人五月| 在线人成精品免费视频| 欧美日本激情| 久久精品中文无码资源站| 综合色在线| 成年女人a毛片免费视频| 成人福利在线看|